Mercurial > repos > public > sbplib
changeset 682:ec0ac76006e2 feature/d1_staggered
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 | 9c3c180c7ed3 |
| files | +scheme/Staggered1DAcoustics.m +scheme/Staggered1DAcousticsVariable.m |
| diffstat | 2 files changed, 19 insertions(+), 8 deletions(-) [+] |
line wrap: on
line diff
diff -r 82024d32e333 -r ec0ac76006e2 +scheme/Staggered1DAcoustics.m --- a/+scheme/Staggered1DAcoustics.m Wed Dec 27 21:16:06 2017 +0100 +++ b/+scheme/Staggered1DAcoustics.m Wed Feb 07 20:07:34 2018 -0800 @@ -118,14 +118,20 @@ % BC on the form Lu - g = 0; + % Need a transformation T such that w = T^{−1}*u is a + % meaningful change of variables and T^T*B*T is block-diagonal, + % For linear acoustics, T = T_C meets these criteria. + A = obj.A; B = obj.B; C = inv(A)*B; - % Diagonalize B - [T, Lambda] = eig(B); + % Diagonalize C and use T_C to diagonalize B + [T, ~] = eig(C); + Lambda = T'*B*T; lambda = diag(Lambda); + % Identify in- and outgoing variables Iplus = lambda > 0; Iminus = lambda < 0; @@ -176,11 +182,11 @@ switch boundary case 'l' - tau = -1*obj.e_l * inv(A) * T * sigma * inv(L*Tin); + tau = -1*obj.e_l * inv(A) * inv(T)' * sigma * inv(L*Tin); closure = obj.Hi*tau*L*obj.e_l'; case 'r' - tau = 1*obj.e_r * inv(A) * T * sigma * inv(L*Tin); + tau = 1*obj.e_r * inv(A) * inv(T)' * sigma * inv(L*Tin); closure = obj.Hi*tau*L*obj.e_r'; end
diff -r 82024d32e333 -r ec0ac76006e2 +scheme/Staggered1DAcousticsVariable.m --- a/+scheme/Staggered1DAcousticsVariable.m Wed Dec 27 21:16:06 2017 +0100 +++ b/+scheme/Staggered1DAcousticsVariable.m Wed Feb 07 20:07:34 2018 -0800 @@ -135,6 +135,10 @@ % BC on the form Lu - g = 0; + % Need a transformation T such that w = T^{−1}*u is a + % meaningful change of variables and T^T*B*T is block-diagonal, + % For linear acoustics, T = T_C meets these criteria. + % Get boundary matrices switch boundary case 'l' @@ -146,8 +150,9 @@ end C = inv(A)*B; - % Diagonalize B - [T, Lambda] = eig(B); + % Diagonalize C and use T_C to diagonalize B + [T, ~] = eig(C); + Lambda = T'*B*T; lambda = diag(Lambda); % Identify in- and outgoing characteristic variables @@ -200,11 +205,11 @@ switch boundary case 'l' - tau = -1*obj.e_l * inv(A) * T * sigma * inv(L*Tin); + tau = -1*obj.e_l * inv(A) * inv(T)' * sigma * inv(L*Tin); closure = obj.Hi*tau*L*obj.e_l'; case 'r' - tau = 1*obj.e_r * inv(A) * T * sigma * inv(L*Tin); + tau = 1*obj.e_r * inv(A) * inv(T)' * sigma * inv(L*Tin); closure = obj.Hi*tau*L*obj.e_r'; end