2 % Thermal Block example similar as described in the book of
3 % A.T. patera and G. Rozza (just one parameter more)
6 % - div ( a(x) grad u(x)) = q(x) on Omega
7 % u(x)) = g_D(x) on Gamma_D
8 % a(x) (grad u(x)) n) = g_N(x) on Gamma_N
9 % s = l(u) linear output functional
11 % where Omega = [0,1]^2 is divided into B1 * B2 blocks
12 % QA := B1*B2. The heat conductivities are given by mu_i:
14 % -------------------------
15 % | ... | ..... | mu_QA |
16 % -------------------------
17 % | .... | ... | ..... |
18 % -------------------------
19 % | mu_1 | ... | mu_B1 |
20 % -------------------------
22 % `Gamma_D =` upper edge
23 % `Gamma_N = boundary(Omega)\Gamma_D`
25 % `a(x) = mu_i(x) if x\in B_i`
28 % `g_N = 1` on lower edge, 0 otherwise
29 % `l(u)` = average over lower edge
31 % S. Langhof and B. Haasdonk 26.2.2011
34 model = lin_stat_model_default;
39 % MD: Mein Vorschlag für simple Basis-Generierungsalgorithmen, verwendet
41 % Funktionspointers spezifizieren kann. Genauso könnte man Simple* Klassen
42 % für die anderen Interfaces, also DetailedModel, ReducedModel und
43 % ReducedData implementieren. Das fände ich aber ehrlich gesagt unschön...
44 % Die Basisgenerierung kann so simpel sein, dass es vielleicht Sinn macht,
45 % die in eine kleine Funktion auszulagern, aber die anderen Klassen haben
46 % meistens komplexere Datenstrukturen oder mehr Methoden, so dass die
47 % "Description" Klasse wieder sehr unübersichtlich und intransparent wird,
48 % wenn man das alles hier hinein packt...
51 model.RB_basisgen = @my_RB_basisgen;
52 model.RB_train_rand_seed = 100;
53 model.RB_train_size = 1000;
54 model.RB_stop_epsilon = 1e-3;
55 model.RB_stop_Nmax = 100;
60 % model.default = '3x3 demo version';
62 % number of blocks in X and Y direction
65 model.number_of_blocks = params.B1*params.B2;
67 % each block has a different heat conductivity
68 % upper left has conductivity 1
71 if ~isfield(params,
'mu_range')
72 params.mu_range = [0.1,10];
74 mu_range = params.mu_range;
75 for p = 1:model.number_of_blocks
76 mu_names = [mu_names,{[
'mu',num2str(p)]}];
77 mu_ranges = [mu_ranges,{mu_range}];
79 model.mu_names = mu_names;
80 model.mu_ranges = mu_ranges;
82 %
default values 1 everywhere
83 model.mus = ones(model.number_of_blocks,1);
84 model.get_mu = @(model) model.mus(:);
85 model.set_mu = @thermalblock_set_mu;
87 % simple snapshots without orthogonalization
88 %model.RB_generation_mode =
'lagrangian';
89 %model.RB_mu_list = {[0.1,1,1,1,1,1,1,1,1],...
90 % [1,0.1,1,1,1,1,1,1,1],...
91 % [1,1,0.1,1,1,1,1,1,1],...
92 % [1,1,1,0.1,1,1,1,1,1],...
93 % [1,1,1,1,0.1,1,1,1,1],...
94 % [1,1,1,1,1,0.1,1,1,1],...
95 % [1,1,1,1,1,1,0.1,1,1],...
96 % [1,1,1,1,1,1,1,0.1,1],...
97 % [1,1,1,1,1,1,1,1,0.1]};
99 model.RB_generation_mode =
'model_RB_basisgen';
101 %%%%%%%%% set data functions
102 model.has_diffusivity = 1;
103 model.has_output_functional = 1;
104 model.has_dirichlet_values = 1;
105 model.has_neumann_values = 1;
107 % zero dirichlet values, i.e. 1 component, Q_dir = 1
108 dirichlet_values_coefficients = @(dummy,params) [0];
109 dirichlet_values_components = @(glob,params) {zeros(size(glob,1),1)};
110 model.dirichlet_values = @(glob,params) ...
111 eval_affine_decomp_general(dirichlet_values_components, ...
112 dirichlet_values_coefficients,glob,params);
114 % 1/0 neumann values depending on edge, non parametric, i.e Q_neu = 1;
115 neumann_values_coefficients = @(dummy,params) 1; % single component
116 model.neumann_values = @(glob,params) ...
117 eval_affine_decomp_general(@thermalblock_neumann_values_components, ...
118 neumann_values_coefficients, glob, params);
120 % diffusion tensor: each row four entries a11,a_21,a_12,a_22.
121 % a11(x)=a22(x) = mu_i
if x in block i, a12=a21 = 0.
122 diffusivity_tensor_coefficients = @(dummy,params) ...
124 model.diffusivity_tensor = @(glob,params) ...
125 eval_affine_decomp_general(...
126 @thermalblock_diffusivity_tensor_components, ...
127 diffusivity_tensor_coefficients, glob,params);
129 % only useful
for detailed simulation or nonlinear outputs:
130 %model.output_functional = @output_functional_boundary_integral;
132 model.operators_output = @fem_operators_output_boundary_integral;
133 % output weight
function: simply 1 on lower edge, 0
else => Ql = 1 component
134 output_function_components = @(glob,model) {1.0*(glob(:,2)<eps)};
135 output_function_coefficients = @(glob,model) 1;
136 model.output_function = @(glob,params) ...
137 eval_affine_decomp_general(...
138 output_function_components,...
139 output_function_coefficients, glob,params);
141 model.output_integral_qdeg = 2; %
143 %%%%%%%%%% discretization settings
145 if ~isfield(params,
'numintervals_per_block')
146 params.numintervals_per_block = 5;
150 model.xnumintervals = params.numintervals_per_block*params.B1;
151 model.ynumintervals = params.numintervals_per_block*params.B2;
152 model.xrange = [0,1];
153 model.yrange = [0,1];
155 if ~isfield(params,'pdeg')
159 if ~isfield(params,'qdeg')
164 model.pdeg = params.pdeg; % degree of polynomial functions
165 model.qdeg = params.qdeg; % quadrature degree
166 model.dimrange = 1; % scalar solution
168 model.boundary_type = @thermalblock_boundary_type;
169 %model.normals = @my_normals;
171 % new plot routine, additional block boundaries are plotted
173 model.compute_output_functional = 1;
174 model.yscale_uicontrols = 0.2;
176 % for error estimation:
177 model.coercivity_alpha = @(model) min(model.mus(:));
178 model.continuity_gamma = @(model) max(model.mus(:));
179 model.enable_error_estimator = 1;
181 % for making work with current dmodel/rmodel/descr decomposition
182 model.descr = []; % model is its own description
183 model.crb_enabled = 0;
185 % local data functions:
188 %%%%%%% auxiliary functions %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
190 function res = thermalblock_boundary_type(glob,params)
191 % returns the boundary type
192 res = zeros(size(glob,1),1);
193 i = find(glob(:,1)<=1e-10);
195 i = find(glob(:,1)>=1-1e-10);
197 i = find(glob(:,2)<=1e-10);
199 i = find(glob(:,2)>=1-1e-10);
202 function model = thermalblock_set_mu(model,mu)
203 if length(mu)~=model.number_of_blocks
204 error('length of mu does not fit to number of blocks!');
208 function comp = thermalblock_neumann_values_components(glob,params)
209 % Neumann values component function
210 res = zeros(size(glob,1),1);
211 i = find(glob(:,2)<eps);
215 function res = thermalblock_diffusivity_tensor_components(glob,params)
216 % diffusivity tensor component
function
218 %
for each point in glob find global block number
219 % xi: range 0 ... B1-1
220 xi = floor(glob(:,1)*params.B1);
221 i = find(xi>=params.B1);
225 % xi: range 0 ... B1-1
226 yi = floor(glob(:,2)*params.B2);
227 i = find(yi>=params.B2);
231 block_index = yi*params.B1+xi+1;
232 zeroblock = zeros(size(glob,1),4);
233 res = cell(1,params.number_of_blocks);
234 for q = 1:params.number_of_blocks;
236 i = find(block_index==q);
245 %
for q=1:params.number_of_blocks;
246 % block = block + res{q};
248 %disp(
'check diffusivity matrix!')
251 function alpha = thermalblock_coercivity_alpha(model)
252 alpha = min(model.get_mu(model))
254 %function [RBinit] = my_init_data_basis(model, ...
258 function res = my_RB_basisgen(model,detailed_data)
259 % simple greedy with orthonormalization
262 %disp('greedy with orthonormalization');
264 sim_data = detailed_simulation(model, detailed_data);
265 n = fem_h10_norm(sim_data.uh);
266 detailed_data.RB = sim_data.uh.dofs / n;
267 reduced_data = gen_reduced_data(model, detailed_data);
268 h10_matrix = sim_data.uh.df_info.h10_inner_product_matrix;
270 mus = rand_uniform(model.RB_train_size,model.mu_ranges);
272 max_err_est_sequence = [];
273 max_mu_sequence = {};
275 while( max_err_est> model.RB_stop_epsilon) && ...
276 (size(detailed_data.RB,2)< model.RB_stop_Nmax)
279 for i = 1:size(mus,2);
280 model = set_mu(model, mus(:,i));
281 rb_sim_data = rb_simulation(model,reduced_data);
282 % disp([
'rb sim: Delta = ',num2str(rb_sim_data.Delta)]);
283 if rb_sim_data.Delta > max_err_est
285 max_err_est = rb_sim_data.Delta;
288 disp([
'max error estimator: ',num2str(max_err_est),...
289 ' for mu = ',num2str(mu_max
')]);
290 max_err_est_sequence = [max_err_est_sequence,max_err_est];
291 max_mu_sequence = [max_mu_sequence,{mu_max}];
292 model = set_mu(model,mu_max);
293 sim_data = detailed_simulation(model,detailed_data);
294 pr = sim_data.uh.dofs - ...
295 detailed_data.RB * (detailed_data.RB' * h10_matrix * ...
297 n = sqrt(max(pr
'*h10_matrix*pr,0));
298 detailed_data.RB = [detailed_data.RB,pr/n];
299 reduced_data = gen_reduced_data(model,detailed_data);
300 disp(['new basis size:
',num2str(size(detailed_data.RB,2))]);
303 RB = detailed_data.RB;
304 RB_info.max_mu_sequence = max_mu_sequence;
305 RB_info.max_err_est_sequence = max_err_est_sequence;
308 function p = thermalblock_plot_sim_data(model,detailed_data,sim_data, ...
313 %if ~isfield(plot_params,'no_lines
')
314 plot_params.no_lines = 1;
316 p = lin_stat_plot_sim_data(model,detailed_data,sim_data,plot_params);
318 if ~isfield(plot_params,'plot_blocks
')
319 plot_params.plot_blocks = 1;
321 if plot_params.plot_blocks
322 X = [0:1/model.B1:1];
323 Y = [0:1/model.B2:1];
325 [zeros(1,model.B1+1);...
326 ones(1,model.B1+1)]);
327 set(l1,'color
',[0,0,0],'linestyle
','-.
');
329 l2 = line([zeros(1,model.B2+1);...
330 ones(1,model.B2+1)],...
332 set(l2,'color
',[0,0,0],'linestyle
','-.
');
333 p = [p(:);l1(:);l2(:)];
Reduced basis implementation for linear stationary PDEs.
function ModelDescr descr = descr_default()
This function initializes the default settings for the ARE model.
A triangular conforming grid in two dimensions.
function model = thermalblock_model_struct(params)
Thermal Block example similar as described in the book of A.T. patera and G. Rozza (just one paramete...
a very simple detailed data implementation gathering several detailed snapshots spanning the reduced ...
function local_model = elliptic_discrete_model(model)
function creating a model with local functions out of a model with global functions. See detailed description for explanation.
function p = plot_sim_data(model, model_data, sim_data, plot_params)
function performing the plot of the simulation results as specified in model.
function res = thermalblock(step, params)
thermalblock example