Mercurial > repos > public > sbplib

comparison +scheme/LaplaceCurvilinearNew.m @ 1331:60c875c18de3 feature/D2_boundary_opt

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with feature/poroelastic for Elastic schemes
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Thu, 10 Mar 2022 16:54:26 +0100
parents 5c5815af4b7a
children
comparison
equal deleted inserted replaced
1330:855871e0b852 1331:60c875c18de3
1 classdef LaplaceCurvilinearNew < scheme.Scheme
2 properties
3 m % Number of points in each direction, possibly a vector
4 h % Grid spacing
5 dim % Number of spatial dimensions
6
7 grid
8
9 order % Order of accuracy for the approximation
10
11 a,b % Parameters of the operator
12 weight % Parameter in front of time derivative (e.g. u_tt in wave equation) here: 1/a.
13
14
15 % Inner products and operators for physical coordinates
16 D % Laplace operator
17 H, Hi % Inner product
18 e_w, e_e, e_s, e_n
19 d_w, d_e, d_s, d_n % Normal derivatives at the boundary
20 H_w, H_e, H_s, H_n % Boundary inner products
21 Dx, Dy % Physical derivatives
22 M % Gradient inner product
23
24 % Metric coefficients
25 J, Ji
26 a11, a12, a22
27 K
28 x_u
29 x_v
30 y_u
31 y_v
32 s_w, s_e, s_s, s_n % Boundary integral scale factors
33
34 % Inner product and operators for logical coordinates
35 H_u, H_v % Norms in the x and y directions
36 Hi_u, Hi_v
37 Hu,Hv % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir.
38 Hiu, Hiv
39 du_w, dv_w
40 du_e, dv_e
41 du_s, dv_s
42 du_n, dv_n
43
44 % Borrowing constants
45 theta_M_u, theta_M_v
46 theta_R_u, theta_R_v
47 theta_H_u, theta_H_v
48
49 % Temporary, only used for nonconforming interfaces but should be removed.
50 lambda
51 end
52
53 methods
54 % Implements a*div(b*grad(u)) as a SBP scheme
55 % TODO: Implement proper H, it should be the real physical quadrature, the logic quadrature may be but in a separate variable (H_logic?)
56
57 function obj = LaplaceCurvilinearNew(g, order, a, b, opSet)
58 default_arg('opSet',@sbp.D2Variable);
59 default_arg('a', 1);
60 default_arg('b', @(x,y) 0*x + 1);
61
62 % assert(isa(g, 'grid.Curvilinear'))
63 if isa(a, 'function_handle')
64 a = grid.evalOn(g, a);
65 end
66
67 % If a is scalar
68 if length(a) == 1
69 a = a*ones(g.N(), 1);
70 end
71 a = spdiag(a);
72
73 if isa(b, 'function_handle')
74 b = grid.evalOn(g, b);
75 end
76
77 % If b is scalar
78 if length(b) == 1
79 % b = b*speye(g.N(), g.N());
80 b = b*ones(g.N(), 1);
81 end
82 b = spdiag(b);
83
84 dim = 2;
85 m = g.size();
86 m_u = m(1);
87 m_v = m(2);
88 m_tot = g.N();
89
90 % 1D operators
91 ops_u = opSet(m_u, {0, 1}, order);
92 ops_v = opSet(m_v, {0, 1}, order);
93
94 h_u = ops_u.h;
95 h_v = ops_v.h;
96
97 I_u = speye(m_u);
98 I_v = speye(m_v);
99
100 D1_u = ops_u.D1;
101 D2_u = ops_u.D2;
102 H_u = ops_u.H;
103 Hi_u = ops_u.HI;
104 e_l_u = ops_u.e_l;
105 e_r_u = ops_u.e_r;
106 d1_l_u = ops_u.d1_l;
107 d1_r_u = ops_u.d1_r;
108
109 D1_v = ops_v.D1;
110 D2_v = ops_v.D2;
111 H_v = ops_v.H;
112 Hi_v = ops_v.HI;
113 e_l_v = ops_v.e_l;
114 e_r_v = ops_v.e_r;
115 d1_l_v = ops_v.d1_l;
116 d1_r_v = ops_v.d1_r;
117
118
119 % Logical operators
120 Du = kr(D1_u,I_v);
121 Dv = kr(I_u,D1_v);
122 obj.Hu = kr(H_u,I_v);
123 obj.Hv = kr(I_u,H_v);
124 obj.Hiu = kr(Hi_u,I_v);
125 obj.Hiv = kr(I_u,Hi_v);
126
127 e_w = kr(e_l_u,I_v);
128 e_e = kr(e_r_u,I_v);
129 e_s = kr(I_u,e_l_v);
130 e_n = kr(I_u,e_r_v);
131 obj.du_w = kr(d1_l_u,I_v);
132 obj.dv_w = (e_w'*Dv)';
133 obj.du_e = kr(d1_r_u,I_v);
134 obj.dv_e = (e_e'*Dv)';
135 obj.du_s = (e_s'*Du)';
136 obj.dv_s = kr(I_u,d1_l_v);
137 obj.du_n = (e_n'*Du)';
138 obj.dv_n = kr(I_u,d1_r_v);
139
140
141 % Metric coefficients
142 coords = g.points();
143 x = coords(:,1);
144 y = coords(:,2);
145
146 x_u = Du*x;
147 x_v = Dv*x;
148 y_u = Du*y;
149 y_v = Dv*y;
150
151 J = x_u.*y_v - x_v.*y_u;
152 a11 = 1./J .* (x_v.^2 + y_v.^2);
153 a12 = -1./J .* (x_u.*x_v + y_u.*y_v);
154 a22 = 1./J .* (x_u.^2 + y_u.^2);
155 lambda = 1/2 * (a11 + a22 - sqrt((a11-a22).^2 + 4*a12.^2));
156
157 K = cell(dim, dim);
158 K{1,1} = spdiag(y_v./J);
159 K{1,2} = spdiag(-y_u./J);
160 K{2,1} = spdiag(-x_v./J);
161 K{2,2} = spdiag(x_u./J);
162 obj.K = K;
163
164 obj.x_u = x_u;
165 obj.x_v = x_v;
166 obj.y_u = y_u;
167 obj.y_v = y_v;
168
169 % Assemble full operators
170 L_12 = spdiag(a12);
171 Duv = Du*b*L_12*Dv;
172 Dvu = Dv*b*L_12*Du;
173
174 Duu = sparse(m_tot);
175 Dvv = sparse(m_tot);
176 ind = grid.funcToMatrix(g, 1:m_tot);
177
178 for i = 1:m_v
179 b_a11 = b*a11;
180 D = D2_u(b_a11(ind(:,i)));
181 p = ind(:,i);
182 Duu(p,p) = D;
183 end
184
185 for i = 1:m_u
186 b_a22 = b*a22;
187 D = D2_v(b_a22(ind(i,:)));
188 p = ind(i,:);
189 Dvv(p,p) = D;
190 end
191
192
193 % Physical operators
194 obj.J = spdiag(J);
195 obj.Ji = spdiag(1./J);
196
197 obj.D = obj.Ji*a*(Duu + Duv + Dvu + Dvv);
198 obj.H = obj.J*kr(H_u,H_v);
199 obj.Hi = obj.Ji*kr(Hi_u,Hi_v);
200
201 obj.e_w = e_w;
202 obj.e_e = e_e;
203 obj.e_s = e_s;
204 obj.e_n = e_n;
205
206 %% normal derivatives
207 I_w = ind(1,:);
208 I_e = ind(end,:);
209 I_s = ind(:,1);
210 I_n = ind(:,end);
211
212 a11_w = spdiag(a11(I_w));
213 a12_w = spdiag(a12(I_w));
214 a11_e = spdiag(a11(I_e));
215 a12_e = spdiag(a12(I_e));
216 a22_s = spdiag(a22(I_s));
217 a12_s = spdiag(a12(I_s));
218 a22_n = spdiag(a22(I_n));
219 a12_n = spdiag(a12(I_n));
220
221 s_w = sqrt((e_w'*x_v).^2 + (e_w'*y_v).^2);
222 s_e = sqrt((e_e'*x_v).^2 + (e_e'*y_v).^2);
223 s_s = sqrt((e_s'*x_u).^2 + (e_s'*y_u).^2);
224 s_n = sqrt((e_n'*x_u).^2 + (e_n'*y_u).^2);
225
226 obj.d_w = -1*(spdiag(1./s_w)*(a11_w*obj.du_w' + a12_w*obj.dv_w'))';
227 obj.d_e = (spdiag(1./s_e)*(a11_e*obj.du_e' + a12_e*obj.dv_e'))';
228 obj.d_s = -1*(spdiag(1./s_s)*(a22_s*obj.dv_s' + a12_s*obj.du_s'))';
229 obj.d_n = (spdiag(1./s_n)*(a22_n*obj.dv_n' + a12_n*obj.du_n'))';
230
231 obj.Dx = spdiag( y_v./J)*Du + spdiag(-y_u./J)*Dv;
232 obj.Dy = spdiag(-x_v./J)*Du + spdiag( x_u./J)*Dv;
233
234 %% Boundary inner products
235 obj.H_w = H_v*spdiag(s_w);
236 obj.H_e = H_v*spdiag(s_e);
237 obj.H_s = H_u*spdiag(s_s);
238 obj.H_n = H_u*spdiag(s_n);
239
240 % Misc.
241 obj.m = m;
242 obj.h = [h_u h_v];
243 obj.order = order;
244 obj.grid = g;
245 obj.dim = dim;
246
247 obj.a = a;
248 obj.weight = inv(a);
249 obj.b = b;
250 obj.a11 = a11;
251 obj.a12 = a12;
252 obj.a22 = a22;
253 obj.s_w = spdiag(s_w);
254 obj.s_e = spdiag(s_e);
255 obj.s_s = spdiag(s_s);
256 obj.s_n = spdiag(s_n);
257
258 obj.theta_M_u = h_u*ops_u.borrowing.M.d1;
259 obj.theta_M_v = h_v*ops_v.borrowing.M.d1;
260
261 obj.theta_R_u = h_u*ops_u.borrowing.R.delta_D;
262 obj.theta_R_v = h_v*ops_v.borrowing.R.delta_D;
263
264 obj.theta_H_u = h_u*ops_u.borrowing.H11;
265 obj.theta_H_v = h_v*ops_v.borrowing.H11;
266
267 % Temporary
268 obj.lambda = lambda;
269 end
270
271
272 % Closure functions return the opertors applied to the own doamin to close the boundary
273 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin.
274 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
275 % type is a string specifying the type of boundary condition if there are several.
276 % data is a function returning the data that should be applied at the boundary.
277 % neighbour_scheme is an instance of Scheme that should be interfaced to.
278 % neighbour_boundary is a string specifying which boundary to interface to.
279 function [closure, penalty] = boundary_condition(obj, boundary, type, parameter)
280 default_arg('type','neumann');
281 default_arg('parameter', []);
282
283 e = obj.getBoundaryOperator('e', boundary);
284 d = obj.getBoundaryOperator('d', boundary);
285 H_b = obj.getBoundaryQuadrature(boundary);
286 s_b = obj.getBoundaryScaling(boundary);
287 [th_H, ~, th_R] = obj.getBoundaryBorrowing(boundary);
288 m = obj.getBoundaryNumber(boundary);
289
290 K = obj.K;
291 J = obj.J;
292 Hi = obj.Hi;
293 a = obj.a;
294 b_b = e'*obj.b*e;
295
296 switch type
297 % Dirichlet boundary condition
298 case {'D','d','dirichlet'}
299 tuning = 1.0;
300
301 sigma = 0*b_b;
302 for i = 1:obj.dim
303 sigma = sigma + e'*J*K{i,m}*K{i,m}*e;
304 end
305 sigma = sigma/s_b;
306 tau = tuning*(1/th_R + obj.dim/th_H)*sigma;
307
308 closure = a*Hi*d*b_b*H_b*e' ...
309 -a*Hi*e*tau*b_b*H_b*e';
310
311 penalty = -a*Hi*d*b_b*H_b ...
312 +a*Hi*e*tau*b_b*H_b;
313
314
315 % Neumann boundary condition. Note that the penalty is for du/dn and not b*du/dn.
316 case {'N','n','neumann'}
317 tau1 = -1;
318 tau2 = 0;
319 tau = (tau1*e + tau2*d)*H_b;
320
321 closure = a*Hi*tau*b_b*d';
322 penalty = -a*Hi*tau*b_b;
323
324
325 % Unknown, boundary condition
326 otherwise
327 error('No such boundary condition: type = %s',type);
328 end
329 end
330
331 % type Struct that specifies the interface coupling.
332 % Fields:
333 % -- tuning: penalty strength, defaults to 1.2
334 % -- interpolation: type of interpolation, default 'none'
335 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type)
336
337 defaultType.coupling = 'sat';
338 defaultType.tuning = 1.0;
339 defaultType.interpolation = 'none';
340 default_struct('type', defaultType);
341
342 switch type.coupling
343 case {'cg', 'CG'}
344 [closure, penalty] = interfaceCG(obj,boundary,neighbour_scheme,neighbour_boundary,type);
345 case {'sat', 'SAT'}
346 switch type.interpolation
347 case {'none', ''}
348 [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type);
349 case {'op','OP'}
350 [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type);
351 otherwise
352 error('Unknown type of interpolation: %s ', type.interpolation);
353 end
354 otherwise
355 error('Unknown type of coupling: %s ', type.coupling);
356 end
357 end
358
359 function [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type)
360 tuning = type.tuning;
361
362 dim = obj.dim;
363 % u denotes the solution in the own domain
364 % v denotes the solution in the neighbour domain
365 u = obj;
366 v = neighbour_scheme;
367
368 % Boundary operators, u
369 e_u = u.getBoundaryOperator('e', boundary);
370 d_u = u.getBoundaryOperator('d', boundary);
371 s_b_u = u.getBoundaryScaling(boundary);
372 [th_H_u, ~, th_R_u] = u.getBoundaryBorrowing(boundary);
373 m_u = u.getBoundaryNumber(boundary);
374
375 % Coefficients, u
376 K_u = u.K;
377 J_u = u.J;
378 b_b_u = e_u'*u.b*e_u;
379
380 % Boundary operators, v
381 e_v = v.getBoundaryOperator('e', neighbour_boundary);
382 d_v = v.getBoundaryOperator('d', neighbour_boundary);
383 s_b_v = v.getBoundaryScaling(neighbour_boundary);
384 [th_H_v, ~, th_R_v] = v.getBoundaryBorrowing(neighbour_boundary);
385 m_v = v.getBoundaryNumber(neighbour_boundary);
386
387 % BUGFIX?!?!?
388 if (strcmp(boundary,'s') && strcmp(neighbour_boundary,'e')) || (strcmp(boundary,'e') && strcmp(neighbour_boundary,'s'))
389 e_v = fliplr(e_v);
390 d_v = fliplr(d_v);
391 s_b_v = rot90(s_b_v,2);
392 end
393
394 % Coefficients, v
395 K_v = v.K;
396 J_v = v.J;
397 b_b_v = e_v'*v.b*e_v;
398
399 %--- Penalty strength tau -------------
400 sigma_u = 0*b_b_u;
401 sigma_v = 0*b_b_v;
402 for i = 1:obj.dim
403 sigma_u = sigma_u + e_u'*J_u*K_u{i,m_u}*K_u{i,m_u}*e_u;
404 sigma_v = sigma_v + e_v'*J_v*K_v{i,m_v}*K_v{i,m_v}*e_v;
405 end
406 sigma_u = sigma_u/s_b_u;
407 sigma_v = sigma_v/s_b_v;
408
409 tau_R_u = 1/th_R_u*sigma_u;
410 tau_R_v = 1/th_R_v*sigma_v;
411
412 tau_H_u = dim*1/th_H_u*sigma_u;
413 tau_H_v = dim*1/th_H_v*sigma_v;
414
415 tau = 1/4*tuning*(b_b_u*(tau_R_u + tau_H_u) + b_b_v*(tau_R_v + tau_H_v));
416 %--------------------------------------
417
418 % Operators/coefficients that are only required from this side
419 Hi = u.Hi;
420 H_b = u.getBoundaryQuadrature(boundary);
421 a = u.a;
422
423 closure = 1/2*a*Hi*d_u*b_b_u*H_b*e_u' ...
424 -1/2*a*Hi*e_u*H_b*b_b_u*d_u' ...
425 -a*Hi*e_u*tau*H_b*e_u';
426
427 penalty = -1/2*a*Hi*d_u*b_b_u*H_b*e_v' ...
428 -1/2*a*Hi*e_u*H_b*b_b_v*d_v' ...
429 +a*Hi*e_u*tau*H_b*e_v';
430 end
431
432 function [closure, penalty] = interfaceCG(obj,boundary,neighbour_scheme,neighbour_boundary,type)
433
434 % There is no penalty, only a closure. And the closure is the same as for Neumann BC
435 e = obj.getBoundaryOperator('e', boundary);
436 d = obj.getBoundaryOperator('d', boundary);
437 H_b = obj.getBoundaryQuadrature(boundary);
438
439 e_v = neighbour_scheme.getBoundaryOperator('e', neighbour_boundary);
440
441 Hi = obj.Hi;
442 a = obj.a;
443 b_b = e'*obj.b*e;
444
445 tau1 = -1;
446 tau2 = 0;
447 tau = (tau1*e + tau2*d)*H_b;
448
449 closure = a*Hi*tau*b_b*d';
450
451 % Zero penalty of correct dimensions
452 penalty = 0*e*e_v';
453 end
454
455 function [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type)
456
457 % TODO: Make this work for curvilinear grids
458 warning('LaplaceCurvilinear: Non-conforming grid interpolation only works for Cartesian grids.');
459 warning('LaplaceCurvilinear: Non-conforming interface uses Virtas penalty strength');
460 warning('LaplaceCurvilinear: Non-conforming interface assumes that b is constant');
461
462 % User can request special interpolation operators by specifying type.interpOpSet
463 default_field(type, 'interpOpSet', @sbp.InterpOpsOP);
464 interpOpSet = type.interpOpSet;
465 tuning = type.tuning;
466
467
468 % u denotes the solution in the own domain
469 % v denotes the solution in the neighbour domain
470 e_u = obj.getBoundaryOperator('e', boundary);
471 d_u = obj.getBoundaryOperator('d', boundary);
472 H_b_u = obj.getBoundaryQuadrature(boundary);
473 I_u = obj.getBoundaryIndices(boundary);
474 [~, gamm_u] = obj.getBoundaryBorrowing(boundary);
475
476 e_v = neighbour_scheme.getBoundaryOperator('e', neighbour_boundary);
477 d_v = neighbour_scheme.getBoundaryOperator('d', neighbour_boundary);
478 H_b_v = neighbour_scheme.getBoundaryQuadrature(neighbour_boundary);
479 I_v = neighbour_scheme.getBoundaryIndices(neighbour_boundary);
480 [~, gamm_v] = neighbour_scheme.getBoundaryBorrowing(neighbour_boundary);
481
482
483 % Find the number of grid points along the interface
484 m_u = size(e_u, 2);
485 m_v = size(e_v, 2);
486
487 Hi = obj.Hi;
488 a = obj.a;
489
490 u = obj;
491 v = neighbour_scheme;
492
493 b1_u = gamm_u*u.lambda(I_u)./u.a11(I_u).^2;
494 b2_u = gamm_u*u.lambda(I_u)./u.a22(I_u).^2;
495 b1_v = gamm_v*v.lambda(I_v)./v.a11(I_v).^2;
496 b2_v = gamm_v*v.lambda(I_v)./v.a22(I_v).^2;
497
498 tau_u = -1./(4*b1_u) -1./(4*b2_u);
499 tau_v = -1./(4*b1_v) -1./(4*b2_v);
500
501 tau_u = tuning * spdiag(tau_u);
502 tau_v = tuning * spdiag(tau_v);
503 beta_u = tau_v;
504
505 % Build interpolation operators
506 intOps = interpOpSet(m_u, m_v, obj.order, neighbour_scheme.order);
507 Iu2v = intOps.Iu2v;
508 Iv2u = intOps.Iv2u;
509
510 closure = a*Hi*e_u*tau_u*H_b_u*e_u' + ...
511 a*Hi*e_u*H_b_u*Iv2u.bad*beta_u*Iu2v.good*e_u' + ...
512 a*1/2*Hi*d_u*H_b_u*e_u' + ...
513 -a*1/2*Hi*e_u*H_b_u*d_u';
514
515 penalty = -a*Hi*e_u*tau_u*H_b_u*Iv2u.good*e_v' + ...
516 -a*Hi*e_u*H_b_u*Iv2u.bad*beta_u*e_v' + ...
517 -a*1/2*Hi*d_u*H_b_u*Iv2u.good*e_v' + ...
518 -a*1/2*Hi*e_u*H_b_u*Iv2u.bad*d_v';
519
520 end
521
522 % Returns the boundary operator op for the boundary specified by the string boundary.
523 % op -- string
524 % boundary -- string
525 function o = getBoundaryOperator(obj, op, boundary)
526 assertIsMember(op, {'e', 'd'})
527 assertIsMember(boundary, {'w', 'e', 's', 'n'})
528
529 o = obj.([op, '_', boundary]);
530 end
531
532 % Returns square boundary quadrature matrix, of dimension
533 % corresponding to the number of boundary points
534 %
535 % boundary -- string
536 function H_b = getBoundaryQuadrature(obj, boundary)
537 assertIsMember(boundary, {'w', 'e', 's', 'n'})
538
539 H_b = obj.(['H_', boundary]);
540 end
541
542 % Returns square boundary quadrature scaling matrix, of dimension
543 % corresponding to the number of boundary points
544 %
545 % boundary -- string
546 function s_b = getBoundaryScaling(obj, boundary)
547 assertIsMember(boundary, {'w', 'e', 's', 'n'})
548
549 s_b = obj.(['s_', boundary]);
550 end
551
552 % Returns the coordinate number corresponding to the boundary
553 %
554 % boundary -- string
555 function m = getBoundaryNumber(obj, boundary)
556 assertIsMember(boundary, {'w', 'e', 's', 'n'})
557
558 switch boundary
559 case {'w', 'e'}
560 m = 1;
561 case {'s', 'n'}
562 m = 2;
563 end
564 end
565
566 % Returns the indices of the boundary points in the grid matrix
567 % boundary -- string
568 function I = getBoundaryIndices(obj, boundary)
569 assertIsMember(boundary, {'w', 'e', 's', 'n'})
570
571 ind = grid.funcToMatrix(obj.grid, 1:prod(obj.m));
572 switch boundary
573 case 'w'
574 I = ind(1,:);
575 case 'e'
576 I = ind(end,:);
577 case 's'
578 I = ind(:,1)';
579 case 'n'
580 I = ind(:,end)';
581 end
582 end
583
584 % Returns borrowing constant gamma
585 % boundary -- string
586 function [theta_H, theta_M, theta_R] = getBoundaryBorrowing(obj, boundary)
587 assertIsMember(boundary, {'w', 'e', 's', 'n'})
588
589 switch boundary
590 case {'w','e'}
591 theta_H = obj.theta_H_u;
592 theta_M = obj.theta_M_u;
593 theta_R = obj.theta_R_u;
594 case {'s','n'}
595 theta_H = obj.theta_H_v;
596 theta_M = obj.theta_M_v;
597 theta_R = obj.theta_R_v;
598 end
599 end
600
601 function N = size(obj)
602 N = prod(obj.m);
603 end
604 end
605 end