Mercurial > repos > public > sbplib

comparison +scheme/Staggered1DAcousticsVariable.m @ 885:18e10217dca9 feature/d1_staggered

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Generalize Staggered1DAcoustics to handle variable coefficients defined by both function handles and vectors.
author Martin Almquist <malmquist@stanford.edu>
date Sun, 04 Nov 2018 12:36:30 -0800
parents 8ed102db8e9c
children 9c74d96fc9e2
comparison
equal deleted inserted replaced
762:675e571e6f4a 885:18e10217dca9
1 classdef Staggered1DAcousticsVariable < scheme.Scheme 1 classdef Staggered1DAcousticsVariable < scheme.Scheme
2 properties 2 properties
3 m % Number of points of primal grid in each direction, possibly a vector 3 m % Number of points of primal grid in each direction, possibly a vector
4 h % Grid spacing 4 h % Grid spacing
5 5
6 % Grids 6 % Grids
7 grid % Total grid object 7 grid % Total grid object
8 grid_primal 8 grid_primal
9 grid_dual 9 grid_dual
10 10
11 order % Order accuracy for the approximation 11 order % Order accuracy for the approximation
12 12
13 H % Combined norm 13 H % Combined norm
37 B_l, B_r 37 B_l, B_r
38 end 38 end
39 39
40 40
41 methods 41 methods
42 % Scheme for A*u_t + B u_x = 0, 42 % Scheme for A*u_t + B u_x = 0,
43 % u = [p, v]; 43 % u = [p, v];
44 % A: Diagonal and A > 0, 44 % A: Diagonal and A > 0,
45 % A = [a1, 0; 45 % A = [a1, 0;
46 % 0, a2] 46 % 0, a2]
47 % 47 %
48 % B = B^T and with diagonal entries = 0. 48 % B = B^T and with diagonal entries = 0.
49 % B = [0 b 49 % B = [0 b
50 % b 0] 50 % b 0]
51 % Here we store p on the primal grid and v on the dual 51 % Here we store p on the primal grid and v on the dual
52 % A and B are cell matrices of function handles 52 % A and B are cell matrices of function handles
53 function obj = Staggered1DAcousticsVariable(g, order, A, B) 53 function obj = Staggered1DAcousticsVariable(g, order, A, B)
54 default_arg('B',{@(x)0*x, @(x)0*x+1; @(x)0*x+1, @(x)0*x}); 54 default_arg('B',{@(x)0*x, @(x)0*x+1; @(x)0*x+1, @(x)0*x});
55 default_arg('A',{@(x)0*x+1, @(x)0*x; @(x)0*x, @(x)0*x+1}); 55 default_arg('A',{@(x)0*x+1, @(x)0*x; @(x)0*x, @(x)0*x+1});
67 obj.grid_primal = g.grids{1}; 67 obj.grid_primal = g.grids{1};
68 obj.grid_dual = g.grids{2}; 68 obj.grid_dual = g.grids{2};
69 x_primal = obj.grid_primal.points(); 69 x_primal = obj.grid_primal.points();
70 x_dual = obj.grid_dual.points(); 70 x_dual = obj.grid_dual.points();
71 71
72 % If coefficients are function handles, evaluate them
73 if isa(A{1,1}, 'function_handle')
74 A{1,1} = A{1,1}(x_primal);
75 A{1,2} = A{1,2}(x_dual);
76 A{2,1} = A{2,1}(x_primal);
77 A{2,2} = A{2,2}(x_dual);
78 end
79
80 if isa(B{1,1}, 'function_handle')
81 B{1,1} = B{1,1}(x_primal);
82 B{1,2} = B{1,2}(x_primal);
83 B{2,1} = B{2,1}(x_dual);
84 B{2,2} = B{2,2}(x_dual);
85 end
86
87 % If coefficents are vectors, turn them into diagonal matrices
88 [m, n] = size(A{1,1});
89 if m==1 || n == 1
90 A{1,1} = spdiag(A{1,1}, 0);
91 A{2,1} = spdiag(A{2,1}, 0);
92 A{1,2} = spdiag(A{1,2}, 0);
93 A{2,2} = spdiag(A{2,2}, 0);
94 end
95
96 [m, n] = size(B{1,1});
97 if m==1 || n == 1
98 B{1,1} = spdiag(B{1,1}, 0);
99 B{2,1} = spdiag(B{2,1}, 0);
100 B{1,2} = spdiag(B{1,2}, 0);
101 B{2,2} = spdiag(B{2,2}, 0);
102 end
103
72 % Boundary matrices 104 % Boundary matrices
73 obj.A_l = [A{1,1}(xl), A{1,2}(xl);.... 105 obj.A_l = full([A{1,1}(1,1), A{1,2}(1,1);....
74 A{2,1}(xl), A{2,2}(xl)]; 106 A{2,1}(1,1), A{2,2}(1,1)]);
75 obj.A_r = [A{1,1}(xr), A{1,2}(xr);.... 107 obj.A_r = full([A{1,1}(end,end), A{1,2}(end,end);....
76 A{2,1}(xr), A{2,2}(xr)]; 108 A{2,1}(end,end), A{2,2}(end,end)]);
77 obj.B_l = [B{1,1}(xl), B{1,2}(xl);.... 109 obj.B_l = full([B{1,1}(1,1), B{1,2}(1,1);....
78 B{2,1}(xl), B{2,2}(xl)]; 110 B{2,1}(1,1), B{2,2}(1,1)]);
79 obj.B_r = [B{1,1}(xr), B{1,2}(xr);.... 111 obj.B_r = full([B{1,1}(end,end), B{1,2}(end,end);....
80 B{2,1}(xr), B{2,2}(xr)]; 112 B{2,1}(end,end), B{2,2}(end,end)]);
81 113
82 % Get operators 114 % Get operators
83 ops = sbp.D1StaggeredUpwind(obj.m, xlim, order); 115 ops = sbp.D1StaggeredUpwind(obj.m, xlim, order);
84 obj.h = ops.h; 116 obj.h = ops.h;
85 117
86 % Build combined operators 118 % Build combined operators
87 H_primal = ops.H_primal; 119 H_primal = ops.H_primal;
88 H_dual = ops.H_dual; 120 H_dual = ops.H_dual;
89 obj.H = blockmatrix.toMatrix( {H_primal, []; [], H_dual } ); 121 obj.H = blockmatrix.toMatrix( {H_primal, []; [], H_dual } );
90 obj.Hi = inv(obj.H); 122 obj.Hi = inv(obj.H);
91 123
92 D1_primal = ops.D1_primal; 124 D1_primal = ops.D1_primal;
93 D1_dual = ops.D1_dual; 125 D1_dual = ops.D1_dual;
94 A11_B12 = spdiag(-1./A{1,1}(x_primal).*B{1,2}(x_primal), 0); 126 A11_B12 = -A{1,1}\B{1,2};
95 A22_B21 = spdiag(-1./A{2,2}(x_dual).*B{2,1}(x_dual), 0); 127 A22_B21 = -A{2,2}\B{2,1};
96 D = {[], A11_B12*D1_primal;... 128 D = {[], A11_B12*D1_primal;...
97 A22_B21*D1_dual, []}; 129 A22_B21*D1_dual, []};
98 obj.D = blockmatrix.toMatrix(D); 130 obj.D = blockmatrix.toMatrix(D);
99 obj.D1_primal = D1_primal; 131 obj.D1_primal = D1_primal;
100 obj.D1_dual = D1_dual; 132 obj.D1_dual = D1_dual;
101 133
102 % Combined boundary operators 134 % Combined boundary operators
103 e_primal_l = ops.e_primal_l; 135 e_primal_l = ops.e_primal_l;
104 e_primal_r = ops.e_primal_r; 136 e_primal_r = ops.e_primal_r;
105 e_dual_l = ops.e_dual_l; 137 e_dual_l = ops.e_dual_l;
106 e_dual_r = ops.e_dual_r; 138 e_dual_r = ops.e_dual_r;
107 e_l = {e_primal_l, [];... 139 e_l = {e_primal_l, [];...
108 [] , e_dual_l}; 140 [] , e_dual_l};
109 e_r = {e_primal_r, [];... 141 e_r = {e_primal_r, [];...
110 [] , e_dual_r}; 142 [] , e_dual_r};
111 obj.e_l = blockmatrix.toMatrix(e_l); 143 obj.e_l = blockmatrix.toMatrix(e_l);
112 obj.e_r = blockmatrix.toMatrix(e_r); 144 obj.e_r = blockmatrix.toMatrix(e_r);
113 145
114 % Pick out first or second component of solution 146 % Pick out first or second component of solution
124 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. 156 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
125 % type is a string specifying the type of boundary condition if there are several. 157 % type is a string specifying the type of boundary condition if there are several.
126 % neighbour_scheme is an instance of Scheme that should be interfaced to. 158 % neighbour_scheme is an instance of Scheme that should be interfaced to.
127 % neighbour_boundary is a string specifying which boundary to interface to. 159 % neighbour_boundary is a string specifying which boundary to interface to.
128 function [closure, penalty] = boundary_condition(obj, boundary, type) 160 function [closure, penalty] = boundary_condition(obj, boundary, type)
129 161
130 default_arg('type','p'); 162 default_arg('type','p');
131 163
132 % type = 'p' => boundary condition for p 164 % type = 'p' => boundary condition for p
133 % type = 'v' => boundary condition for v 165 % type = 'v' => boundary condition for v
134 % No other types implemented yet 166 % No other types implemented yet
200 A = sparse(A); 232 A = sparse(A);
201 T = sparse(T); 233 T = sparse(T);
202 sigma = sparse(sigma); 234 sigma = sparse(sigma);
203 L = sparse(L); 235 L = sparse(L);
204 Tin = sparse(Tin); 236 Tin = sparse(Tin);
205 237
206 switch boundary 238 switch boundary
207 case 'l' 239 case 'l'
208 tau = -1*obj.e_l * inv(A) * inv(T)' * sigma * inv(L*Tin); 240 tau = -1*obj.e_l * inv(A) * inv(T)' * sigma * inv(L*Tin);
209 closure = obj.Hi*tau*L*obj.e_l'; 241 closure = obj.Hi*tau*L*obj.e_l';
210 242
211 case 'r' 243 case 'r'
212 tau = 1*obj.e_r * inv(A) * inv(T)' * sigma * inv(L*Tin); 244 tau = 1*obj.e_r * inv(A) * inv(T)' * sigma * inv(L*Tin);
213 closure = obj.Hi*tau*L*obj.e_r'; 245 closure = obj.Hi*tau*L*obj.e_r';
214 246
215 end 247 end
216 248
217 penalty = -obj.Hi*tau; 249 penalty = -obj.Hi*tau;
218 250
219 end 251 end
220 252
221 % Uses the hat variable method for the interface coupling, 253 % Uses the hat variable method for the interface coupling,
222 % see hypsyst_varcoeff.pdf in the hypsyst repository. 254 % see hypsyst_varcoeff.pdf in the hypsyst repository.
223 % Notation: u for left side of interface, v for right side. 255 % Notation: u for left side of interface, v for right side.
224 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) 256 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
225 257
226 % Get boundary matrices 258 % Get boundary matrices
227 switch boundary 259 switch boundary
228 case 'l' 260 case 'l'
247 e_v = neighbour_scheme.e_l; 279 e_v = neighbour_scheme.e_l;
248 end 280 end
249 C_u = inv(A_u)*B_u; 281 C_u = inv(A_u)*B_u;
250 C_v = inv(A_v)*B_v; 282 C_v = inv(A_v)*B_v;
251 283
252 % Diagonalize C 284 % Diagonalize C
253 [T_u, ~] = eig(C_u); 285 [T_u, ~] = eig(C_u);
254 Lambda_u = inv(T_u)*C_u*T_u; 286 Lambda_u = inv(T_u)*C_u*T_u;
255 lambda_u = diag(Lambda_u); 287 lambda_u = diag(Lambda_u);
256 S_u = inv(T_u); 288 S_u = inv(T_u);
257 289
277 Sm_v = S_v(Im_v,:); 309 Sm_v = S_v(Im_v,:);
278 310
279 % Create S_tilde and T_tilde 311 % Create S_tilde and T_tilde
280 Stilde = [Sp_u; Sm_v]; 312 Stilde = [Sp_u; Sm_v];
281 Ttilde = inv(Stilde); 313 Ttilde = inv(Stilde);
282 Ttilde_p = Ttilde(:,1); 314 Ttilde_p = Ttilde(:,1);
283 Ttilde_m = Ttilde(:,2); 315 Ttilde_m = Ttilde(:,2);
284 316
285 % Solve for penalty parameters theta_1,2 317 % Solve for penalty parameters theta_1,2
286 LHS = Ttilde_m*Sm_v*Tm_u; 318 LHS = Ttilde_m*Sm_v*Tm_u;
287 RHS = B_u*Tm_u; 319 RHS = B_u*Tm_u;
288 th_u = RHS./LHS; 320 th_u = RHS./LHS;
289 TH_u = diag(th_u); 321 TH_u = diag(th_u);
290 322
291 LHS = Ttilde_p*Sp_u*Tp_v; 323 LHS = Ttilde_p*Sp_u*Tp_v;
292 RHS = -B_v*Tp_v; 324 RHS = -B_v*Tp_v;
293 th_v = RHS./LHS; 325 th_v = RHS./LHS;
294 TH_v = diag(th_v); 326 TH_v = diag(th_v);
295 327
296 % Construct penalty matrices 328 % Construct penalty matrices
297 Z_u = TH_u*Ttilde_m*Sm_v; 329 Z_u = TH_u*Ttilde_m*Sm_v;
298 Z_u = sparse(Z_u); 330 Z_u = sparse(Z_u);
311 penalty = penalty_v; 343 penalty = penalty_v;
312 case 'r' 344 case 'r'
313 closure = closure_u; 345 closure = closure_u;
314 penalty = penalty_u; 346 penalty = penalty_u;
315 end 347 end
316 348
317 end 349 end
318 350
319 function N = size(obj) 351 function N = size(obj)
320 N = 2*obj.m+1; 352 N = 2*obj.m+1;
321 end 353 end
322 354
323 end 355 end