1 function model = advection_fv_output_opt_model(params)
2 %
function model = advection_fv_output_model(params);
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
10 % example of a model not as
class but as structure plus function pointers
11 % goal: library works with structure- or
class-models
13 % possible fields of params:
14 % coarse_factor: coarsening factor between 1 and 8
15 % time-steps, gridsize are scaled down with this.
17 % The model is adapted to work with optimization.
19 % M. Dihlmann 03.02.2010
22 if (nargin ==1) & (isfield(params,
'coarse_factor'))
23 coarse_factor = params.coarse_factor;
28 % specification of the time information
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;
39 % disp(
'model with debugging turned on.');
42 model.error_estimation = 1;
45 %model.grid_initfile =
'rectangle_triagrid.mat';
46 model.gridtype =
'triagrid';
47 model.xnumintervals = 256/coarse_factor;
48 model.ynumintervals = 128/coarse_factor;
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;
57 model.dim_U = model.xnumintervals * model.ynumintervals * 2;
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;
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;
74 %model.dirichlet_values_coefficients_ptr = @(params) 1;
75 %model.dirichlet_values_components_ptr = @(glob,params) {ones(1, ...
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;
98 % perhaps these are redundant later...
99 model.divclean_mode = 0;
100 model.flux_quad_degree = 1;
101 model.flux_linear = 1;
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;
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;
122 model.operators_ptr = @fv_operators_implicit_explicit;
123 %
for use of
new data
function access, use now
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;
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
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;
147 %model.sbox_ymax = 0.5;
148 model.sbox_ymax = 0.5;
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';
161 model.RB_stop_epsilon = 1e-2;
162 model.RB_stop_timeout = 24*60*60; % 1 minute
163 model.RB_stop_Nmax = 4;
165 model.get_inner_product_matrix = @(model_data) model_data.W;
168 model.init_values_algorithm = @fv_init_values;
169 model.rb_operators = @lin_evol_rb_operators;
170 model.rb_simulation = @lin_evol_rb_simulation;
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;
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
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;
207 %
return; % temporary end
209 %model.L_I_inv_norm_bound = 1; % bounds
for implicit/
explicit operator
210 %model.L_E_norm_bound = 1;
212 %model.neumann_values_ptr = 0;
214 %model.name_init_values =
'decomp_function_ptr';
216 %model.init_values_coefficients_ptr = @my_init_values_coefficients;
217 %model.init_values_components_ptr = @my_init_values_components;
219 %model.name_diffusive_num_flux =
'none';
222 %name_diffusive_num_flux =
'gradient';
223 % precomputed, divcleaned velocity field
224 % name_flux =
'gdl2';
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;
231 %model.rb_init_values = @rb_init_values;
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;
240 %model.use_velocity_matrixfile = 1;
241 %model.divclean_mode =
'none'; % file is already cleaned by optimization
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'];%
249 %model.lxf_lambda = 1.0194e+003;
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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
257 %
new function syntax: global coordinates as rows in glob_coord,
258 % time and parameter variables in params
260 % dirichlet-data: convex combination of equidistant cones, decaying
262 function res = my_dirichlet_values_coefficients_t0(params);
264 res = my_dirichlet_values_coefficients(params);
266 function res = my_dirichlet_values_coefficients(params);
267 % res is a vector of coefficients
268 Q_0 = params.cone_number;
270 max_pos = params.cone_weight * (Q_0-1)+1;
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));
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;
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);
285 res{q} = 1-min(sqrt((glob(:,1)-cone_pos_x(q)).^2+...
286 (glob(:,2)-1).^2)*(Q_0+1),1);
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);
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)...
299 resy = [zeros(size(glob,1),1),...
300 -1*((4-glob(:,1).^2))];
303 %
function res = my_velocity_coefficients2(params);
306 %
function res = my_velocity_components2(glob,params);
307 %res = {repmat([0;-1],1,size(glob,2))};
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...
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.
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...
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.