1 function model = thermalblock_dd_model(params)
2 %
function model = thermalblock_dd_model(params)
4 % small model
for the `\mu`-dependent PDE
6 % `- \nabla (k(x;\mu) \nabla u(x;\mu)) = h(x;\mu)`, in `\Omega`,
7 % `u(x;\mu) = 0`, on `\Gamma`.
9 % - `k(x;\mu)` is the block-wise constant heatcoefficient
10 % - `u(x;\mu)` is the solution
11 % - `h(x;\mu)` is the source
function
19 if ~isfield(params,
'numintervals')
20 params.numintervals = 16;
22 if ~isfield(params,
'pdeg')
25 if ~isfield(params,'detailed_train_savepath')
26 params.detailed_train_savepath = '/dom_dec';
29 model = lin_stat_model_default;
30 model.decomp_mode = 0; % default
33 model.mu_names = {
'mu_1',
'mu_2',
'mu_3',
'mu_f'};
34 model.mu_ranges = {[0.1,10],[0.1,10],[0.1,10],[0,1]};
36 model.RB_generation_mode =
'lagrangian';
37 model.detailed_train_savepath = params.detailed_train_savepath;
38 model.filecache_ignore_fields_in_model = {};
40 model.RB_mu_list = {[0.1,0.1,0.1,1],[0.1,0.1,10,1],[0.1,10,0.1,1], ...
41 [0.1,10,10,1],[10,0.1,10,1],[10,10,0.1,1],[0.1, ...
42 5.05,5.05,0],[10,0.1,0.1,1],[0.1,5.05,10,1], ...
43 [10,0.1,10,0],[10,10,0.1,0],[5.05,0.1,0.1,1], ...
44 [0.1,10,5.05,1],[0.1,10,10,0],[5.05,10,10,1], ...
45 [10,0.1,5.05,1],[0.1,10,5.05,0],[10,0.1,0.1,0], ...
46 [0.1,0.1,10,0],[0.1,10,0.1,0],[10,5.05,0.1,1], ...
47 [5.05,0.1,5.05,1],[5.05,5.05,0.1,1],[10,5.05, ...
48 0.1,0],[5.05,10,10,0],[0.1,5.05,0.1,1]};
50 %
default parameter values:
51 model = set_mu(model,[1,1,1,0.35]);
53 % data
function specification:
54 model.has_diffusivity = 1;
55 model.has_advection = 0;
56 model.has_reaction = 0;
58 model.has_output_functional = 0;
59 model.has_dirichlet_values = 1;
60 model.has_neumann_values = 0;
61 model.has_robin_values = 0;
66 source_coefficients = @(grid,eindices,loc,params) 2*[params.mu_f; ...
68 source_components = @my_source_components;
70 model.source = @(grid,eindices,loc,params) ...
71 eval_affine_decomp_general(source_components,source_coefficients, ...
72 grid,eindices,loc,params);
75 diffusivity_tensor_coefficients = @(grid,eincices,loc,params) ...
76 [params.mu_1; params.mu_2; params.mu_3; 1];
77 diffusivity_tensor_components = @my_diffusivity_tensor_components;
79 model.diffusivity_tensor = @(grid,eindices,loc,params) ...
80 eval_affine_decomp_general(diffusivity_tensor_components, ...
81 diffusivity_tensor_coefficients, ...
82 grid,eindices,loc,params);
85 dirichlet_coefficients = @(grid,eindices,face_index,lloc,params) 0;
86 dirichlet_components = @(grid,eindices,face_index,lloc,params) ...
87 { zeros(size(eindices,1),1) };
89 model.dirichlet_values = @(grid,eindices,face_index,lloc,params) ...
90 eval_affine_decomp_general(dirichlet_components, ...
91 dirichlet_coefficients,grid, ...
92 eindices,face_index,lloc,params);
95 model.gridtype =
'triagrid';
96 model.xnumintervals = params.numintervals;
97 model.ynumintervals = params.numintervals;
101 % discretization parameters:
102 model.pdeg = params.pdeg;
103 model.qdeg = params.pdeg * 3;
106 %
for error estimation:
107 model.coercivity_alpha = @my_coercivity_alpha;
108 model.continuity_gamma = @my_continuity_gamma;
110 %%%%%%%%%%%% auxiliary functions:
112 function res = my_source_components(grid,eindices,loc,params)
113 glob = local2global(grid,eindices,loc,params);
114 res = { exp( -params.gamma_1 * ( (glob(:,1) - 0.5).^2 + ...
115 (glob(:,2) - 0.5).^2 ) ), ...
116 exp( -params.gamma_2 * ( (glob(:,1) - 0.875).^2 + ...
117 (glob(:,2) - 0.875).^2 ) ) };
119 function res = my_diffusivity_tensor_components(grid,eindices,loc,params)
120 CX = grid.CX(eindices);
121 CY = grid.CY(eindices);
122 res = {[(CX<0.5) & (CY<0.5), zeros(length(CX),2), (CX<0.5) & (CY<0.5)], ...
123 [(CX>0.5) & (CY<0.5), zeros(length(CX),2), (CX>0.5) & (CY<0.5)], ...
124 [(CX<0.5) & (CY>0.5), zeros(length(CX),2), (CX<0.5) & (CY>0.5)], ...
125 [(CX>0.5) & (CY>0.5), zeros(length(CX),2), (CX>0.5) & (CY>0.5)] ...
128 function res = my_coercivity_alpha(model)
129 mus = model.get_mu(model);
130 res = min([1;mus(1:end-1)]);
132 function res = my_continuity_gamma(model)
133 mus = model.get_mu(model);
134 res = max([1;mus(1:end-1)]);