rbmatlab  1.16.09
 All Classes Namespaces Files Functions Variables Modules Pages
stokes_model_default.m
1 function model = stokes_model_default
2 %function model = stokes_model_default
3 %
4 % Default settings for Stokes models.
5 
6 % IM 14.06.2013
7 
8 
9 model = [];
10 model.rb_problem_type = 'stokes';
11 
12 model.has_volume_integral_rhs = 0;
13 model.has_volume_integral_matrix = 0;
14 model.has_boundary_integral_rhs = 0;
15 model.has_boundary_integral_matrix = 0;
16 
17 
18 model.infsup_method = 'rbf_interpolant'; %=scm
19 
20 model.has_reaction = 0;
21 model.has_nonlinearity = 0;
22 model.has_dirichlet_values = 0;
23 model.has_geometry_transformation = 0;
24 
25 model.mu_names = {};
26 model.mu_ranges = {};
27 model.RB_mu_list = {};
28 model.set_mu = @set_mu_default;
29 model.get_mu = @get_mu_default;
30 model.get_reynolds_number = @(params) ...
31  params.U_mean(params)*params.L(params)/params.kinematic_viscosity(params);
32 
33 model.debug = 0;
34 model.verbose = 1;
35 model.decomp_mode = 0;
36 model.disable_caching = 1;
37 model.RB_generation_mode = 'model_RB_basisgen'; % default
38 
39 model.gen_model_data = @stokes_gen_model_data;
40 model.operators = @Fem.Assembly.operator;
41 model.detailed_simulation = @stokes_detailed_simulation;
42 model.get_inner_product_matrix = @stokes_get_inner_product_matrix;
43 model.gen_detailed_data = @stokes_gen_detailed_data;
44 model.get_dofs_from_sim_data = @(sim_data) sim_data.uh.dofs;
45 model.rb_init_data_basis = @rb_empty_init_data_basis;
46 model.error_algorithm = @stokes_relative_error_algorithm;
47 model.get_estimator_from_sim_data = @(sim_data) sim_data.Delta;
48 model.get_rb_from_detailed_data = @(detailed_data) copy(detailed_data.RB);
49 model.get_rb_size = @(model, detailed_data) size(detailed_data.RB, 2);
50 model.gen_reduced_data = @stokes_gen_reduced_data;
51 %model.gen_reduced_data = @stokes_enrich_reduced_data;
52 model.rb_simulation = @stokes_rb_simulation;
53 model.infsup_constant = @stokes_infsup_constant;
54 
55 model.orthonormalize = @stokes_orthonormalize;
56 model.plot_sim_data = @stokes_plot_sim_data;
57 model.RB_basisgen = @stokes_rb_basisgen;
58 model.RB_extension_algorithm = @stokes_rb_supremizer_extension;
59 model.rb_reconstruction = @stokes_rb_reconstruction;
60 model.reduced_data_subset = @stokes_reduced_data_subset;
61 model.save_detailed_simulations = @save_detailed_simulations;
62 model.set_rb_in_detailed_data = @stokes_set_rb_in_detailed_data;
63 
64 model.SCMparams.NumberOfParameters = 500;
65 model.SCMparams.type='infsup';
66 
67 end
68 
69 
70 %% auxiliary functions
71 function detailed_data = stokes_gen_detailed_data(model, detailed_data)
72 
73 detailed_data.RB_info.start_time = datestr(now, '-mm_dd_yy-HH_MM_SS');
74 
75 % assemble components once
76 if ~detailed_data.operators.componentsAssembled
77  detailed_data.operators.assemble_components(model, detailed_data);
78 end
79 
80 if strncmp(model.RB_error_indicator, 'estimator', 9)
81 
82  % compute stability constants and forward to rb_basis_generation
83  if ~isfield(detailed_data, 'infsup_constant')
84  if strcmp(model.infsup_method, 'scm')
85  model.SCMparams.dirichlet_gids = detailed_data.bc_info.dirichlet_gids;
86  detailed_data.infsup_constant = SCM(model, model.SCMparams, detailed_data);
87  elseif strcmp(model.infsup_method, 'rbf_interpolant')
88  detailed_data.infsup_constant = stokes_infsup_rbf_interpolant(model, detailed_data);
89  end
90  end
91 
92  if model.has_nonlinearity && ~isfield(detailed_data, 'rho_sqr')
93  detailed_data.rho_sqr = stokes_sobolev_embedding_constant(model, detailed_data);
94  end
95 end
96 
97 % use rbmatlab functionality
98 detailed_data.SRDW = StokesReducedDataWrapper([]);
99 detailed_data = rb_basis_generation(model,detailed_data);
100 
101 if ~exist([rbmatlabtemp, '/data/RB'], 'dir')
102  mkdir([rbmatlabtemp, '/data/RB']);
103 end
104 save([rbmatlabtemp, '/data/RB/stokes', detailed_data.RB_info.start_time, ...
105  '_data.mat'], 'model', 'detailed_data');
106 end
107 
108 function W = stokes_get_inner_product_matrix(detailed_data)
109 
110 W = VecMat.CompositeMatrixDefault(2);
111 
112 W = set(W, detailed_data.df_info.values{1}.l2_inner_product_matrix ...
113  + detailed_data.df_info.values{1}.h10_inner_product_matrix, 1);
114 W = set(W, detailed_data.df_info.values{2}.l2_inner_product_matrix, 2);
115 end
116 
117 
118 function RB = rb_empty_init_data_basis(model, detailed_data)
119 RB = VecMat.CompositeMatrixDefault(2);
120 for i = 1:2
121  RB = set(RB, zeros(detailed_data.df_info.values{i}.ndofs, 0), i);
122 end
123 end
124 
125 function detailed_data = stokes_set_rb_in_detailed_data(detailed_data, ...
126  RB)
127 if isempty(RB)
128  detailed_data.RB = rb_empty_init_data_basis([], detailed_data);
129 else
130  detailed_data.RB = copy(RB);
131 end
132 end
133 
134 function res = stokes_relative_error_algorithm(uh, uN, W, model)
135 
136 if strcmp(model.RB_error_indicator, 'error_relative')
137 
138  n = size(get(W, 1), 1);
139  vh = uh(1:n);
140  ph = uh((n+1):end);
141  eN = uh - uN;
142  evN = eN(1:n);
143  epN = eN((n+1):end);
144 
145  res = (evN' * get(W, 1) * evN) / (vh' * get(W, 1) * vh);
146  res = res + (epN' * get(W, 2) * epN) / (ph' * get(W, 2) * ph);
147 elseif strcmp(model.RB_error_indicator, 'error_velocity')
148 
149  n = size(get(W, 1), 1);
150  vh = uh(1:n);
151  vN = uN(1:n);
152  res = (vh - vN)' * get(W, 1) * (vh - vN);
153 else
154 
155  res = (uh - uN)' * W * (uh - uN);
156 end
157 
158 if res < 0
159  res = 0;
160 end
161 res = sqrt(res);
162 end
163 
164 
165 function RB = stokes_rb_basisgen(model, detailed_data)
166 %function RB = stokes_rb_basisgen(model, detailed_data)
167 %
168 % Performs generation of basis by collecting snapshots for parameter list
169 % without orthogonalization.
170 
171 RB = VecMat.CompositeMatrixDefault(2);
172 
173 % loop over RB_mu_list
174 for m = 1:length(model.RB_mu_list)
175 
176  model = set_mu(model, model.RB_mu_list{m});
177  [RBext, ~] = model.RB_extension_algorithm(model, detailed_data);
178 
179  % save basis vectors
180  RB = [RB, RBext];
181 end
182 
183 end
184 
185 
186 function Xon = stokes_orthonormalize(X, W, epsilon)
187 %function Xon = stokes_orthonormalize(X, W, epsilon)
188 %
189 % For seperate orthonormalization of velocity and pressure basis vectors.
190 
191 if nargin < 4
192  epsilon = 1e-09;
193 end
194 
196 
197 if isa(X, 'VecMat.ICompositeMatrix')
198 
199  if size(X, 2) == 0
200  return
201  end
202 
203  Xv = get(X, 1);
204  Xp = get(X, 2);
205 
206  %Xvon = orthonormalize_gram_schmidt(Xv, get(W, 1), epsilon);
207  %Xpon = orthonormalize_gram_schmidt(Xp, get(W, 2), epsilon);
208  Xvon = VecMat.gram_schmidt_reiterate(Xv, get(W, 1), epsilon);
209  Xpon = VecMat.gram_schmidt_reiterate(Xp, get(W, 2), epsilon);
210 
211  Xon = set(Xon, Xvon, 1);
212  Xon = set(Xon, Xpon, 2);
213 
214 else
215  warning('stokes_model_default:stokes_orthonormalize', 'input not supported');
216 end
217 
218 end
function p = stokes_plot_sim_data(model, model_data, sim_data, plot_params)
Visualizes solution of stokes problem.
For block-diagonal matrices.
function sim_data = stokes_detailed_simulation(model, model_data)
Performs a detailed simulation of a (Navier-)stokes problem.
function save_detailed_simulations(model, model_data, M, savepath)
perform loop over detailed simulations and save results or check consistency with existing saved resu...
function detailed_data = rb_basis_generation(model, detailed_data)
reduced basis construction with different methods