Mercurial > repos > public > sbplib

comparison +scheme/elasticDilationVariable.m @ 674:dd84b8862aa8 feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
First implementation of elastic dilation variable. Constant coeff is stable.
author Martin Almquist <malmquist@stanford.edu>
date Tue, 16 Jan 2018 17:41:55 -0800
parents
children 90bf651abc7c
comparison
equal deleted inserted replaced
673:9e1d2351f539 674:dd84b8862aa8
1 classdef elasticDilationVariable < scheme.Scheme
2 properties
3 m % Number of points in each direction, possibly a vector
4 h % Grid spacing
5
6 grid
7 dim
8
9 order % Order accuracy for the approximation
10
11 A % Variable coefficient lambda of the operator (as diagonal matrix here)
12 RHO % Density (as diagonal matrix here)
13
14 D % Total operator
15 D1 % First derivatives
16 D2 % Second derivatives
17 Div % Divergence operator used for BC
18 H, Hi % Inner products
19 phi % Borrowing constant for (d1 - e^T*D1) from R
20 H11 % First element of H
21 e_l, e_r
22 d1_l, d1_r % Normal derivatives at the boundary
23 E % E{i}^T picks out component i
24
25 H_boundary % Boundary inner products
26
27 A_boundary_l % Variable coefficient at boundaries
28 A_boundary_r %
29 end
30
31 methods
32 % Implements the shear part of the elastic wave equation, i.e.
33 % rho u_{i,tt} = d_i a d_j u_j + d_j a d_j u_i
34 % where a = mu.
35
36 function obj = elasticDilationVariable(g ,order, a_fun, rho_fun, opSet)
37 default_arg('opSet',@sbp.D2Variable);
38 default_arg('a_fun', @(x,y) 0*x+1);
39 default_arg('rho_fun', @(x,y) 0*x+1);
40 dim = 2;
41
42 assert(isa(g, 'grid.Cartesian'))
43
44 a = grid.evalOn(g, a_fun);
45 rho = grid.evalOn(g, rho_fun);
46 m = g.size();
47 m_tot = g.N();
48
49 h = g.scaling();
50 L = (m-1).*h;
51
52 % 1D operators
53 ops = cell(dim,1);
54 for i = 1:dim
55 ops{i} = opSet(m(i), {0, L(i)}, order);
56 end
57
58 % Borrowing constants
59 beta = ops{1}.borrowing.R.delta_D;
60 obj.H11 = ops{1}.borrowing.H11;
61 obj.phi = beta/obj.H11;
62
63 I = cell(dim,1);
64 D1 = cell(dim,1);
65 D2 = cell(dim,1);
66 H = cell(dim,1);
67 Hi = cell(dim,1);
68 e_l = cell(dim,1);
69 e_r = cell(dim,1);
70 d1_l = cell(dim,1);
71 d1_r = cell(dim,1);
72
73 for i = 1:dim
74 I{i} = speye(m(i));
75 D1{i} = ops{i}.D1;
76 D2{i} = ops{i}.D2;
77 H{i} = ops{i}.H;
78 Hi{i} = ops{i}.HI;
79 e_l{i} = ops{i}.e_l;
80 e_r{i} = ops{i}.e_r;
81 d1_l{i} = ops{i}.d1_l;
82 d1_r{i} = ops{i}.d1_r;
83 end
84
85 %====== Assemble full operators ========
86 A = spdiag(a);
87 obj.A = A;
88 RHO = spdiag(rho);
89 obj.RHO = RHO;
90
91
92 obj.D1 = cell(dim,1);
93 obj.D2 = cell(dim,1);
94 obj.e_l = cell(dim,1);
95 obj.e_r = cell(dim,1);
96 obj.d1_l = cell(dim,1);
97 obj.d1_r = cell(dim,1);
98
99 % D1
100 obj.D1{1} = kron(D1{1},I{2});
101 obj.D1{2} = kron(I{1},D1{2});
102
103 % Boundary operators
104 obj.e_l{1} = kron(e_l{1},I{2});
105 obj.e_l{2} = kron(I{1},e_l{2});
106 obj.e_r{1} = kron(e_r{1},I{2});
107 obj.e_r{2} = kron(I{1},e_r{2});
108
109 obj.d1_l{1} = kron(d1_l{1},I{2});
110 obj.d1_l{2} = kron(I{1},d1_l{2});
111 obj.d1_r{1} = kron(d1_r{1},I{2});
112 obj.d1_r{2} = kron(I{1},d1_r{2});
113
114 % D2
115 for i = 1:dim
116 obj.D2{i} = sparse(m_tot);
117 end
118 ind = grid.funcToMatrix(g, 1:m_tot);
119
120 for i = 1:m(2)
121 D = D2{1}(a(ind(:,i)));
122 p = ind(:,i);
123 obj.D2{1}(p,p) = D;
124 end
125
126 for i = 1:m(1)
127 D = D2{2}(a(ind(i,:)));
128 p = ind(i,:);
129 obj.D2{2}(p,p) = D;
130 end
131
132 % Quadratures
133 obj.H = kron(H{1},H{2});
134 obj.Hi = inv(obj.H);
135 obj.H_boundary = cell(dim,1);
136 obj.H_boundary{1} = H{2};
137 obj.H_boundary{2} = H{1};
138
139 % Boundary coefficient matrices and quadratures
140 obj.A_boundary_l = cell(dim,1);
141 obj.A_boundary_r = cell(dim,1);
142 for i = 1:dim
143 obj.A_boundary_l{i} = obj.e_l{i}'*A*obj.e_l{i};
144 obj.A_boundary_r{i} = obj.e_r{i}'*A*obj.e_r{i};
145 end
146
147 % E{i}^T picks out component i.
148 E = cell(dim,1);
149 I = speye(m_tot,m_tot);
150 for i = 1:dim
151 e = sparse(dim,1);
152 e(i) = 1;
153 E{i} = kron(I,e);
154 end
155 obj.E = E;
156
157 % Differentiation matrix D (without SAT)
158 D2 = obj.D2;
159 D1 = obj.D1;
160 D = sparse(dim*m_tot,dim*m_tot);
161 d = @kroneckerDelta; % Kronecker delta
162 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta
163 for i = 1:dim
164 for j = 1:dim
165 D = D + E{i}*inv(RHO)*( d(i,j)*D2{i}*E{j}' +...
166 db(i,j)*D1{i}*A*D1{j}*E{j}' ...
167 );
168 end
169 end
170 obj.D = D;
171 %=========================================%
172
173 % Divergence operator for BC
174 Div = cell(dim,1);
175 for i = 1:dim
176 Div{i} = sparse(m_tot,dim*m_tot);
177 for j = 1:dim
178 Div{i} = Div{i} + d(i,j)*(obj.e_l{i}*obj.d1_l{i}' + obj.e_r{i}*obj.d1_r{i}')*E{j}' ...
179 + db(i,j)*obj.D1{j}*E{j}';
180 end
181 end
182 obj.Div = Div;
183
184 % Misc.
185 obj.m = m;
186 obj.h = h;
187 obj.order = order;
188 obj.grid = g;
189 obj.dim = dim;
190
191 end
192
193
194 % Closure functions return the operators applied to the own domain to close the boundary
195 % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other doamin.
196 % Here penalty{i,j} enforces data component j on solution component i
197 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
198 % type is a string specifying the type of boundary condition if there are several.
199 % data is a function returning the data that should be applied at the boundary.
200 % neighbour_scheme is an instance of Scheme that should be interfaced to.
201 % neighbour_boundary is a string specifying which boundary to interface to.
202 function [closure, penalty] = boundary_condition(obj, boundary, type, parameter)
203 default_arg('type','free');
204 default_arg('parameter', []);
205
206 delta = @kroneckerDelta; % Kronecker delta
207 delta_b = @(i,j) 1-delta(i,j); % Logical not of Kronecker delta
208
209 % j is the coordinate direction of the boundary
210 % nj: outward unit normal component.
211 % nj = -1 for west, south, bottom boundaries
212 % nj = 1 for east, north, top boundaries
213 [j, nj] = obj.get_boundary_number(boundary);
214 switch nj
215 case 1
216 e = obj.e_r;
217 d = obj.d1_r;
218 case -1
219 e = obj.e_l;
220 d = obj.d1_l;
221 end
222
223 E = obj.E;
224 Hi = obj.Hi;
225 H_gamma = obj.H_boundary{j};
226 A = obj.A;
227 RHO = obj.RHO;
228 D1 = obj.D1;
229
230 phi = obj.phi;
231 H11 = obj.H11;
232 h = obj.h;
233
234 switch type
235 % Dirichlet boundary condition
236 case {'D','d','dirichlet'}
237 error('Not implemented');
238 tuning = 1.2;
239 phi = obj.phi;
240
241 closures = cell(obj.dim,1);
242 penalties = cell(obj.dim,obj.dim);
243 % Loop over components
244 for i = 1:obj.dim
245 H_gamma_i = obj.H_boundary{i};
246 sigma_ij = tuning*delta(i,j)*2/(gamma*h(j)) +...
247 tuning*delta_b(i,j)*(2/(H11*h(j)) + 2/(H11*h(j)*phi));
248
249 ci = E{i}*inv(RHO)*nj*Hi*...
250 ( (e{j}*H_gamma*e{j}'*A*e{j}*d{j}')'*E{i}' + ...
251 delta(i,j)*(e{j}*H_gamma*e{j}'*A*e{j}*d{j}')'*E{j}' ...
252 ) ...
253 - sigma_ij*E{i}*inv(RHO)*Hi*A*e{j}*H_gamma*e{j}'*E{i}';
254
255 cj = E{j}*inv(RHO)*nj*Hi*...
256 ( delta_b(i,j)*(e{j}*H_gamma*e{j}'*A*D1{i})'*E{i}' ...
257 );
258
259 if isempty(closures{i})
260 closures{i} = ci;
261 else
262 closures{i} = closures{i} + ci;
263 end
264
265 if isempty(closures{j})
266 closures{j} = cj;
267 else
268 closures{j} = closures{j} + cj;
269 end
270
271 pii = - E{i}*inv(RHO)*nj*Hi*...
272 ( (H_gamma*e{j}'*A*e{j}*d{j}')' + ...
273 delta(i,j)*(H_gamma*e{j}'*A*e{j}*d{j}')' ...
274 ) ...
275 + sigma_ij*E{i}*inv(RHO)*Hi*A*e{j}*H_gamma;
276
277 pji = - E{j}*inv(RHO)*nj*Hi*...
278 ( delta_b(i,j)*(H_gamma*e{j}'*A*D1{i})' );
279
280 % Dummies
281 pij = - 0*E{i}*e{j};
282 pjj = - 0*E{j}*e{j};
283
284 if isempty(penalties{i,i})
285 penalties{i,i} = pii;
286 else
287 penalties{i,i} = penalties{i,i} + pii;
288 end
289
290 if isempty(penalties{j,i})
291 penalties{j,i} = pji;
292 else
293 penalties{j,i} = penalties{j,i} + pji;
294 end
295
296 if isempty(penalties{i,j})
297 penalties{i,j} = pij;
298 else
299 penalties{i,j} = penalties{i,j} + pij;
300 end
301
302 if isempty(penalties{j,j})
303 penalties{j,j} = pjj;
304 else
305 penalties{j,j} = penalties{j,j} + pjj;
306 end
307 end
308 [rows, cols] = size(closures{1});
309 closure = sparse(rows, cols);
310 for i = 1:obj.dim
311 closure = closure + closures{i};
312 end
313 penalty = penalties;
314
315 % Free boundary condition
316 case {'F','f','Free','free'}
317 closures = cell(obj.dim,1);
318 penalties = cell(obj.dim,obj.dim);
319
320 % Divergence operator
321 Div = obj.Div{j};
322
323 % Loop over components
324 %for i = 1:obj.dim
325 closure = -nj*E{j}*inv(RHO)*Hi*e{j} ...
326 *H_gamma*e{j}'*A*e{j}*e{j}'*Div;
327 penalty = nj*E{j}*inv(RHO)*Hi*e{j} ...
328 *H_gamma*e{j}'*A*e{j};
329 %end
330 % [rows, cols] = size(closures{1});
331 % closure = sparse(rows, cols);
332 % for i = 1:obj.dim
333 % closure = closure + closures{i};
334 % for j = 1:obj.dim
335 % if i~=j
336 % [rows cols] = size(penalties{j,j});
337 % penalties{i,j} = sparse(rows,cols);
338 % end
339 % end
340 % end
341 % penalty = penalties;
342
343
344 % Unknown boundary condition
345 otherwise
346 error('No such boundary condition: type = %s',type);
347 end
348 end
349
350 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
351 % u denotes the solution in the own domain
352 % v denotes the solution in the neighbour domain
353 tuning = 1.2;
354 % tuning = 20.2;
355 error('Interface not implemented');
356 end
357
358 % Returns the coordinate number and outward normal component for the boundary specified by the string boundary.
359 function [j, nj] = get_boundary_number(obj, boundary)
360
361 switch boundary
362 case {'w','W','west','West', 'e', 'E', 'east', 'East'}
363 j = 1;
364 case {'s','S','south','South', 'n', 'N', 'north', 'North'}
365 j = 2;
366 otherwise
367 error('No such boundary: boundary = %s',boundary);
368 end
369
370 switch boundary
371 case {'w','W','west','West','s','S','south','South'}
372 nj = -1;
373 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
374 nj = 1;
375 end
376 end
377
378 % Returns the coordinate number and outward normal component for the boundary specified by the string boundary.
379 function [return_op] = get_boundary_operator(obj, op, boundary)
380
381 switch boundary
382 case {'w','W','west','West', 'e', 'E', 'east', 'East'}
383 j = 1;
384 case {'s','S','south','South', 'n', 'N', 'north', 'North'}
385 j = 2;
386 otherwise
387 error('No such boundary: boundary = %s',boundary);
388 end
389
390 switch op
391 case 'e'
392 switch boundary
393 case {'w','W','west','West','s','S','south','South'}
394 return_op = obj.e_l{j};
395 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
396 return_op = obj.e_r{j};
397 end
398 case 'd'
399 switch boundary
400 case {'w','W','west','West','s','S','south','South'}
401 return_op = obj.d_l{j};
402 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
403 return_op = obj.d_r{j};
404 end
405 otherwise
406 error(['No such operator: operatr = ' op]);
407 end
408
409 end
410
411 function N = size(obj)
412 N = prod(obj.m);
413 end
414 end
415 end