Mercurial > repos > public > sbplib

comparison +scheme/ViscoElastic2d.m @ 1307:fcca6ad8b102 feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add diffOp for viscoElastic
author Martin Almquist <malmquist@stanford.edu>
date Sun, 19 Jul 2020 20:30:16 -0700
parents
children 5016f3f3badb
comparison
equal deleted inserted replaced
1306:633757e582e5 1307:fcca6ad8b102
1 classdef ViscoElastic2d < scheme.Scheme
2
3 % Discretizes the visco-elastic wave equation in curvilinear coordinates.
4 % Assumes fully compatible operators.
5
6 properties
7 m % Number of points in each direction, possibly a vector
8 h % Grid spacing
9
10 grid
11 dim
12
13 order % Order of accuracy for the approximation
14
15 % Diagonal matrices for variable coefficients
16 % J, Ji
17 RHO % Density
18 C % Elastic stiffness tensor
19 ETA % Effective viscosity, used in strain rate eq
20
21 D % Total operator
22 Delastic % Elastic operator (momentum balance)
23 Dviscous % Acts on viscous strains in momentum balance
24 DstrainRate % Acts on u and gamma, returns strain rate gamma_t
25
26 D1, D1Tilde % Physical derivatives
27 sigma % Cell matrix of physical stress operators
28
29 % Inner products
30 H
31
32 % Boundary inner products (for scalar field)
33 H_w, H_e, H_s, H_n
34
35 % Restriction operators
36 Eu, Egamma % Pick out all components of u/gamma
37 eU, eGamma % Pick out one specific component
38
39 % Bundary restriction ops
40 e_scalar_w, e_scalar_e, e_scalar_s, e_scalar_n
41
42 n_w, n_e, n_s, n_n % Physical normals
43
44 elasticObj
45
46 end
47
48 methods
49
50 % The coefficients can either be function handles or grid functions
51 function obj = ViscoElastic2d(g, order, rho, C, eta)
52 default_arg('rho', @(x,y) 0*x+1);
53 default_arg('eta', @(x,y) 0*x+1);
54 dim = 2;
55
56 C_default = cell(dim,dim,dim,dim);
57 for i = 1:dim
58 for j = 1:dim
59 for k = 1:dim
60 for l = 1:dim
61 C_default{i,j,k,l} = @(x,y) 0*x ;
62 end
63 end
64 end
65 end
66 default_arg('C', C_default);
67
68 assert(isa(g, 'grid.Curvilinear'));
69
70 if isa(rho, 'function_handle')
71 rho = grid.evalOn(g, rho);
72 end
73
74 if isa(eta, 'function_handle')
75 eta = grid.evalOn(g, eta);
76 end
77
78 C_mat = cell(dim,dim,dim,dim);
79 for i = 1:dim
80 for j = 1:dim
81 for k = 1:dim
82 for l = 1:dim
83 if isa(C{i,j,k,l}, 'function_handle')
84 C{i,j,k,l} = grid.evalOn(g, C{i,j,k,l});
85 end
86 C_mat{i,j,k,l} = spdiag(C{i,j,k,l});
87 end
88 end
89 end
90 end
91 obj.C = C_mat;
92
93 elasticObj = scheme.Elastic2dCurvilinearAnisotropic(g, order, rho, C);
94
95 % Construct a pair of first derivatives
96 K = elasticObj.K;
97 for i = 1:dim
98 for j = 1:dim
99 K{i,j} = spdiag(K{i,j});
100 end
101 end
102 J = elasticObj.J;
103 Ji = elasticObj.Ji;
104 D_ref = elasticObj.refObj.D1;
105 D1 = cell(dim, 1);
106 D1Tilde = cell(dim, 1);
107 for i = 1:dim
108 D1{i} = 0*D_ref{i};
109 D1Tilde{i} = 0*D_ref{i};
110 for j = 1:dim
111 D1{i} = D1{i} + K{i,j}*D_ref{j};
112 D1Tilde{i} = D1Tilde{i} + Ji*D_ref{j}*J*K{i,j};
113 end
114 end
115 obj.D1 = D1;
116 obj.D1Tilde = D1Tilde;
117
118 eU = elasticObj.E;
119
120 % Storage order for gamma: 11-12-21-22.
121 I = speye(g.N(), g.N());
122 eGamma = cell(dim, dim);
123 e = cell(dim, dim);
124 e{1,1} = [1;0;0;0];
125 e{1,2} = [0;1;0;0];
126 e{2,1} = [0;0;1;0];
127 e{2,2} = [0;0;0;1];
128 for i = 1:dim
129 for j = 1:dim
130 eGamma{i,j} = kron(I, e{i,j});
131 end
132 end
133
134 % Store u first, then gamma
135 mU = dim*g.N();
136 mGamma = dim^2*g.N();
137 Iu = speye(mU, mU);
138 Igamma = speye(mGamma, mGamma);
139
140 Eu = cell2mat({Iu, sparse(mU, mGamma)})';
141 Egamma = cell2mat({sparse(mGamma, mU), Igamma})';
142
143 for i = 1:dim
144 eU{i} = Eu*eU{i};
145 end
146 for i = 1:dim
147 for j = 1:dim
148 eGamma{i,j} = Egamma*eGamma{i,j};
149 end
150 end
151
152 obj.eGamma = eGamma;
153 obj.eU = eU;
154 obj.Egamma = Egamma;
155 obj.Eu = Eu;
156
157 % Build stress operator
158 sigma = cell(dim, dim);
159 C = obj.C;
160 for i = 1:dim
161 for j = 1:dim
162 sigma{i,j} = spalloc(g.N(), (dim^2 + dim)*g.N(), order^2*g.N());
163 for k = 1:dim
164 for l = 1:dim
165 sigma{i,j} = sigma{i,j} + C{i,j,k,l}*(D1{k}*eU{l}' - eGamma{k,l}');
166 end
167 end
168 end
169 end
170
171 % Elastic operator
172 Delastic = Eu*elasticObj.D*Eu';
173
174 % Add viscous strains to momentum balance
175 RHOi = spdiag(1./rho);
176 Dviscous = spalloc((dim^2 + dim)*g.N(), (dim^2 + dim)*g.N(), order^2*(dim^2 + dim)*g.N());
177 for i = 1:dim
178 for j = 1:dim
179 for k = 1:dim
180 for l = 1:dim
181 Dviscous = Dviscous - eU{j}*RHOi*D1Tilde{i}*C{i,j,k,l}*eGamma{k,l}';
182 end
183 end
184 end
185 end
186
187 ETA = spdiag(eta);
188 DstrainRate = 0*Delastic;
189 for i = 1:dim
190 for j = 1:dim
191 DstrainRate = DstrainRate + eGamma{i,j}*(ETA\sigma{i,j});
192 end
193 end
194
195 obj.D = Delastic + Dviscous + DstrainRate;
196 obj.Delastic = Delastic;
197 obj.Dviscous = Dviscous;
198 obj.DstrainRate = DstrainRate;
199 obj.sigma = sigma;
200
201 %---- Set remaining object properties ------
202 obj.RHO = elasticObj.RHO;
203 obj.ETA = ETA;
204 obj.H = elasticObj.H;
205
206 obj.n_w = elasticObj.n_w;
207 obj.n_e = elasticObj.n_e;
208 obj.n_s = elasticObj.n_s;
209 obj.n_n = elasticObj.n_n;
210
211 obj.H_w = elasticObj.H_w;
212 obj.H_e = elasticObj.H_e;
213 obj.H_s = elasticObj.H_s;
214 obj.H_n = elasticObj.H_n;
215
216 obj.e_scalar_w = elasticObj.e_scalar_w;
217 obj.e_scalar_e = elasticObj.e_scalar_e;
218 obj.e_scalar_s = elasticObj.e_scalar_s;
219 obj.e_scalar_n = elasticObj.e_scalar_n;
220
221 % Misc.
222 obj.elasticObj = elasticObj;
223 obj.m = elasticObj.m;
224 obj.h = elasticObj.h;
225
226 obj.order = order;
227 obj.grid = g;
228 obj.dim = dim;
229
230 end
231
232
233 % Closure functions return the operators applied to the own domain to close the boundary
234 % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other doamin.
235 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
236 % bc is a cell array of component and bc type, e.g. {1, 'd'} for Dirichlet condition
237 % on the first component. Can also be e.g.
238 % {'normal', 'd'} or {'tangential', 't'} for conditions on
239 % tangential/normal component.
240 % data is a function returning the data that should be applied at the boundary.
241 % neighbour_scheme is an instance of Scheme that should be interfaced to.
242 % neighbour_boundary is a string specifying which boundary to interface to.
243
244 % For displacement bc:
245 % bc = {comp, 'd', dComps},
246 % where
247 % dComps = vector of components with displacement BC. Default: 1:dim.
248 % In this way, we can specify one BC at a time even though the SATs depend on all BC.
249 function [closure, penalty] = boundary_condition(obj, boundary, bc, tuning)
250 default_arg('tuning', 1.0);
251 assert( iscell(bc), 'The BC type must be a 2x1 or 3x1 cell array' );
252
253 component = bc{1};
254 type = bc{2};
255 dim = obj.dim;
256
257 n = obj.getNormal(boundary);
258 H_gamma = obj.getBoundaryQuadratureForScalarField(boundary);
259 e = obj.getBoundaryOperatorForScalarField('e', boundary);
260
261 H = obj.H;
262 RHO = obj.RHO;
263 C = obj.C;
264 Eu = obj.Eu;
265 eU = obj.eU;
266 eGamma = obj.eGamma;
267
268 switch type
269 case {'F','f','Free','free','traction','Traction','t','T'}
270 [closure, penalty] = obj.elasticObj.boundary_condition(boundary, bc, tuning);
271 closure = Eu*closure*Eu';
272 penalty = Eu*penalty;
273
274 j = component;
275 for i = 1:dim
276 for k = 1:dim
277 for l = 1:dim
278 closure = closure + eU{j}*( (RHO*H)\(C{i,j,k,l}*e*H_gamma*n{i}*e'*eGamma{k,l}') );
279 end
280 end
281 end
282 end
283
284 end
285
286 function [closure, penalty] = displacementBCNormalTangential(obj, boundary, bc, tuning)
287 disp('WARNING: displacementBCNormalTangential is only guaranteed to work for displacement BC on one component and traction on the other');
288 u = obj;
289
290 component = bc{1};
291 type = bc{2};
292
293 switch component
294 case 'n'
295 en = u.getBoundaryOperator('en', boundary);
296 tau_n = u.getBoundaryOperator('tau_n', boundary);
297 N = u.getNormal(boundary);
298 case 't'
299 en = u.getBoundaryOperator('et', boundary);
300 tau_n = u.getBoundaryOperator('tau_t', boundary);
301 N = u.getTangent(boundary);
302 end
303
304 % Operators
305 e = u.getBoundaryOperatorForScalarField('e', boundary);
306 h11 = u.getBorrowing(boundary);
307 n = u.getNormal(boundary);
308
309 C = u.C;
310 Ji = u.Ji;
311 s = spdiag(u.(['s_' boundary]));
312 m_tot = u.grid.N();
313
314 Hi = u.E{1}*inv(u.H)*u.E{1}' + u.E{2}*inv(u.H)*u.E{2}';
315 RHOi = u.E{1}*inv(u.RHO)*u.E{1}' + u.E{2}*inv(u.RHO)*u.E{2}';
316
317 H_gamma = u.getBoundaryQuadratureForScalarField(boundary);
318 dim = u.dim;
319
320 % Preallocate
321 [~, m_int] = size(H_gamma);
322 closure = sparse(dim*m_tot, dim*m_tot);
323 penalty = sparse(dim*m_tot, m_int);
324
325 % Term 1: The symmetric term
326 Z = sparse(m_int, m_int);
327 for i = 1:dim
328 for j = 1:dim
329 for l = 1:dim
330 for k = 1:dim
331 Z = Z + n{i}*N{j}*e'*Ji*C{i,j,k,l}*e*n{k}*N{l};
332 end
333 end
334 end
335 end
336
337 Z = -tuning*dim*1/h11*s*Z;
338 closure = closure + en*H_gamma*Z*en';
339 penalty = penalty - en*H_gamma*Z;
340
341 % Term 2: The symmetrizing term
342 closure = closure + tau_n*H_gamma*en';
343 penalty = penalty - tau_n*H_gamma;
344
345 % Multiply all normal component terms by inverse of density x quadraure
346 closure = RHOi*Hi*closure;
347 penalty = RHOi*Hi*penalty;
348 end
349
350 % type Struct that specifies the interface coupling.
351 % Fields:
352 % -- tuning: penalty strength, defaults to 1.0
353 % -- interpolation: type of interpolation, default 'none'
354 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type)
355
356 defaultType.tuning = 1.0;
357 defaultType.interpolation = 'none';
358 defaultType.type = 'standard';
359 default_struct('type', defaultType);
360
361 switch type.type
362 case 'standard'
363 [closure, penalty] = obj.refObj.interface(boundary,neighbour_scheme.refObj,neighbour_boundary,type);
364 case 'frictionalFault'
365 [closure, penalty] = obj.interfaceFrictionalFault(boundary,neighbour_scheme,neighbour_boundary,type);
366 end
367
368 end
369
370 function [closure, penalty] = interfaceFrictionalFault(obj,boundary,neighbour_scheme,neighbour_boundary,type)
371 tuning = type.tuning;
372
373 % u denotes the solution in the own domain
374 % v denotes the solution in the neighbour domain
375
376 u = obj;
377 v = neighbour_scheme;
378
379 % Operators, u side
380 e_u = u.getBoundaryOperatorForScalarField('e', boundary);
381 en_u = u.getBoundaryOperator('en', boundary);
382 tau_n_u = u.getBoundaryOperator('tau_n', boundary);
383 h11_u = u.getBorrowing(boundary);
384 n_u = u.getNormal(boundary);
385
386 C_u = u.C;
387 Ji_u = u.Ji;
388 s_u = spdiag(u.(['s_' boundary]));
389 m_tot_u = u.grid.N();
390
391 % Operators, v side
392 e_v = v.getBoundaryOperatorForScalarField('e', neighbour_boundary);
393 en_v = v.getBoundaryOperator('en', neighbour_boundary);
394 tau_n_v = v.getBoundaryOperator('tau_n', neighbour_boundary);
395 h11_v = v.getBorrowing(neighbour_boundary);
396 n_v = v.getNormal(neighbour_boundary);
397
398 C_v = v.C;
399 Ji_v = v.Ji;
400 s_v = spdiag(v.(['s_' neighbour_boundary]));
401 m_tot_v = v.grid.N();
402
403 % Operators that are only required for own domain
404 Hi = u.E{1}*inv(u.H)*u.E{1}' + u.E{2}*inv(u.H)*u.E{2}';
405 RHOi = u.E{1}*inv(u.RHO)*u.E{1}' + u.E{2}*inv(u.RHO)*u.E{2}';
406
407 % Shared operators
408 H_gamma = u.getBoundaryQuadratureForScalarField(boundary);
409 dim = u.dim;
410
411 % Preallocate
412 [~, m_int] = size(H_gamma);
413 closure = sparse(dim*m_tot_u, dim*m_tot_u);
414 penalty = sparse(dim*m_tot_u, dim*m_tot_v);
415
416 % Continuity of normal displacement, term 1: The symmetric term
417 Z_u = sparse(m_int, m_int);
418 Z_v = sparse(m_int, m_int);
419 for i = 1:dim
420 for j = 1:dim
421 for l = 1:dim
422 for k = 1:dim
423 Z_u = Z_u + n_u{i}*n_u{j}*e_u'*Ji_u*C_u{i,j,k,l}*e_u*n_u{k}*n_u{l};
424 Z_v = Z_v + n_v{i}*n_v{j}*e_v'*Ji_v*C_v{i,j,k,l}*e_v*n_v{k}*n_v{l};
425 end
426 end
427 end
428 end
429
430 Z = -tuning*dim*( 1/(4*h11_u)*s_u*Z_u + 1/(4*h11_v)*s_v*Z_v );
431 closure = closure + en_u*H_gamma*Z*en_u';
432 penalty = penalty + en_u*H_gamma*Z*en_v';
433
434 % Continuity of normal displacement, term 2: The symmetrizing term
435 closure = closure + 1/2*tau_n_u*H_gamma*en_u';
436 penalty = penalty + 1/2*tau_n_u*H_gamma*en_v';
437
438 % Continuity of normal traction
439 closure = closure - 1/2*en_u*H_gamma*tau_n_u';
440 penalty = penalty + 1/2*en_u*H_gamma*tau_n_v';
441
442 % Multiply all normal component terms by inverse of density x quadraure
443 closure = RHOi*Hi*closure;
444 penalty = RHOi*Hi*penalty;
445
446 % ---- Tangential tractions are imposed just like traction BC ------
447 closure = closure + obj.boundary_condition(boundary, {'t', 't'});
448
449 end
450
451
452 % Returns h11 for the boundary specified by the string boundary.
453 % op -- string
454 function h11 = getBorrowing(obj, boundary)
455 assertIsMember(boundary, {'w', 'e', 's', 'n'})
456
457 switch boundary
458 case {'w','e'}
459 h11 = obj.refObj.h11{1};
460 case {'s', 'n'}
461 h11 = obj.refObj.h11{2};
462 end
463 end
464
465 % Returns the outward unit normal vector for the boundary specified by the string boundary.
466 % n is a cell of diagonal matrices for each normal component, n{1} = n_1, n{2} = n_2.
467 function n = getNormal(obj, boundary)
468 assertIsMember(boundary, {'w', 'e', 's', 'n'})
469
470 n = obj.(['n_' boundary]);
471 end
472
473 % Returns the unit tangent vector for the boundary specified by the string boundary.
474 % t is a cell of diagonal matrices for each normal component, t{1} = t_1, t{2} = t_2.
475 function t = getTangent(obj, boundary)
476 assertIsMember(boundary, {'w', 'e', 's', 'n'})
477
478 t = obj.(['tangent_' boundary]);
479 end
480
481 % Returns the boundary operator op for the boundary specified by the string boundary.
482 % op -- string
483 function o = getBoundaryOperator(obj, op, boundary)
484 assertIsMember(boundary, {'w', 'e', 's', 'n'})
485 assertIsMember(op, {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2', 'en', 'et', 'tau_n', 'tau_t'})
486
487 o = obj.([op, '_', boundary]);
488
489 end
490
491 % Returns the boundary operator op for the boundary specified by the string boundary.
492 % op -- string
493 function o = getBoundaryOperatorForScalarField(obj, op, boundary)
494 assertIsMember(boundary, {'w', 'e', 's', 'n'})
495 assertIsMember(op, {'e'})
496
497 switch op
498
499 case 'e'
500 o = obj.(['e_scalar', '_', boundary]);
501 end
502
503 end
504
505 % Returns the boundary operator T_ij (cell format) for the boundary specified by the string boundary.
506 % Formula: tau_i = T_ij u_j
507 % op -- string
508 function T = getBoundaryTractionOperator(obj, boundary)
509 assertIsMember(boundary, {'w', 'e', 's', 'n'})
510
511 T = obj.(['T', '_', boundary]);
512 end
513
514 % Returns square boundary quadrature matrix, of dimension
515 % corresponding to the number of boundary unknowns
516 %
517 % boundary -- string
518 function H = getBoundaryQuadrature(obj, boundary)
519 assertIsMember(boundary, {'w', 'e', 's', 'n'})
520
521 H = obj.getBoundaryQuadratureForScalarField(boundary);
522 I_dim = speye(obj.dim, obj.dim);
523 H = kron(H, I_dim);
524 end
525
526 % Returns square boundary quadrature matrix, of dimension
527 % corresponding to the number of boundary grid points
528 %
529 % boundary -- string
530 function H_b = getBoundaryQuadratureForScalarField(obj, boundary)
531 assertIsMember(boundary, {'w', 'e', 's', 'n'})
532
533 H_b = obj.(['H_', boundary]);
534 end
535
536 function N = size(obj)
537 N = (obj.dim + obj.dim^2)*prod(obj.m);
538 end
539 end
540 end