Mercurial > repos > public > sbplib

comparison +scheme/elasticShearVariable.m @ 664:8e6dfd22fc59 feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Free BC now yields symmetric negative semidefinite discretization.
author Martin Almquist <malmquist@stanford.edu>
date Fri, 15 Dec 2017 16:23:10 -0800
parents b45ec2b28cc2
children ed853945ee99
comparison
equal deleted inserted replaced
663:b45ec2b28cc2 664:8e6dfd22fc59
2 properties 2 properties
3 m % Number of points in each direction, possibly a vector 3 m % Number of points in each direction, possibly a vector
4 h % Grid spacing 4 h % Grid spacing
5 5
6 grid 6 grid
7 dim
7 8
8 order % Order accuracy for the approximation 9 order % Order accuracy for the approximation
9 10
10 a % Variable coefficient lambda of the operator 11 A % Variable coefficient lambda of the operator (as diagonal matrix here)
11 rho % Density 12 RHO % Density (as diagonal matrix here)
12 13
13 D % Total operator 14 D % Total operator
14 D1 % First derivatives 15 D1 % First derivatives
15 D2 % Second derivatives 16 D2 % Second derivatives
16 H, Hi % Inner products 17 H, Hi % Inner products
17 e_l, e_r 18 e_l, e_r
18 d_l, d_r % Normal derivatives at the boundary 19 d1_l, d1_r % Normal derivatives at the boundary
20 E % E{i}^T picks out component i
21
19 H_boundary % Boundary inner products 22 H_boundary % Boundary inner products
23
24 A_boundary_l % Variable coefficient at boundaries
25 A_boundary_r %
20 end 26 end
21 27
22 methods 28 methods
23 % Implements the shear part of the elastic wave equation, i.e. 29 % Implements the shear part of the elastic wave equation, i.e.
24 % rho u_{i,tt} = d_i a d_j u_j + d_j a d_j u_i 30 % rho u_{i,tt} = d_i a d_j u_j + d_j a d_j u_i
54 e_r = cell(dim,1); 60 e_r = cell(dim,1);
55 d1_l = cell(dim,1); 61 d1_l = cell(dim,1);
56 d1_r = cell(dim,1); 62 d1_r = cell(dim,1);
57 63
58 for i = 1:dim 64 for i = 1:dim
59 I{i} = speye{m(i)); 65 I{i} = speye(m(i));
60 D1{i} = ops{i}.D1; 66 D1{i} = ops{i}.D1;
61 D2{i} = ops{i}.D2; 67 D2{i} = ops{i}.D2;
62 H{i} = ops{i}.H; 68 H{i} = ops{i}.H;
63 Hi{i} = ops{i}.HI; 69 Hi{i} = ops{i}.HI;
64 e_l{i} = ops{i}.e_l; 70 e_l{i} = ops{i}.e_l;
67 d1_r{i} = ops{i}.d1_r; 73 d1_r{i} = ops{i}.d1_r;
68 end 74 end
69 75
70 %====== Assemble full operators ======== 76 %====== Assemble full operators ========
71 A = spdiag(a); 77 A = spdiag(a);
78 obj.A = A;
72 RHO = spdiag(rho); 79 RHO = spdiag(rho);
80 obj.RHO = RHO;
81
73 82
74 obj.D1 = cell(dim,1); 83 obj.D1 = cell(dim,1);
75 obj.D2 = cell(dim,1); 84 obj.D2 = cell(dim,1);
76 obj.e_l = cell(dim,1); 85 obj.e_l = cell(dim,1);
77 obj.e_r = cell(dim,1); 86 obj.e_r = cell(dim,1);
78 obj.d1_l = cell(dim,1); 87 obj.d1_l = cell(dim,1);
79 obj.d1_r = cell(dim,1); 88 obj.d1_r = cell(dim,1);
80 89
81 % D1 90 % D1
82 obj.D1{1} = kron{D1{1},I(2)}; 91 obj.D1{1} = kron(D1{1},I{2});
83 obj.D1{2} = kron{I{1},D1(2)}; 92 obj.D1{2} = kron(I{1},D1{2});
84 93
85 % Boundary operators 94 % Boundary operators
86 obj.e_l{1} = kron{e_l{1},I(2)}; 95 obj.e_l{1} = kron(e_l{1},I{2});
87 obj.e_l{2} = kron{I{1},e_l(2)}; 96 obj.e_l{2} = kron(I{1},e_l{2});
88 obj.e_r{1} = kron{e_r{1},I(2)}; 97 obj.e_r{1} = kron(e_r{1},I{2});
89 obj.e_r{2} = kron{I{1},e_r(2)}; 98 obj.e_r{2} = kron(I{1},e_r{2});
90 99
91 obj.d1_l{1} = kron{d1_l{1},I(2)}; 100 obj.d1_l{1} = kron(d1_l{1},I{2});
92 obj.d1_l{2} = kron{I{1},d1_l(2)}; 101 obj.d1_l{2} = kron(I{1},d1_l{2});
93 obj.d1_r{1} = kron{d1_r{1},I(2)}; 102 obj.d1_r{1} = kron(d1_r{1},I{2});
94 obj.d1_r{2} = kron{I{1},d1_r(2)}; 103 obj.d1_r{2} = kron(I{1},d1_r{2});
95 104
96 % D2 105 % D2
97 for i = 1:dim 106 for i = 1:dim
98 obj.D2{i} = sparse(m_tot); 107 obj.D2{i} = sparse(m_tot);
99 end 108 end
111 obj.D2{2}(p,p) = D; 120 obj.D2{2}(p,p) = D;
112 end 121 end
113 122
114 % Quadratures 123 % Quadratures
115 obj.H = kron(H{1},H{2}); 124 obj.H = kron(H{1},H{2});
125 obj.Hi = inv(obj.H);
116 obj.H_boundary = cell(dim,1); 126 obj.H_boundary = cell(dim,1);
117 obj.H_boundary{1} = H{2}; 127 obj.H_boundary{1} = H{2};
118 obj.H_boundary{2} = H{1}; 128 obj.H_boundary{2} = H{1};
119 129
120 % Boundary coefficient matrices and quadratures 130 % Boundary coefficient matrices and quadratures
121 obj.A_boundary_l = cell(dim,1); 131 obj.A_boundary_l = cell(dim,1);
122 obj.A_boundary_r = cell(dim,1); 132 obj.A_boundary_r = cell(dim,1);
123 for i = 1:dim 133 for i = 1:dim
124 obj.A_boundary_l{i} = e_l{i}'*A*e_l{i}; 134 obj.A_boundary_l{i} = obj.e_l{i}'*A*obj.e_l{i};
125 obj.A_boundary_r{i} = e_r{i}'*A*e_r{i}; 135 obj.A_boundary_r{i} = obj.e_r{i}'*A*obj.e_r{i};
126 end 136 end
127 137
128 % E{i}^T picks out component i. 138 % E{i}^T picks out component i.
129 E = cell(dim,1); 139 E = cell(dim,1);
130 I = speye{mtot,mtot}; 140 I = speye(m_tot,m_tot);
131 for i = 1:dim 141 for i = 1:dim
132 e = sparse(dim,1); 142 e = sparse(dim,1);
133 e(i) = 1; 143 e(i) = 1;
134 E{i} = kron(e,I) 144 E{i} = kron(I,e);
135 end 145 end
136 obj.E = E; 146 obj.E = E;
137 147
138 % Differentiation matrix D (without SAT) 148 % Differentiation matrix D (without SAT)
139 D = 0; 149 D2 = obj.D2;
150 D1 = obj.D1;
151 D = sparse(dim*m_tot,dim*m_tot);
140 d = @kroneckerDelta; % Kronecker delta 152 d = @kroneckerDelta; % Kronecker delta
141 db = @(i,j) 1-dij(i,j); % Logical not of Kronecker delta 153 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta
142 for i = 1:dim 154 for i = 1:dim
143 for j = 1:dim 155 for j = 1:dim
144 D = D + E{i}*inv(rho)*( d(i,j)*D2{i}*E{j}' +... 156 D = D + E{i}*inv(RHO)*( d(i,j)*D2{i}*E{j}' +...
145 db(i,j)*D1{j}*A*D1{i}*E{j}' + ... 157 db(i,j)*D1{j}*A*D1{i}*E{j}' + ...
146 D2{j}*E{i}' ... 158 D2{j}*E{i}' ...
147 ); 159 );
148 end 160 end
149 end 161 end
150 obj.D = D; 162 obj.D = D;
151 %=========================================% 163 %=========================================%
152
153
154 164
155 % Misc. 165 % Misc.
156 obj.m = m; 166 obj.m = m;
157 obj.h = h; 167 obj.h = h;
158 obj.order = order; 168 obj.order = order;
159 obj.grid = g; 169 obj.grid = g;
160 170 obj.dim = dim;
161 obj.a = a;
162 obj.b = b;
163 171
164 % obj.gamm_u = h_u*ops_u.borrowing.M.d1; 172 % obj.gamm_u = h_u*ops_u.borrowing.M.d1;
165 % obj.gamm_v = h_v*ops_v.borrowing.M.d1; 173 % obj.gamm_v = h_v*ops_v.borrowing.M.d1;
166 end 174 end
167 175
176 function [closure, penalty] = boundary_condition(obj, boundary, type, parameter) 184 function [closure, penalty] = boundary_condition(obj, boundary, type, parameter)
177 default_arg('type','free'); 185 default_arg('type','free');
178 default_arg('parameter', []); 186 default_arg('parameter', []);
179 187
180 delta = @kroneckerDelta; % Kronecker delta 188 delta = @kroneckerDelta; % Kronecker delta
181 delta_b = @(i,j) 1-dij(i,j); % Logical not of Kronecker delta 189 delta_b = @(i,j) 1-delta(i,j); % Logical not of Kronecker delta
182 190
183 % j is the coordinate direction of the boundary 191 % j is the coordinate direction of the boundary
184 % nj: outward unit normal component. 192 % nj: outward unit normal component.
185 % nj = -1 for west, south, bottom boundaries 193 % nj = -1 for west, south, bottom boundaries
186 % nj = 1 for east, north, top boundaries 194 % nj = 1 for east, north, top boundaries
187 [j, nj] = obj.get_boundary_number(boundary); 195 [j, nj] = obj.get_boundary_number(boundary);
188 switch nj 196 switch nj
189 case 1 197 case 1
190 e = obj.e_r; 198 e = obj.e_r;
191 d = obj.d_r; 199 d = obj.d1_r;
192 case -1 200 case -1
193 e = obj.e_l; 201 e = obj.e_l;
194 d = obj.d_l; 202 d = obj.d1_l;
195 end 203 end
196 204
197 E = obj.E; 205 E = obj.E;
198 Hi = obj.Hi; 206 Hi = obj.Hi;
199 H_gamma = obj.H_boundary{j}; 207 H_gamma = obj.H_boundary{j};
200 A = obj.A; 208 A = obj.A;
209 RHO = obj.RHO;
210 D1 = obj.D1;
201 211
202 switch type 212 switch type
203 % Dirichlet boundary condition 213 % Dirichlet boundary condition
204 case {'D','d','dirichlet'} 214 case {'D','d','dirichlet'}
205 error('Dirichlet not implemented') 215 error('Dirichlet not implemented')
218 penalty = -obj.a*obj.Hi*tau; 228 penalty = -obj.a*obj.Hi*tau;
219 229
220 230
221 % Free boundary condition 231 % Free boundary condition
222 case {'F','f','Free','free'} 232 case {'F','f','Free','free'}
223 closure = 0; 233 closure = sparse(obj.dim*obj.grid.N,obj.dim*obj.grid.N);
224 penalty = 0; 234 penalty = sparse(obj.dim*obj.grid.N,obj.dim*obj.grid.N);
225 % Loop over components 235 % Loop over components
226 for i = 1:3 236 for i = 1:obj.dim
227 closure = closure + E{i}*(-sign)*Hi*e{j}*H_gamma*(... 237 closure = closure + E{i}*inv(RHO)*(-nj)*Hi*e{j}*H_gamma*(...
228 e{j}'*A*e{j}*d{j}'*E{i}' + ... 238 e{j}'*A*e{j}*d{j}'*E{i}' + ...
229 delta(i,j)*e{j}'*A*e{i}*d{i}*E{j}' + ... 239 delta(i,j)*e{j}'*A*e{i}*d{i}'*E{j}' + ...
230 delta_b(i,j)*A*D1{i}*E{j}' ... 240 delta_b(i,j)*e{j}'*A*D1{i}*E{j}' ...
231 ); 241 );
232 penalty = penalty - E{i}*(-sign)*Hi*e{j}*H_gamma; 242 penalty = penalty - E{i}*inv(RHO)*(-nj)*Hi*e{j}*H_gamma*e{j}'*E{j}';
233 end 243 end
244
234 245
235 % Unknown boundary condition 246 % Unknown boundary condition
236 otherwise 247 otherwise
237 error('No such boundary condition: type = %s',type); 248 error('No such boundary condition: type = %s',type);
238 end 249 end
244 tuning = 1.2; 255 tuning = 1.2;
245 % tuning = 20.2; 256 % tuning = 20.2;
246 error('Interface not implemented'); 257 error('Interface not implemented');
247 end 258 end
248 259
249 % Ruturns the coordinate number and outward normal component for the boundary specified by the string boundary. 260 % Returns the coordinate number and outward normal component for the boundary specified by the string boundary.
250 function [j, nj] = get_boundary_number(obj, boundary) 261 function [j, nj] = get_boundary_number(obj, boundary)
251 262
252 switch boundary 263 switch boundary
253 case {'w','W','west','West', 'e', 'E', 'east', 'East'} 264 case {'w','W','west','West', 'e', 'E', 'east', 'East'}
254 j = 1; 265 j = 1;
264 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} 275 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
265 nj = 1; 276 nj = 1;
266 end 277 end
267 end 278 end
268 279
280 % Returns the coordinate number and outward normal component for the boundary specified by the string boundary.
281 function [return_op] = get_boundary_operator(obj, op, boundary)
282
283 switch boundary
284 case {'w','W','west','West', 'e', 'E', 'east', 'East'}
285 j = 1;
286 case {'s','S','south','South', 'n', 'N', 'north', 'North'}
287 j = 2;
288 otherwise
289 error('No such boundary: boundary = %s',boundary);
290 end
291
292 switch op
293 case 'e'
294 switch boundary
295 case {'w','W','west','West','s','S','south','South'}
296 return_op = obj.e_l{j};
297 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
298 return_op = obj.e_r{j};
299 end
300 case 'd'
301 switch boundary
302 case {'w','W','west','West','s','S','south','South'}
303 return_op = obj.d_l{j};
304 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
305 return_op = obj.d_r{j};
306 end
307 otherwise
308 error(['No such operator: operatr = ' op]);
309 end
310
311 end
312
269 function N = size(obj) 313 function N = size(obj)
270 N = prod(obj.m); 314 N = prod(obj.m);
271 end 315 end
272 end 316 end
273 end 317 end