rbmatlab  1.16.09
 All Classes Namespaces Files Functions Variables Modules Pages
advection_fv_output_opt_model.m
1 function model = advection_fv_output_opt_model(params)
2 %function model = advection_fv_output_model(params);
3 %
4 % model for the new advection model with functional output
5 % two overlapping velocity fields and boundary value modification
6 % are the parameters. Output functional is average concentration
7 % in subdomain. Pure explicit discretization of convection with ldg
8 % discretization
9 %
10 % example of a model not as class but as structure plus function pointers
11 % goal: library works with structure- or class-models
12 %
13 % possible fields of params:
14 % coarse_factor: coarsening factor between 1 and 8
15 % time-steps, gridsize are scaled down with this.
16 %
17 % The model is adapted to work with optimization.
18 
19 % M. Dihlmann 03.02.2010
20 
21 model = [];
22 if (nargin ==1) & (isfield(params,'coarse_factor'))
23  coarse_factor = params.coarse_factor;
24 else
25  coarse_factor = 1;
26 end;
27 
28 % specification of the time information
29 model.T = 1;
30 model.nt = 2048/coarse_factor;
31 %model.save_time_indices = 0:16/coarse_factor:model.nt;
32 model.save_time_indices = 0:8/coarse_factor:model.nt;
33 disp('nt to be adjusted later!');
34 model.dt = model.T/model.nt;
35 model.verbose = 9;
36 model.debug = 0;
37 %model.debug = 1; %
38 %if model.debug
39 % disp('model with debugging turned on.');
40 %end;
41 model.axis_equal = 1;
42 model.error_estimation = 1;
43 
44 % grid information
45 %model.grid_initfile = 'rectangle_triagrid.mat';
46 model.gridtype = 'triagrid';
47 model.xnumintervals = 256/coarse_factor;
48 model.ynumintervals = 128/coarse_factor;
49 model.xrange = [0,2];
50 model.yrange = [0,1];
51 % set completely dirichlet boundary
52 model.bnd_rect_corner1 = [0,0]-eps;
53 model.bnd_rect_corner2 = [2,1]+eps;
54 %model.opts.bnd_rect_index = [-1];
55 model.element_quadrature = @triaquadrature;
56 model.intersection_quadrature = @intervalquadrature;
57 model.dim_U = model.xnumintervals * model.ynumintervals * 2;
58 
59 % Main global function pointers expected in every model
60 model.gen_model_data = @lin_evol_gen_model_data;
61 model.gen_detailed_data = @lin_evol_gen_detailed_data;
62 model.gen_reduced_data = @lin_evol_gen_reduced_data;
63 model.detailed_simulation = @lin_evol_detailed_simulation;
64 model.rb_simulation = @lin_evol_rb_simulation;
65 model.rb_reconstruction = @lin_evol_rb_reconstruction;
66 model.plot_sim_data = @lin_evol_plot_sim_data;
67 
68 
69 % Dirichlet and Initial data
70 model.dirichlet_values_ptr = @dirichlet_values_affine_decomposed;
71 model.dirichlet_values_coefficients_ptr = @my_dirichlet_values_coefficients;
72 model.dirichlet_values_components_ptr = @my_dirichlet_values_components;
73 % dummy one dirichlet_values:
74 %model.dirichlet_values_coefficients_ptr = @(params) 1;
75 %model.dirichlet_values_components_ptr = @(glob,params) {ones(1, ...
76 % size(glob,2))};
77 model.rb_init_data_basis = @RB_init_data_basis;
78 model.set_rb_in_detailed_data = @(detailed_data,RB) ...
79  setfield(detailed_data,'RB',RB);%new
80  model.get_rb_size = @(model,detailed_data)size(detailed_data.RB,2);%new
81 model.init_values_ptr = @init_values_affine_decomposed;
82 model.init_values_coefficients_ptr = @my_dirichlet_values_coefficients_t0;
83 model.init_values_components_ptr = @my_dirichlet_values_components;
84 model.cone_number = 3; % number of components = cones
85 model.cone_range = [0,1]; % x-range of cone-support
86 model.cone_weight = 1; % leftmost cone with full weight
87 model.velocity_ptr = @velocity_affine_decomposed;
88 model.velocity_coefficients_ptr = @my_velocity_coefficients;
89 model.velocity_components_ptr = @my_velocity_components;
90 %model.vx_weight = 1.0;
91 model.vx_weight = 0.75;
92 %model.vy_weight = 1.0;
93 model.vy_weight = 0;
94 model.conv_flux_ptr = @conv_flux_linear;
95 
96 
97 
98 % perhaps these are redundant later...
99 model.divclean_mode = 0;
100 model.flux_quad_degree = 1;
101 model.flux_linear = 1;
102 
103 model.init_values_algorithm = @disc_init_values;
104 model.init_values_qdeg = 0;
105 model.pdeg = 0; % FV schemes with piecewise constant ansatz functions
106 model.evaluate_basis = @fv_evaluate_basis;
107 model.l2project = @l2project;
108 model.local_mass_matrix = @fv_local_mass_matrix_tria; % triangular grid
109 model.mass_matrix = @fv_mass_matrix;
110 model.ndofs_per_element = 1;
111 model.plot = @fv_plot; % then also plot_sequence will work
112 
113 % for use of old datafunctions, use the following:
114 %model.operators_ptr = @fv_operators_implicit_explicit;
115 %model.operators_conv_explicit = @fv_operators_conv_explicit;
116 %model.operators_conv_implicit = @fv_operators_conv_implicit;
117 %model.operators_diff_explicit = @fv_operators_diff_explicit;
118 %model.operators_diff_implicit = @fv_operators_diff_implicit;
119 %model.operators_neumann_explicit = @fv_operators_neumann_explicit_old;
120 %model.operators_neumann_implicit = @fv_operators_neumann_implicit;
121 
122 model.operators_ptr = @fv_operators_implicit_explicit;
123 % for use of new data function access, use now
124 model.operators_conv_explicit = @fv_operators_conv_explicit_engquist_osher;
125 model.operators_conv_implicit = @fv_operators_zero;
126 model.operators_diff_explicit = @fv_operators_zero;
127 model.operators_diff_implicit = @fv_operators_zero;
128 model.operators_neumann_implicit = @fv_operators_zero;
129 model.operators_neumann_explicit = @fv_operators_zero;
130 model.operators_output = @fv_operators_output;
131 %model.operators_diff_implicit = @fv_operators_diff_implicit_gradient;
132 %model.operators_neumann_explicit = @fv_operators_neumann_explicit;
133 
134 %model.flux_linear = 1;
135 model.data_const_in_time = 0; % time varying data...
136 %model.diffusivity_ptr = ...;
137 %model.neumann_value_ptr = @...
138 model.affinely_decomposed = 1; % data functions allow affine
139  % parameter decomposition
140 
141 model.compute_output_functional = 1; % turn on computation of output
142 model.output_function_ptr = @output_function_box_mean;
143 %model.sbox_xmin = 1;
144 model.sbox_xmin = 1;
145 model.sbox_xmax = 2;
146 model.sbox_ymin = 0;
147 %model.sbox_ymax = 0.5;
148 model.sbox_ymax = 0.5;
149 
150 %model.disc_element_mean_ptr = @fv_element_mean;
151 
152 % RB information:
153 model.mu_names = {'cone_weight','vx_weight','vy_weight'};
154 model.mu_ranges = {[0 1],[0 1],[0,1]};
155 %model.mu_initial_values = [0.5, 0.5, 0.5]; %Initial values for parameter optimization
156 %model.mu_set_for_opt = [0, 1, 1] %1: Parameter should be optmized, 2: no optimization for this parameter
157 model.RB_numintervals = [1,1,1];
158 model.RB_generation_mode = 'greedy_uniform_fixed';
159 model.RB_error_indicator = 'estimator';
160 model.RB_extension_algorithm = @RB_extension_PCA_fixspace;
161 model.RB_stop_epsilon = 1e-2;
162 model.RB_stop_timeout = 24*60*60; % 1 minute
163 model.RB_stop_Nmax = 4;
164 
165 model.get_inner_product_matrix = @(model_data) model_data.W;
166 
167 model.rb_init_values = @rb_init_values_default;
168 model.init_values_algorithm = @fv_init_values;
169 model.rb_operators = @lin_evol_rb_operators;
170 model.rb_simulation = @lin_evol_rb_simulation;
171 model.rb_reconstruction = @rb_reconstruction_default;
172 model.reduced_data_subset = @lin_evol_reduced_data_subset;
173 model.error_algorithm = @fv_error;
174 model.error_norm = 'l2';
175 model.L_I_inv_norm_bound = 1;
176 model.L_E_norm_bound = 1;
177 model.get_estimator_from_sim_data = @(sim_data) sim_data.Delta(end);
178 model.get_dofs_from_sim_data = @(sim_data) sim_data.U;
179 model.filecache_ignore_fields_in_model = {'N','Nmax','mu_ranges'};
180 model.filecache_ignore_fields_in_detailed_data = {'RB_info'};
181 model.get_rb_from_detailed_data = @(detailed_data) detailed_data.RB;
182 model.set_rb_in_detailed_data = @(detailed_data,RB) ...
183  setfield(detailed_data,'RB',RB);
184 model.PCA_fixspace = @PCA_fixspace;
185 
186 
187 
188 
189 
190 %Optimization
191 
192 model.optimization.init_params = [0.5,0.5,0.5];
193 model.optimization.params_to_optimize = [0,1,1];
194 model.optimization.opt_mode = 'detailed'; %'reduced'
195 model.optimization.optimizer = @detailed_grid_search;
196 %model.opt_method='grid-search';
197 %model.opt_param=[] %Parameters for optimization method
198 model.optimization.opt_params.grid_density = [6,6,6]; %Density of the search grid. F.Ex. "3" means using 3 parameters per dimension
199 model.optimization.objective_function = @lin_evol_get_output_detailed;
200 
201 % set parameter setting function
202 model.set_time = @set_time;
203 model.set_mu = @set_mu_default;%@set_mu_lin_evol_opt;
204 model.get_mu = @get_mu_default;
205 
206 
207 %return; % temporary end
208 
209 %model.L_I_inv_norm_bound = 1; % bounds for implicit/explicit operator
210 %model.L_E_norm_bound = 1;
211 
212 %model.neumann_values_ptr = 0;
213 
214 %model.name_init_values = 'decomp_function_ptr';
215 % implement below!!
216 %model.init_values_coefficients_ptr = @my_init_values_coefficients;
217 %model.init_values_components_ptr = @my_init_values_components;
218 
219 %model.name_diffusive_num_flux = 'none';
220 
221 
222 %name_diffusive_num_flux = 'gradient';
223  % precomputed, divcleaned velocity field
224 % name_flux = 'gdl2';
225 %
226 % lambda = 2.0e-7; % v = - lambda * grad p
227 % name_convective_num_flux = 'lax-friedrichs';
228 % inner_product_matrix_algorithm = @fv_inner_product_matrix;
229 % verbose = 5;
230 %
231 %model.rb_init_values = @rb_init_values;
232 %
233 %% further method pointers, which are specific to model:
234 model.orthonormalize = @model_orthonormalize_gram_schmidt;
235 model.inner_product = @fv_inner_product;
236 %model.PCA = @model_PCA_fixspace;
237 %model.cached_detailed_simulation = @cached_detailed_simulation;
238 %model.save_detailed_simulations = @save_detailed_simulations;
239 %model.load_detailed_simulation = @load_detailed_simulation;
240 %model.use_velocity_matrixfile = 1;
241 %model.divclean_mode = 'none'; % file is already cleaned by optimization
242 
243 % set matrix-file name for optional generation or reading of file
244 %model.velocity_matrixfile = ['vel_', model.name_flux,'_',...
245 % num2str( model.xnumintervals),'x',...
246 % num2str( model.ynumintervals),...
247 % '_l',num2str( model.lambda),'.mat'];%
248 %
249 %model.lxf_lambda = 1.0194e+003;
250 
251 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
252 %%%%% auxiliary functions used as pointers above:
253 %%%%% check later, if they are of general interest to be exported
254 %%%%% as standalone functions
255 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
256 
257 % new function syntax: global coordinates as rows in glob_coord,
258 % time and parameter variables in params
259 
260 % dirichlet-data: convex combination of equidistant cones, decaying
261 % in time
262 function res = my_dirichlet_values_coefficients_t0(params);
263 params.t = 0;
264 res = my_dirichlet_values_coefficients(params);
265 
266 function res = my_dirichlet_values_coefficients(params);
267 % res is a vector of coefficients
268 Q_0 = params.cone_number;
269 res = zeros(1,Q_0);
270 max_pos = params.cone_weight * (Q_0-1)+1;
271 t = params.t;
272 for q = 1:Q_0
273  res(q) = (1-(max_pos-q))*(1-t) * ((max_pos>=q) && (max_pos < q+1)) ...
274  + (1+(max_pos-q))*(1-t) * ((max_pos>=q-1) && (max_pos < q));
275 end;
276 
277 function res = my_dirichlet_values_components(glob,params);
278 % res is a cell-array of row-vectors of global evaluations
279 Q_0 = params.cone_number;
280 res = cell(1,Q_0);
281 %delta_cone = 1/(Q_0+1);
282 delta_cone = (params.cone_range(2)-params.cone_range(1))/(Q_0+1);
283 cone_pos_x = delta_cone * (1:Q_0)+ params.cone_range(1);
284 for q = 1:Q_0
285  res{q} = 1-min(sqrt((glob(:,1)-cone_pos_x(q)).^2+...
286  (glob(:,2)-1).^2)*(Q_0+1),1);
287 end;
288 
289 % velocity field: overlap of y- and x parabolic profiles
290 function res = my_velocity_coefficients(params);
291 res = [params.vx_weight, params.vy_weight]*(1-params.t);
292 
293 function res = my_velocity_components(glob,params);
294 resx = [5*(1-glob(:,2).^2),...
295  zeros(size(glob,1),1)];
296 %resy = [zeros(1,size(glob,2));...
297 % -2*((1-glob(1,:).^2 .*(glob(1,:)>=0) & (glob(1,:)<=1)...
298 % ))];
299 resy = [zeros(size(glob,1),1),...
300  -1*((4-glob(:,1).^2))];
301 res = {resx, resy};
302 
303 %function res = my_velocity_coefficients2(params);
304 %res = 1;
305 
306 %function res = my_velocity_components2(glob,params);
307 %res = {repmat([0;-1],1,size(glob,2))};
308 
function [ sim_data , tictoc ] = load_detailed_simulation(m, savepath, params)
load single trajectory of previously saved results.
function res = intervalquadrature(poldeg, func, varargin)
integration of function func over reference interval == unit interval. by Gaussian quadrature exactly...
function r = verbose(level, message, messageId)
This function displays messages depending on a message-id and/or a level. Aditionally you can set/res...
Definition: verbose.m:17
function [ RBext , dummy ] = RB_extension_PCA_fixspace(model, detailed_data)
function computing a RB basis extension for given parameters by the POD-Greedy algorithm.
function [ L_I_diff , bdir_I_diff ] = fv_operators_diff_implicit_gradient(model, model_data, U, NU_ind)
computes diffusion contributions to finite volume time evolution matrices, or their Frechet derivati...
function [ L_E_conv , bdir_E_conv ] = fv_operators_conv_explicit_engquist_osher(model, model_data, U, NU_ind)
computes convection contribution to finite volume time evolution matrices, or their Frechet derivati...
function Udirichlet = dirichlet_values(model, X, Y)
UDIRICHLET = DIRICHLET_VALUES([X],[Y], MODEL) Examples dirichlet_values([0,1,2],[1,1,1],struct(name_dirichlet_values, homogeneous, ... c_dir, 1)) dirichlet_values([0:0.1:1],[0],struct(name_dirichlet_values, xstripes, ... c_dir, [0 1 2], ... dir_borders, [0.3 0.6])) dirichlet_values([0:0.1:1],[0],struct(name_dirichlet_values, box, ... c_dir, 1, ... dir_box_xrange, [0.3 0.6], ... dir_box_yrange, [-0.1 0.1]))
function [ L_E_neu , b_E_neu ] = fv_operators_neumann_explicit(model, model_data, U, NU_ind)
computes a neumann contribution matrix for finite volume time evolution operators, or their Frechet derivative
function rb_sim_data = rb_reconstruction_default(model, detailed_data, rb_sim_data)
(trivial) function computing a detailed reconstruction by linear combination of the coefficients in t...
function p = fv_plot(gridbase grid, dofs, params)
routine plotting a single fv function of fv_functions.
Definition: fv_plot.m:17
function a0 = rb_init_values_default(model, detailed_data)
function computing the reduced basis initial values. If the decomposition mode is coefficients...
function [ v , l2norm ] = fv_operators_output(model, model_data)
function returning components, coefficients, and complete operator for a linear output functional on ...
function p = plot_sequence(varargin)
plotting a sequence of data slices on polygonal 2d grid (constructed from params if empty) and provid...
Definition: plot_sequence.m:17
function y = lin_evol_get_output_detailed(model, varargin)
lin_evol_get_output_detailed(model, varargin)
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 Umean = fv_element_mean(model, model_data, U, I)
function computing the element averages of a discrete function U in the grid elements with indices I...
function reduced_data_subset = lin_evol_reduced_data_subset(model, reduced_data)
method which modifies reduced_data, which is the data, that will be passed to the online-simulation a...
function [ flux , lambda ] = conv_flux_linear(glob, U, params)
function computing the convective flux of a convection problem.