Mercurial > repos > public > sbplib

comparison +scheme/Staggered1DAcousticsVariable.m @ 682:ec0ac76006e2 feature/d1_staggered

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Modify BC in Staggered 1D acoustics to ensure that only physically meaningful changes of variables are used. See hypsyst varcoeff sec 3.
author Martin Almquist <malmquist@stanford.edu>
date Wed, 07 Feb 2018 20:07:34 -0800
parents 82024d32e333
children 8ed102db8e9c
comparison
equal deleted inserted replaced
672:82024d32e333 682:ec0ac76006e2
133 % type = 'v' => boundary condition for v 133 % type = 'v' => boundary condition for v
134 % No other types implemented yet 134 % No other types implemented yet
135 135
136 % BC on the form Lu - g = 0; 136 % BC on the form Lu - g = 0;
137 137
138 % Need a transformation T such that w = T^{−1}*u is a
139 % meaningful change of variables and T^T*B*T is block-diagonal,
140 % For linear acoustics, T = T_C meets these criteria.
141
138 % Get boundary matrices 142 % Get boundary matrices
139 switch boundary 143 switch boundary
140 case 'l' 144 case 'l'
141 A = obj.A_l; 145 A = obj.A_l;
142 B = obj.B_l; 146 B = obj.B_l;
144 A = obj.A_r; 148 A = obj.A_r;
145 B = obj.B_r; 149 B = obj.B_r;
146 end 150 end
147 C = inv(A)*B; 151 C = inv(A)*B;
148 152
149 % Diagonalize B 153 % Diagonalize C and use T_C to diagonalize B
150 [T, Lambda] = eig(B); 154 [T, ~] = eig(C);
155 Lambda = T'*B*T;
151 lambda = diag(Lambda); 156 lambda = diag(Lambda);
152 157
153 % Identify in- and outgoing characteristic variables 158 % Identify in- and outgoing characteristic variables
154 Iplus = lambda > 0; 159 Iplus = lambda > 0;
155 Iminus = lambda < 0; 160 Iminus = lambda < 0;
198 L = sparse(L); 203 L = sparse(L);
199 Tin = sparse(Tin); 204 Tin = sparse(Tin);
200 205
201 switch boundary 206 switch boundary
202 case 'l' 207 case 'l'
203 tau = -1*obj.e_l * inv(A) * T * sigma * inv(L*Tin); 208 tau = -1*obj.e_l * inv(A) * inv(T)' * sigma * inv(L*Tin);
204 closure = obj.Hi*tau*L*obj.e_l'; 209 closure = obj.Hi*tau*L*obj.e_l';
205 210
206 case 'r' 211 case 'r'
207 tau = 1*obj.e_r * inv(A) * T * sigma * inv(L*Tin); 212 tau = 1*obj.e_r * inv(A) * inv(T)' * sigma * inv(L*Tin);
208 closure = obj.Hi*tau*L*obj.e_r'; 213 closure = obj.Hi*tau*L*obj.e_r';
209 214
210 end 215 end
211 216
212 penalty = -obj.Hi*tau; 217 penalty = -obj.Hi*tau;