Mercurial > repos > public > sbplib

comparison +scheme/LaplaceCurvilinearNewCorner.m @ 1055:b828e589d540 feature/poroelastic

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