Mercurial > repos > public > sbplib

comparison +scheme/Elastic2dVariableAnisotropicUpwind.m @ 1227:02dfe3a56742 feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add Upwind ElasticAnisotropic schemes. Seem to work really well!
author Martin Almquist <malmquist@stanford.edu>
date Sat, 16 Nov 2019 14:26:06 -0800
parents
children b5025bd67be1
comparison
equal deleted inserted replaced
1214:e1d4cb8b5309 1227:02dfe3a56742
1 classdef Elastic2dVariableAnisotropicUpwind < scheme.Scheme
2
3 % Discretizes the elastic wave equation:
4 % rho u_{i,tt} = dj C_{ijkl} dk u_j
5 % opSet should be cell array of opSets, one per dimension. This
6 % is useful if we have periodic BC in one direction.
7 % Assumes fully compatible operators
8
9 properties
10 m % Number of points in each direction, possibly a vector
11 h % Grid spacing
12
13 grid
14 dim
15
16 order % Order of accuracy for the approximation
17
18 % Diagonal matrices for variable coefficients
19 RHO, RHOi, RHOi_kron % Density
20 C % Elastic stiffness tensor
21
22 D % Total operator
23 Dp, Dm % First derivatives
24
25 % Boundary operators in cell format, used for BC
26 T_w, T_e, T_s, T_n
27
28 % Traction operators
29 tau_w, tau_e, tau_s, tau_n % Return vector field
30 tau1_w, tau1_e, tau1_s, tau1_n % Return scalar field
31 tau2_w, tau2_e, tau2_s, tau2_n % Return scalar field
32
33 % Inner products
34 H, Hi, Hi_kron, H_1D
35
36 % Boundary inner products (for scalar field)
37 H_w, H_e, H_s, H_n
38
39 % Boundary restriction operators
40 e_w, e_e, e_s, e_n % Act on vector field, return vector field at boundary
41 e1_w, e1_e, e1_s, e1_n % Act on vector field, return scalar field at boundary
42 e2_w, e2_e, e2_s, e2_n % Act on vector field, return scalar field at boundary
43 e_scalar_w, e_scalar_e, e_scalar_s, e_scalar_n; % Act on scalar field, return scalar field
44
45 % E{i}^T picks out component i
46 E
47
48 % Borrowing constants of the form gamma*h, where gamma is a dimensionless constant.
49 h11 % First entry in norm matrix
50
51 end
52
53 methods
54
55 % The coefficients can either be function handles or grid functions
56 % optFlag -- if true, extra computations are performed, which may be helpful for optimization.
57 function obj = Elastic2dVariableAnisotropicUpwind(g, order, rho, C, opSet, optFlag)
58 default_arg('rho', @(x,y) 0*x+1);
59 default_arg('opSet',{@sbp.D1Upwind, @sbp.D1Upwind});
60 default_arg('optFlag', false);
61 dim = 2;
62
63 C_default = cell(dim,dim,dim,dim);
64 for i = 1:dim
65 for j = 1:dim
66 for k = 1:dim
67 for l = 1:dim
68 C_default{i,j,k,l} = @(x,y) 0*x + 1;
69 end
70 end
71 end
72 end
73 default_arg('C', C_default);
74 assert(isa(g, 'grid.Cartesian'))
75
76 if isa(rho, 'function_handle')
77 rho = grid.evalOn(g, rho);
78 end
79
80 C_mat = cell(dim,dim,dim,dim);
81 for i = 1:dim
82 for j = 1:dim
83 for k = 1:dim
84 for l = 1:dim
85 if isa(C{i,j,k,l}, 'function_handle')
86 C{i,j,k,l} = grid.evalOn(g, C{i,j,k,l});
87 end
88 C_mat{i,j,k,l} = spdiag(C{i,j,k,l});
89 end
90 end
91 end
92 end
93 obj.C = C_mat;
94
95 m = g.size();
96 m_tot = g.N();
97 lim = g.lim;
98 if isempty(lim)
99 x = g.x;
100 lim = cell(length(x),1);
101 for i = 1:length(x)
102 lim{i} = {min(x{i}), max(x{i})};
103 end
104 end
105
106 % 1D operators
107 ops = cell(dim,1);
108 h = zeros(dim,1);
109 for i = 1:dim
110 ops{i} = opSet{i}(m(i), lim{i}, order);
111 h(i) = ops{i}.h;
112 end
113
114 % Borrowing constants
115 for i = 1:dim
116 obj.h11{i} = ops{i}.H(1,1);
117 end
118
119 I = cell(dim,1);
120 Dp = cell(dim,1);
121 Dm = cell(dim,1);
122 H = cell(dim,1);
123 Hi = cell(dim,1);
124 e_0 = cell(dim,1);
125 e_m = cell(dim,1);
126 d1_0 = cell(dim,1);
127 d1_m = cell(dim,1);
128
129 for i = 1:dim
130 I{i} = speye(m(i));
131 Dp{i} = ops{i}.Dp;
132 Dm{i} = ops{i}.Dm;
133 H{i} = ops{i}.H;
134 Hi{i} = ops{i}.HI;
135 e_0{i} = ops{i}.e_l;
136 e_m{i} = ops{i}.e_r;
137 d1_0{i} = (ops{i}.e_l' * Dm{i})';
138 d1_m{i} = (ops{i}.e_r' * Dm{i})';
139 end
140
141 %====== Assemble full operators ========
142 I_dim = speye(dim, dim);
143 RHO = spdiag(rho);
144 obj.RHO = RHO;
145 obj.RHOi = inv(RHO);
146 obj.RHOi_kron = kron(obj.RHOi, I_dim);
147
148 obj.Dp = cell(dim,1);
149 obj.Dm = cell(dim,1);
150
151 % D1
152 obj.Dp{1} = kron(Dp{1},I{2});
153 obj.Dp{2} = kron(I{1},Dp{2});
154 obj.Dm{1} = kron(Dm{1},I{2});
155 obj.Dm{2} = kron(I{1},Dm{2});
156
157 % Boundary restriction operators
158 e_l = cell(dim,1);
159 e_r = cell(dim,1);
160 e_l{1} = kron(e_0{1}, I{2});
161 e_l{2} = kron(I{1}, e_0{2});
162 e_r{1} = kron(e_m{1}, I{2});
163 e_r{2} = kron(I{1}, e_m{2});
164
165 e_scalar_w = e_l{1};
166 e_scalar_e = e_r{1};
167 e_scalar_s = e_l{2};
168 e_scalar_n = e_r{2};
169
170 e_w = kron(e_scalar_w, I_dim);
171 e_e = kron(e_scalar_e, I_dim);
172 e_s = kron(e_scalar_s, I_dim);
173 e_n = kron(e_scalar_n, I_dim);
174
175 % E{i}^T picks out component i.
176 E = cell(dim,1);
177 I = speye(m_tot,m_tot);
178 for i = 1:dim
179 e = sparse(dim,1);
180 e(i) = 1;
181 E{i} = kron(I,e);
182 end
183 obj.E = E;
184
185 e1_w = (e_scalar_w'*E{1}')';
186 e1_e = (e_scalar_e'*E{1}')';
187 e1_s = (e_scalar_s'*E{1}')';
188 e1_n = (e_scalar_n'*E{1}')';
189
190 e2_w = (e_scalar_w'*E{2}')';
191 e2_e = (e_scalar_e'*E{2}')';
192 e2_s = (e_scalar_s'*E{2}')';
193 e2_n = (e_scalar_n'*E{2}')';
194
195 % Quadratures
196 obj.H = kron(H{1},H{2});
197 obj.Hi = inv(obj.H);
198 obj.H_w = H{2};
199 obj.H_e = H{2};
200 obj.H_s = H{1};
201 obj.H_n = H{1};
202 obj.H_1D = {H{1}, H{2}};
203
204 % Differentiation matrix D (without SAT)
205 Dp = obj.Dp;
206 Dm = obj.Dm;
207 D = sparse(dim*m_tot,dim*m_tot);
208 for i = 1:dim
209 for j = 1:dim
210 for k = 1:dim
211 for l = 1:dim
212 D = D + E{j}*Dp{i}*C_mat{i,j,k,l}*Dm{k}*E{l}';
213 end
214 end
215 end
216 end
217 D = obj.RHOi_kron*D;
218 obj.D = D;
219 %=========================================%'
220
221 % Numerical traction operators for BC.
222 %
223 % Formula at boundary j: % tau^{j}_i = sum_l T^{j}_{il} u_l
224 %
225 T_l = cell(dim,1);
226 T_r = cell(dim,1);
227 tau_l = cell(dim,1);
228 tau_r = cell(dim,1);
229
230 D1 = obj.Dm;
231
232 % Boundary j
233 for j = 1:dim
234 T_l{j} = cell(dim,dim);
235 T_r{j} = cell(dim,dim);
236 tau_l{j} = cell(dim,1);
237 tau_r{j} = cell(dim,1);
238
239 [~, n_l] = size(e_l{j});
240 [~, n_r] = size(e_r{j});
241
242 % Traction component i
243 for i = 1:dim
244 tau_l{j}{i} = sparse(dim*m_tot, n_l);
245 tau_r{j}{i} = sparse(dim*m_tot, n_r);
246
247 % Displacement component l
248 for l = 1:dim
249 T_l{j}{i,l} = sparse(m_tot, n_l);
250 T_r{j}{i,l} = sparse(m_tot, n_r);
251
252 % Derivative direction k
253 for k = 1:dim
254 T_l{j}{i,l} = T_l{j}{i,l} ...
255 - (e_l{j}'*C_mat{j,i,k,l}*D1{k})';
256 T_r{j}{i,l} = T_r{j}{i,l} ...
257 + (e_r{j}'*C_mat{j,i,k,l}*D1{k})';
258 end
259 tau_l{j}{i} = tau_l{j}{i} + (T_l{j}{i,l}'*E{l}')';
260 tau_r{j}{i} = tau_r{j}{i} + (T_r{j}{i,l}'*E{l}')';
261 end
262 end
263 end
264
265 % Traction tensors, T_ij
266 obj.T_w = T_l{1};
267 obj.T_e = T_r{1};
268 obj.T_s = T_l{2};
269 obj.T_n = T_r{2};
270
271 % Restriction operators
272 obj.e_w = e_w;
273 obj.e_e = e_e;
274 obj.e_s = e_s;
275 obj.e_n = e_n;
276
277 obj.e1_w = e1_w;
278 obj.e1_e = e1_e;
279 obj.e1_s = e1_s;
280 obj.e1_n = e1_n;
281
282 obj.e2_w = e2_w;
283 obj.e2_e = e2_e;
284 obj.e2_s = e2_s;
285 obj.e2_n = e2_n;
286
287 obj.e_scalar_w = e_scalar_w;
288 obj.e_scalar_e = e_scalar_e;
289 obj.e_scalar_s = e_scalar_s;
290 obj.e_scalar_n = e_scalar_n;
291
292 % First component of traction
293 obj.tau1_w = tau_l{1}{1};
294 obj.tau1_e = tau_r{1}{1};
295 obj.tau1_s = tau_l{2}{1};
296 obj.tau1_n = tau_r{2}{1};
297
298 % Second component of traction
299 obj.tau2_w = tau_l{1}{2};
300 obj.tau2_e = tau_r{1}{2};
301 obj.tau2_s = tau_l{2}{2};
302 obj.tau2_n = tau_r{2}{2};
303
304 % Traction vectors
305 obj.tau_w = (e_w'*e1_w*obj.tau1_w')' + (e_w'*e2_w*obj.tau2_w')';
306 obj.tau_e = (e_e'*e1_e*obj.tau1_e')' + (e_e'*e2_e*obj.tau2_e')';
307 obj.tau_s = (e_s'*e1_s*obj.tau1_s')' + (e_s'*e2_s*obj.tau2_s')';
308 obj.tau_n = (e_n'*e1_n*obj.tau1_n')' + (e_n'*e2_n*obj.tau2_n')';
309
310 % Kroneckered norms and coefficients
311 obj.Hi_kron = kron(obj.Hi, I_dim);
312
313 % Misc.
314 obj.m = m;
315 obj.h = h;
316 obj.order = order;
317 obj.grid = g;
318 obj.dim = dim;
319
320 end
321
322
323 % Closure functions return the operators applied to the own domain to close the boundary
324 % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other doamin.
325 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
326 % bc is a cell array of component and bc type, e.g. {1, 'd'} for Dirichlet condition
327 % on the first component. Can also be e.g.
328 % {'normal', 'd'} or {'tangential', 't'} for conditions on
329 % tangential/normal component.
330 % data is a function returning the data that should be applied at the boundary.
331 % neighbour_scheme is an instance of Scheme that should be interfaced to.
332 % neighbour_boundary is a string specifying which boundary to interface to.
333
334 % For displacement bc:
335 % bc = {comp, 'd', dComps},
336 % where
337 % dComps = vector of components with displacement BC. Default: 1:dim.
338 % In this way, we can specify one BC at a time even though the SATs depend on all BC.
339 function [closure, penalty] = boundary_condition(obj, boundary, bc, tuning)
340 default_arg('tuning', 1.0);
341
342 assert( iscell(bc), 'The BC type must be a 2x1 or 3x1 cell array' );
343 comp = bc{1};
344 type = bc{2};
345 if ischar(comp)
346 comp = obj.getComponent(comp, boundary);
347 end
348
349 e = obj.getBoundaryOperatorForScalarField('e', boundary);
350 tau = obj.getBoundaryOperator(['tau' num2str(comp)], boundary);
351 T = obj.getBoundaryTractionOperator(boundary);
352 h11 = obj.getBorrowing(boundary);
353 H_gamma = obj.getBoundaryQuadratureForScalarField(boundary);
354 nu = obj.getNormal(boundary);
355
356 E = obj.E;
357 Hi = obj.Hi;
358 RHOi = obj.RHOi;
359 C = obj.C;
360
361 dim = obj.dim;
362 m_tot = obj.grid.N();
363
364 % Preallocate
365 [~, col] = size(tau);
366 closure = sparse(dim*m_tot, dim*m_tot);
367 penalty = sparse(dim*m_tot, col);
368
369 j = comp;
370 switch type
371
372 % Dirichlet boundary condition
373 case {'D','d','dirichlet','Dirichlet','displacement','Displacement'}
374
375 if numel(bc) >= 3
376 dComps = bc{3};
377 else
378 dComps = 1:dim;
379 end
380
381 % Loops over components that Dirichlet penalties end up on
382 % Y: symmetrizing part of penalty
383 % Z: symmetric part of penalty
384 % X = Y + Z.
385
386 % Nonsymmetric part goes on all components to
387 % yield traction in discrete energy rate
388 for i = 1:dim
389 Y = T{j,i}';
390 X = e*Y;
391 closure = closure + E{i}*RHOi*Hi*X'*e*H_gamma*(e'*E{j}' );
392 penalty = penalty - E{i}*RHOi*Hi*X'*e*H_gamma;
393 end
394
395 % Symmetric part only required on components with displacement BC.
396 % (Otherwise it's not symmetric.)
397 for i = dComps
398 Z = sparse(m_tot, m_tot);
399 for l = 1:dim
400 for k = 1:dim
401 Z = Z + nu(l)*C{l,i,k,j}*nu(k);
402 end
403 end
404 Z = -tuning*dim/h11*Z;
405 X = Z;
406 closure = closure + E{i}*RHOi*Hi*X'*e*H_gamma*(e'*E{j}' );
407 penalty = penalty - E{i}*RHOi*Hi*X'*e*H_gamma;
408 end
409
410 % Free boundary condition
411 case {'F','f','Free','free','traction','Traction','t','T'}
412 closure = closure - E{j}*RHOi*Hi*e*H_gamma*tau';
413 penalty = penalty + E{j}*RHOi*Hi*e*H_gamma;
414
415 % Unknown boundary condition
416 otherwise
417 error('No such boundary condition: type = %s',type);
418 end
419 end
420
421 % type Struct that specifies the interface coupling.
422 % Fields:
423 % -- tuning: penalty strength, defaults to 1.0
424 % -- interpolation: type of interpolation, default 'none'
425 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type)
426
427 defaultType.tuning = 1.0;
428 defaultType.interpolation = 'none';
429 default_struct('type', defaultType);
430
431 switch type.interpolation
432 case {'none', ''}
433 [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type);
434 case {'op','OP'}
435 [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type);
436 otherwise
437 error('Unknown type of interpolation: %s ', type.interpolation);
438 end
439 end
440
441 function [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type)
442 tuning = type.tuning;
443
444 % u denotes the solution in the own domain
445 % v denotes the solution in the neighbour domain
446
447 u = obj;
448 v = neighbour_scheme;
449
450 % Operators, u side
451 e_u = u.getBoundaryOperatorForScalarField('e', boundary);
452 tau_u = u.getBoundaryOperator('tau', boundary);
453 h11_u = u.getBorrowing(boundary);
454 nu_u = u.getNormal(boundary);
455
456 E_u = u.E;
457 C_u = u.C;
458 m_tot_u = u.grid.N();
459
460 % Operators, v side
461 e_v = v.getBoundaryOperatorForScalarField('e', neighbour_boundary);
462 tau_v = v.getBoundaryOperator('tau', neighbour_boundary);
463 h11_v = v.getBorrowing(neighbour_boundary);
464 nu_v = v.getNormal(neighbour_boundary);
465
466 E_v = v.E;
467 C_v = v.C;
468 m_tot_v = v.grid.N();
469
470 % Operators that are only required for own domain
471 Hi = u.Hi_kron;
472 RHOi = u.RHOi_kron;
473 e_kron = u.getBoundaryOperator('e', boundary);
474 T_u = u.getBoundaryTractionOperator(boundary);
475
476 % Shared operators
477 H_gamma = u.getBoundaryQuadratureForScalarField(boundary);
478 H_gamma_kron = u.getBoundaryQuadrature(boundary);
479 dim = u.dim;
480
481 % Preallocate
482 [~, m_int] = size(H_gamma);
483 closure = sparse(dim*m_tot_u, dim*m_tot_u);
484 penalty = sparse(dim*m_tot_u, dim*m_tot_v);
485
486 % ---- Continuity of displacement ------
487
488 % Y: symmetrizing part of penalty
489 % Z: symmetric part of penalty
490 % X = Y + Z.
491
492 % Loop over components to couple across interface
493 for j = 1:dim
494
495 % Loop over components that penalties end up on
496 for i = 1:dim
497 Y = 1/2*T_u{j,i}';
498 Z_u = sparse(m_int, m_int);
499 Z_v = sparse(m_int, m_int);
500 for l = 1:dim
501 for k = 1:dim
502 Z_u = Z_u + e_u'*nu_u(l)*C_u{l,i,k,j}*nu_u(k)*e_u;
503 Z_v = Z_v + e_v'*nu_v(l)*C_v{l,i,k,j}*nu_v(k)*e_v;
504 end
505 end
506 Z = -tuning*dim*( 1/(4*h11_u)*Z_u + 1/(4*h11_v)*Z_v );
507 X = Y + Z*e_u';
508 closure = closure + E_u{i}*X'*H_gamma*e_u'*E_u{j}';
509 penalty = penalty - E_u{i}*X'*H_gamma*e_v'*E_v{j}';
510 end
511 end
512
513 % ---- Continuity of traction ------
514 closure = closure - 1/2*e_kron*H_gamma_kron*tau_u';
515 penalty = penalty - 1/2*e_kron*H_gamma_kron*tau_v';
516
517 % ---- Multiply by inverse of density x quadraure ----
518 closure = RHOi*Hi*closure;
519 penalty = RHOi*Hi*penalty;
520
521 end
522
523 function [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type)
524 error('Non-conforming interfaces not implemented yet.');
525 end
526
527 % Returns the component number that is the tangential/normal component
528 % at the specified boundary
529 function comp = getComponent(obj, comp_str, boundary)
530 assertIsMember(comp_str, {'normal', 'tangential'});
531 assertIsMember(boundary, {'w', 'e', 's', 'n'});
532
533 switch boundary
534 case {'w', 'e'}
535 switch comp_str
536 case 'normal'
537 comp = 1;
538 case 'tangential'
539 comp = 2;
540 end
541 case {'s', 'n'}
542 switch comp_str
543 case 'normal'
544 comp = 2;
545 case 'tangential'
546 comp = 1;
547 end
548 end
549 end
550
551 % Returns h11 for the boundary specified by the string boundary.
552 % op -- string
553 function h11 = getBorrowing(obj, boundary)
554 assertIsMember(boundary, {'w', 'e', 's', 'n'})
555
556 switch boundary
557 case {'w','e'}
558 h11 = obj.h11{1};
559 case {'s', 'n'}
560 h11 = obj.h11{2};
561 end
562 end
563
564 % Returns the outward unit normal vector for the boundary specified by the string boundary.
565 function nu = getNormal(obj, boundary)
566 assertIsMember(boundary, {'w', 'e', 's', 'n'})
567
568 switch boundary
569 case 'w'
570 nu = [-1,0];
571 case 'e'
572 nu = [1,0];
573 case 's'
574 nu = [0,-1];
575 case 'n'
576 nu = [0,1];
577 end
578 end
579
580 % Returns the boundary operator op for the boundary specified by the string boundary.
581 % op -- string
582 function o = getBoundaryOperator(obj, op, boundary)
583 assertIsMember(boundary, {'w', 'e', 's', 'n'})
584 assertIsMember(op, {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2'})
585
586 switch op
587 case {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2'}
588 o = obj.([op, '_', boundary]);
589 end
590
591 end
592
593 % Returns the boundary operator op for the boundary specified by the string boundary.
594 % op -- string
595 function o = getBoundaryOperatorForScalarField(obj, op, boundary)
596 assertIsMember(boundary, {'w', 'e', 's', 'n'})
597 assertIsMember(op, {'e'})
598
599 switch op
600
601 case 'e'
602 o = obj.(['e_scalar', '_', boundary]);
603 end
604
605 end
606
607 % Returns the boundary operator T_ij (cell format) for the boundary specified by the string boundary.
608 % Formula: tau_i = T_ij u_j
609 % op -- string
610 function T = getBoundaryTractionOperator(obj, boundary)
611 assertIsMember(boundary, {'w', 'e', 's', 'n'})
612
613 T = obj.(['T', '_', boundary]);
614 end
615
616 % Returns square boundary quadrature matrix, of dimension
617 % corresponding to the number of boundary unknowns
618 %
619 % boundary -- string
620 function H = getBoundaryQuadrature(obj, boundary)
621 assertIsMember(boundary, {'w', 'e', 's', 'n'})
622
623 H = obj.getBoundaryQuadratureForScalarField(boundary);
624 I_dim = speye(obj.dim, obj.dim);
625 H = kron(H, I_dim);
626 end
627
628 % Returns square boundary quadrature matrix, of dimension
629 % corresponding to the number of boundary grid points
630 %
631 % boundary -- string
632 function H_b = getBoundaryQuadratureForScalarField(obj, boundary)
633 assertIsMember(boundary, {'w', 'e', 's', 'n'})
634
635 H_b = obj.(['H_', boundary]);
636 end
637
638 function N = size(obj)
639 N = obj.dim*prod(obj.m);
640 end
641 end
642 end