Mercurial > repos > public > sbplib

comparison +scheme/Staggered1DAcousticsVariable.m @ 761:8ed102db8e9c feature/d1_staggered

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add discont interface in 1D acoustics staggered, using the hat variable interface coupling.
author Martin Almquist <malmquist@stanford.edu>
date Mon, 18 Jun 2018 16:36:24 -0700
parents ec0ac76006e2
children 18e10217dca9
comparison
equal deleted inserted replaced
760:408fb46f7266 761:8ed102db8e9c
216 216
217 penalty = -obj.Hi*tau; 217 penalty = -obj.Hi*tau;
218 218
219 end 219 end
220 220
221 % Uses the hat variable method for the interface coupling,
222 % see hypsyst_varcoeff.pdf in the hypsyst repository.
223 % Notation: u for left side of interface, v for right side.
221 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) 224 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
222 225
223 error('Staggered1DAcoustics, interface not implemented'); 226 % Get boundary matrices
227 switch boundary
228 case 'l'
229 A_v = obj.A_l;
230 B_v = obj.B_l;
231 Hi_v = obj.Hi;
232 e_v = obj.e_l;
233
234 A_u = neighbour_scheme.A_r;
235 B_u = neighbour_scheme.B_r;
236 Hi_u = neighbour_scheme.Hi;
237 e_u = neighbour_scheme.e_r;
238 case 'r'
239 A_u = obj.A_r;
240 B_u = obj.B_r;
241 Hi_u = obj.Hi;
242 e_u = obj.e_r;
243
244 A_v = neighbour_scheme.A_l;
245 B_v = neighbour_scheme.B_l;
246 Hi_v = neighbour_scheme.Hi;
247 e_v = neighbour_scheme.e_l;
248 end
249 C_u = inv(A_u)*B_u;
250 C_v = inv(A_v)*B_v;
251
252 % Diagonalize C
253 [T_u, ~] = eig(C_u);
254 Lambda_u = inv(T_u)*C_u*T_u;
255 lambda_u = diag(Lambda_u);
256 S_u = inv(T_u);
257
258 [T_v, ~] = eig(C_v);
259 Lambda_v = inv(T_v)*C_v*T_v;
260 lambda_v = diag(Lambda_v);
261 S_v = inv(T_v);
262
263 % Identify in- and outgoing characteristic variables
264 Im_u = lambda_u < 0;
265 Ip_u = lambda_u > 0;
266 Im_v = lambda_v < 0;
267 Ip_v = lambda_v > 0;
268
269 Tp_u = T_u(:,Ip_u);
270 Tm_u = T_u(:,Im_u);
271 Sp_u = S_u(Ip_u,:);
272 Sm_u = S_u(Im_u,:);
273
274 Tp_v = T_v(:,Ip_v);
275 Tm_v = T_v(:,Im_v);
276 Sp_v = S_v(Ip_v,:);
277 Sm_v = S_v(Im_v,:);
278
279 % Create S_tilde and T_tilde
280 Stilde = [Sp_u; Sm_v];
281 Ttilde = inv(Stilde);
282 Ttilde_p = Ttilde(:,1);
283 Ttilde_m = Ttilde(:,2);
284
285 % Solve for penalty parameters theta_1,2
286 LHS = Ttilde_m*Sm_v*Tm_u;
287 RHS = B_u*Tm_u;
288 th_u = RHS./LHS;
289 TH_u = diag(th_u);
290
291 LHS = Ttilde_p*Sp_u*Tp_v;
292 RHS = -B_v*Tp_v;
293 th_v = RHS./LHS;
294 TH_v = diag(th_v);
295
296 % Construct penalty matrices
297 Z_u = TH_u*Ttilde_m*Sm_v;
298 Z_u = sparse(Z_u);
299
300 Z_v = TH_v*Ttilde_p*Sp_u;
301 Z_v = sparse(Z_v);
302
303 closure_u = Hi_u*e_u*Z_u*e_u';
304 penalty_u = -Hi_u*e_u*Z_u*e_v';
305 closure_v = Hi_v*e_v*Z_v*e_v';
306 penalty_v = -Hi_v*e_v*Z_v*e_u';
224 307
225 switch boundary 308 switch boundary
226 % Upwind coupling 309 case 'l'
227 case {'l','left'} 310 closure = closure_v;
228 tau = -1*obj.e_l; 311 penalty = penalty_v;
229 closure = obj.Hi*tau*obj.e_l'; 312 case 'r'
230 penalty = -obj.Hi*tau*neighbour_scheme.e_r'; 313 closure = closure_u;
231 case {'r','right'} 314 penalty = penalty_u;
232 tau = 0*obj.e_r;
233 closure = obj.Hi*tau*obj.e_r';
234 penalty = -obj.Hi*tau*neighbour_scheme.e_l';
235 end 315 end
236 316
237 end 317 end
238 318
239 function N = size(obj) 319 function N = size(obj)
240 N = obj.m; 320 N = 2*obj.m+1;
241 end 321 end
242 322
243 end 323 end
244 324
245 methods(Static) 325 methods(Static)