Mercurial > repos > public > sbplib

comparison +scheme/Staggered1DAcoustics.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 993aac771efd
children
comparison
equal deleted inserted replaced
672:82024d32e333 682:ec0ac76006e2
116 % type = 'characteristic' => incoming characteristc = data 116 % type = 'characteristic' => incoming characteristc = data
117 % No other types implemented yet 117 % No other types implemented yet
118 118
119 % BC on the form Lu - g = 0; 119 % BC on the form Lu - g = 0;
120 120
121 % Need a transformation T such that w = T^{−1}*u is a
122 % meaningful change of variables and T^T*B*T is block-diagonal,
123 % For linear acoustics, T = T_C meets these criteria.
124
121 A = obj.A; 125 A = obj.A;
122 B = obj.B; 126 B = obj.B;
123 C = inv(A)*B; 127 C = inv(A)*B;
124 128
125 % Diagonalize B 129 % Diagonalize C and use T_C to diagonalize B
126 [T, Lambda] = eig(B); 130 [T, ~] = eig(C);
131 Lambda = T'*B*T;
127 lambda = diag(Lambda); 132 lambda = diag(Lambda);
133
128 134
129 % Identify in- and outgoing variables 135 % Identify in- and outgoing variables
130 Iplus = lambda > 0; 136 Iplus = lambda > 0;
131 Iminus = lambda < 0; 137 Iminus = lambda < 0;
132 138
174 L = sparse(L); 180 L = sparse(L);
175 Tin = sparse(Tin); 181 Tin = sparse(Tin);
176 182
177 switch boundary 183 switch boundary
178 case 'l' 184 case 'l'
179 tau = -1*obj.e_l * inv(A) * T * sigma * inv(L*Tin); 185 tau = -1*obj.e_l * inv(A) * inv(T)' * sigma * inv(L*Tin);
180 closure = obj.Hi*tau*L*obj.e_l'; 186 closure = obj.Hi*tau*L*obj.e_l';
181 187
182 case 'r' 188 case 'r'
183 tau = 1*obj.e_r * inv(A) * T * sigma * inv(L*Tin); 189 tau = 1*obj.e_r * inv(A) * inv(T)' * sigma * inv(L*Tin);
184 closure = obj.Hi*tau*L*obj.e_r'; 190 closure = obj.Hi*tau*L*obj.e_r';
185 191
186 end 192 end
187 193
188 penalty = -obj.Hi*tau; 194 penalty = -obj.Hi*tau;