changeset 672:82024d32e333 feature/d1_staggered

Correct char BC in staggered acoustics 1D variable coeff.
author Martin Almquist <malmquist@stanford.edu>
date Wed, 27 Dec 2017 21:16:06 +0100
parents 993aac771efd
children ec0ac76006e2
files +scheme/Staggered1DAcousticsVariable.m
diffstat 1 files changed, 24 insertions(+), 5 deletions(-) [+]
line wrap: on
line diff
--- a/+scheme/Staggered1DAcousticsVariable.m	Wed Dec 27 17:18:51 2017 +0100
+++ b/+scheme/Staggered1DAcousticsVariable.m	Wed Dec 27 21:16:06 2017 +0100
@@ -144,6 +144,7 @@
                     A = obj.A_r;
                     B = obj.B_r;
             end
+            C = inv(A)*B;
 
             % Diagonalize B
             [T, Lambda] = eig(B);
@@ -155,15 +156,11 @@
 
             switch boundary
             case 'l'
-                Iout = Iminus;
                 Iin = Iplus;
             case 'r'
-                Iout = Iplus;
                 Iin = Iminus;
             end
-
             Tin = T(:,Iin);
-            Tout = T(:,Iout);
 
             switch type
             case 'p'
@@ -171,7 +168,21 @@
             case 'v'
                 L = [0, 1];
             case 'characteristic'
-                L = Tin';
+                % Diagonalize C
+                [T_C, Lambda_C] = eig(C);
+                lambda_C = diag(Lambda_C);
+                % Identify in- and outgoing characteristic variables
+                Iplus = lambda_C > 0;
+                Iminus = lambda_C < 0;
+
+                switch boundary
+                case 'l'
+                    Iin_C = Iplus;
+                case 'r'
+                    Iin_C = Iminus;
+                end
+                T_C_inv = inv(T_C);
+                L = T_C_inv(Iin_C,:);
             otherwise
                 error('Boundary condition not implemented.');
             end
@@ -179,6 +190,14 @@
             % Penalty parameters
             sigma = [0; 0];
             sigma(Iin) = lambda(Iin);
+
+            % Sparsify
+            A = sparse(A);
+            T = sparse(T);
+            sigma = sparse(sigma);
+            L = sparse(L);
+            Tin = sparse(Tin);
+            
             switch boundary
             case 'l'
                 tau = -1*obj.e_l * inv(A) * T * sigma * inv(L*Tin);