rbmatlab  1.16.09
 All Classes Namespaces Files Functions Variables Modules Pages
elliptic_debug_model.m
1 function model = elliptic_debug_model(params);
2 %function model = elliptic_debug_model(params);
3 %
4 % small example of a model, i.e. a structure describing the data
5 % functions and geometry information of a general elliptic equation consisting
6 % of diffusion, convection, reaction equation:
7 %
8 % - div ( a(x) grad u(x)) + div (b(x) u(x)) + c(x) u(x) = f(x) on Omega
9 % u(x)) = g_D(x) on Gamma_D
10 % a(x) (grad u(x)) n) = g_N(x) on Gamma_N
11 % alpha(x) u(x) + beta(x) a(x) (grad u(x)) n) = g_R(x) on Gamma_R
12 %
13 % s = l(u) linear output functional
14 %
15 % Here, we denote the functions as
16 % u: scalar 'solution' (if known, for validation purpose)
17 % f: scalar 'source'
18 % a: tensor valued 'diffusivity_tensor'
19 % b: vector valued 'velocity'
20 % c: scalar 'reaction'
21 % g_D: scalar 'dirichlet_values'
22 % g_N: scalar 'neumann_values'
23 % g_R: scalar 'robin_values'
24 % alpha: scalar 'robin_alpha'
25 % beta: scalar 'robin_beta'
26 %
27 % Each function allows the evaluation in many points
28 % simultaneuously by
29 %
30 % model.source(glob)
31 % or model.source(glob,params)
32 %
33 % where glob is a n times 2 matrix of row-wise points. The result
34 % is a n times 1 vector of resulting values of f.
35 %
36 % model.diffusivity_tensor(glob)
37 % or model.diffusivity_tensor(glob,params)
38 %
39 % where glob is a n times 2 matrix of row-wise points. The result
40 % is a n times 4 matrix of resulting values of diffusivity, where
41 % the values of a are sorted in matlab-order as [a_11, a_21, a_12, a_22];
42 %
43 % additionally, the model has a function, which determines, whether
44 % a point lies on a Dirichlet, Neumann or Robin boundary:
45 %
46 % model.boundary_type(glob)
47 % 0 no boundary (inner edge or point)
48 % -1 indicates Dirichlet-boundary
49 % -2 indicates Neumann-boundary
50 % -3 indicates Robin-boundary
51 %
52 % Additionally, the normals in a boundary point can be requested by
53 %
54 % model.normals(glob)
55 % Here, glob are assumed to be boundary points lying ON THE
56 % INTERIOR of an edge, such that the outer unit normal is well-defined.
57 %
58 % The latter 2 methods boundary_type() and normals() need not be
59 % implemented for grid-based methods, as the normals simply can be
60 % obtained by the grid. The methods are only required, if using
61 % meshless methods with data functions that use normals.
62 %
63 % The output functional must be specified
64 % s = model.output_functional(model, model_data, u)
65 % where u is a dof vector of the discrete solution. model_data is
66 % assumed to contain the grid
67 %
68 % possible fields of params:
69 % numintervals: the unit square is divided into numintervals
70 % intervals per side. default is 10;
71 %
72 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
73 %
74 % The data functions given in this model are a parametrized model
75 % to be used for debugging purpose. The solution for all parameters
76 % is identical, but the data functions change.
77 % The model parameters are mu_a in [-1,1], mu_b in [0,1] and mu_c
78 % in [0,1]. This model should be used for validating any new
79 % elliptic solver BY PARAMETER VARIATION.
80 % by mu_b=mu_c = 0 the advection and reaction terms can be turned
81 % off. by params.all_dirichlet_boundary = 1, all boundary types can
82 % be set to dirichlet to turn off Neumann and Robin boundaries.
83 %
84 % The data functions
85 % are
86 %
87 % a = [1,0;0,1] + mu_a [0, 1/2; 1/2, 0];
88 % b = [1/2;1/2] (x-y) mu_b;
89 % c = (1-x)mu_c
90 %
91 % the exact solution is assumed to be u(x) = sin(k_x x) sin(k_y y)
92 %
93 % The computational domain is the unit interval. Dirichlet
94 % boundaries assigned right and upper. Neumann assigned left and
95 % Robin assigned at bottom.
96 %
97 % The source term, the boundary values are then specified by the PDE.
98 %
99 % An output functional is simply the average over the complete domain
100 
101 % B. Haasdonk 24.1.2011
102 
103 if nargin == 0
104  params = [];
105 end;
106 
107 if ~isfield(params,'numintervals')
108  params.numintervals = 10; % 2 x 2 grid! with 8 triangles
109 end;
110 
111 if ~isfield(params,'pdeg')
112  params.pdeg = 1;
113 end;
114 
115 if ~isfield(params,'qdeg')
116  params.qdeg = 2;
117 end;
118 
119 model = [];
120 model = lin_stat_model_default;
121 
122 model.kx = 4;
123 model.ky = 4;
124 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
125 % Play with the following parameters for debugging a discretization
126 %%% initial setting: simple poisson problem:
127 %model.mu_a = 0;
128 %model.mu_b = 0;
129 %model.mu_c = 0;
130 %model.all_dirichlet_boundary = 1;
131 %%% more complex volume integrals
132 %model.mu_a = 1;
133 %model.mu_b = 1;
134 %model.mu_c = 1;
135 %model.all_dirichlet_boundary = 1;
136 %%% robin and neumann boundary conditions and complex volume integrals
137 model.mu_a = 1;
138 model.mu_b = 1;
139 model.mu_c = 1;
140 model.all_dirichlet_boundary = 0;
141 
142 % parameter settings for rb-methods:
143 model.mu_names ={'mu_a','mu_b','mu_c'};
144 model.mu_ranges ={[0,1],[0,1],[0,1]};
145 model.RB_mu_list = {[0,0,0],[0,0,1],[0,1,0],[0,1,1],[1,0,0],...
146  [1,0,1],[1,1,0],[1,1,1]};
147 model.RB_generation_mode = 'lagrangian';
148 
149 %%%%%%% data function specification:
150 
151 % if any of the following flags is set, the corresponding function
152 % pointers must be provided below.
153 
154 model.has_reaction = 1;
155 %model.has_reaction = 0;
156 model.has_source = 1;
157 model.has_reaction = 1;
158 %model.has_reaction = 0;
159 model.has_advection = 1;
160 %model.has_advection = 0;
161 model.has_diffusivity = 1;
162 model.has_output_functional = 1;
163 model.has_dirichlet_values = 1;
164 model.has_neumann_values = 1;
165 %model.has_neumann_values = 0;
166 %model.has_robin_values = 1;
167 model.has_robin_values = 1;
168 
169 % solution is known: for validation purpose
170 model.solution = @(glob,params) ...
171  sin(glob(:,1)*params.kx).* ...
172  sin(glob(:,2)*params.ky);
173 % gradient: each row one gradient
174 model.solution_gradient = @(glob,params) ...
175  [params.kx*cos(glob(:,1)*params.kx).* sin(glob(:,2)*params.ky), ...
176  params.ky*sin(glob(:,1)*params.kx).* cos(glob(:,2)*params.ky) ...
177  ];
178 % hessian: each row: u_xx u_yx, u_xy, u_yy
179 model.solution_hessian = @(glob,params) ...
180  [-params.kx*params.kx*sin(glob(:,1)*params.kx).* sin(glob(:,2)*params.ky), ...
181  params.kx * params.ky * cos(glob(:,1)*params.kx).* cos(glob(:,2)*params.ky) ...
182  params.kx* params.ky * cos(glob(:,1)*params.kx).* cos(glob(:,2)*params.ky), ...
183  -params.ky*params.ky * sin(glob(:,1)*params.kx).* sin(glob(:,2)*params.ky) ...
184  ];
185 
186 % params is assumed to be the model
187 % c = (1-x) mu_c + x (1-mu_c)
188 reaction_coefficients = @(dummy,params) ...
189  [params.mu_c; 1-params.mu_c];
190 reaction_components = @(glob,params) ...
191  {1-glob(:,1), glob(:,1)};
192 model.reaction = @(glob,params) ...
193  eval_affine_decomp_general(reaction_components, ...
194  reaction_coefficients,glob,params);
195 
196 % b = [1/2; 1/2] *x * mu_b, single component
197 velocity_coefficients = @(dummy,params) ...
198  [params.mu_b];
199 velocity_components = @(glob,params) ...
200  {ones(size(glob))*0.5.*repmat([glob(:,1)-glob(:,2)],1,2)};
201 model.velocity = @(glob,params) ...
202  eval_affine_decomp_general(velocity_components, ...
203  velocity_coefficients, glob, params);
204 
205 % a = eye(2) + mu_a [0,1/2; 1/2, 0], two components
206 diffusivity_tensor_coefficients = @(dummy,params) ...
207  [1; params.mu_a];
208 diffusivity_tensor_components = @(glob,params) ...
209  {...
210  [ones(size(glob,1),1),...
211  zeros(size(glob,1),1),...
212  zeros(size(glob,1),1),...
213  ones(size(glob,1),1)],...
214  0.5 * [zeros(size(glob,1),1),...
215  ones(size(glob,1),1),...
216  ones(size(glob,1),1),...
217  zeros(size(glob,1),1)] ...
218  };
219 model.diffusivity_tensor = @(glob,params) ...
220  eval_affine_decomp_general(diffusivity_tensor_components, ...
221  diffusivity_tensor_coefficients, glob,params);
222 
223 % dirichlet values:
224 % simply exact solution, not parameter dependent, 1 component!
225 dirichlet_values_coefficients = @(dummy,params) [1];
226 dirichlet_values_components = @(glob,params) ...
227  {model.solution(glob,params)};
228 model.dirichlet_values = @(glob,params) ...
229  eval_affine_decomp_general(dirichlet_values_components, ...
230  dirichlet_values_coefficients,glob,params);
231 
232 % neumann_values:
233 % simply exact solution, multiplied with diffusion-tensor and
234 % normal, hence coefficients and component number coincide
235 neumann_values_coefficients = @(dummy,params) ...
236  diffusivity_tensor_coefficients(dummy,params);
237 model.neumann_values = @(glob,params) ...
238  eval_affine_decomp_general(@my_neumann_values_components, ...
239  neumann_values_coefficients, glob, params);
240 
241 % robin_values:
242 % simply weighted sum of dirichlet and neumann values
243 model.robin_alpha = @(glob,params) ones(size(glob,1),1);
244 model.robin_beta = @(glob,params) ones(size(glob,1),1);
245 % first dirichlet coefficients, then neumann coefficients
246 model.robin_values = @(glob,params) ...
247  eval_affine_decomp_general(@my_robin_values_components, ...
248  @my_robin_values_coefficients,glob,params);
249 
250 % source terms: obtain from data functions
251 model.source = @(glob,params) ...
252  eval_affine_decomp_general(@my_source_components, ...
253  @my_source_coefficients,glob,params);
254 
255 % output functional, e.g. average over unit-square:
256 model.output_functional = @output_functional_volume_integral;
257 model.output_function = @output_function_box_mean;
258 model.sbox_xmin = 0;
259 model.sbox_ymin = 0;
260 model.sbox_xmax = 1;
261 model.sbox_ymax = 1;
262 model.output_integral_qdeg = 2; %
263 
264 %%%%%%% geometry specification and discretization:
265 
266 model.gridtype = 'triagrid';
267 model.xnumintervals = params.numintervals;
268 model.ynumintervals = params.numintervals;
269 model.xrange = [0,1];
270 model.yrange = [0,1];
271 
272 % numerics settings:
273 model.pdeg = params.pdeg; % degree of polynomial functions
274 model.qdeg = params.qdeg; % quadrature degree
275 model.dimrange = 1; % scalar solution
276 
277 % The following 2 methods boundary_type() and normals() need not be
278 % implemented for grid-based methods, as the normals simply can be
279 % obtained by the grid. The methods are only required, if using
280 % meshless methods with data functions that use normals.
281 if model.all_dirichlet_boundary == 1
282  model.boundary_type = @all_dirichlet_boundary_type;
283 else
284  model.boundary_type = @mixed_boundary_type;
285 end;
286 model.normals = @my_normals;
287 
288 % additional settings for reduced basis approach:
289 model.mu_names = {'mu_a','mu_b','mu_c'};
290 model.mu_ranges = {[-1,1],[0,1],[0,1]};
291 
292 model.coercivity_alpha = @(model) 1;
293 model.enable_error_estimator = 0;
294 
295 %%%%%%% auxiliary functions:
296 
297 % all edges of unit square are dirichlet, other inner
298 function res = all_dirichlet_boundary_type(glob,params)
299 res = zeros(size(glob,1),1);
300 i = find(glob(:,1)<=1e-10);
301 i = [i, find(glob(:,1)>=1-1e-10)];
302 i = [i, find(glob(:,2)<=1e-10)];
303 i = [i, find(glob(:,2)>=1-1e-10)];
304 res(i) = -1;
305 
306 % right and upper edges of unit square are dirichlet,
307 % left neumann, lower robin
308 function res = mixed_boundary_type(glob,params)
309 res = zeros(size(glob,1),1);
310 i = find(glob(:,1)<=1e-10);
311 res(i) = -2;
312 i = find(glob(:,2)<=1e-10);
313 res(i) = -3;
314 i = find(glob(:,1)>=1-1e-10);
315 res(i) = -1;
316 i = find(glob(:,2)>=1-1e-10);
317 res(i) = -1;
318 
319 function res = my_normals(glob,params);
320 res = zeros(size(glob,1),2); % each row one normal
321 i = find(glob(:,1)>1-1e-10);
322 res(i,1)= 1.0;
323 i = find(glob(:,1)<1e-10);
324 res(i,1)= -1.0;
325 i = find(glob(:,2)>1-1e-10);
326 res(i,2)= 1.0;
327 i = find(glob(:,2)<1e-10);
328 res(i,2)= -1.0;
329 % remove diagonal normals
330 i = find (sum(abs(res),2)>1.5);
331 res(i,1)= 0;
332 
333 % neumann boundary value components:
334 % multiplication of diffusivity components with solution gradient
335 function res = my_neumann_values_components(glob,params);
336 diff_comp = params.diffusivity_tensor(glob,params);
337 u_grad = params.solution_gradient(glob,params);
338 normals = params.normals(glob,params);
339 res = cell(1,length(diff_comp));
340 for q = 1:length(diff_comp)
341  % Agradu: each row one Agradu value
342  if ~isempty(diff_comp{q})
343  Agradu = [diff_comp{q}(:,1).* u_grad(:,1) + ...
344  diff_comp{q}(:,3).* u_grad(:,2), ...
345  diff_comp{q}(:,2).* u_grad(:,1) + ...
346  diff_comp{q}(:,4).* u_grad(:,2)];
347  res{q} = Agradu(:,1).* normals(:,1) + Agradu(:,2).* normals(:,2);
348  else
349  res{q} = zeros(size(glob,1),1);
350  end;
351 end;
352 
353 function res = my_robin_values_coefficients(dummy,params)
354 res = [params.dirichlet_values(dummy,params); ...
355  params.neumann_values(dummy,params)];
356 
357 function res = my_robin_values_components(glob,params)
358 dir_comp = params.dirichlet_values(glob,params);
359 neu_comp = params.neumann_values(glob,params);
360 robin_alpha = params.robin_alpha(glob,params);
361 robin_beta = params.robin_beta(glob,params);
362 for q = 1:length(dir_comp)
363  dir_comp{q} = dir_comp{q}.*robin_alpha;
364 end;
365 for q = 1:length(neu_comp)
366  neu_comp{q} = neu_comp{q}.*robin_beta;
367 end;
368 res = [dir_comp, neu_comp];
369 
370 % source term: diffusivity, then advection, then reaction contributions
371 function res = my_source_coefficients(dummy,params);
372 res = [-params.diffusivity_tensor(dummy,params);...
373  params.velocity(dummy,params);...
374  params.reaction(dummy,params)];
375 
376 % source term: diffusivity, then advection, then reaction
377 % contribution. Assume that a(x) has vanishing derivative,
378 % i.e piecewise constant!!!
379 % velocity assumed to be divergence free!
380 function res = my_source_components(glob,params);
381 u = params.solution(glob,params);
382 u_grad = params.solution_gradient(glob,params);
383 % hessian: each row: u_xx u_yx, u_xy, u_yy
384 u_hessian = params.solution_hessian(glob,params);
385 diff_comp = params.diffusivity_tensor(glob,params);
386 vel_comp = params.velocity(glob,params);
387 reac_comp = params.reaction(glob,params);
388 source_diff_comp = cell(1,length(diff_comp));
389 source_adv_comp = cell(1,length(vel_comp));
390 source_reac_comp = cell(1,length(reac_comp));
391 for q = 1:length(reac_comp)
392  source_reac_comp{q} = reac_comp{q}.* u;
393 end;
394 for q = 1:length(vel_comp)
395  source_adv_comp{q} = vel_comp{q}(:,1).* u_grad(:,1) + ...
396  vel_comp{q}(:,2).* u_grad(:,2);
397 end;
398 for q = 1:length(diff_comp)
399  source_diff_comp{q} = ...
400  diff_comp{q}(:,1).* u_hessian(:,1) + ...
401  diff_comp{q}(:,2).* u_hessian(:,2) + ...
402  diff_comp{q}(:,3).* u_hessian(:,3) + ...
403  diff_comp{q}(:,4).* u_hessian(:,4);
404 end;
405 res = [source_diff_comp, source_adv_comp, source_reac_comp];
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]))