Mercurial > repos > public > sbplib

changeset 671:993aac771efd feature/d1_staggered

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Correct characteristic BC in staggered acoustics 1D constant coeff.
author Martin Almquist <malmquist@stanford.edu>
date Wed, 27 Dec 2017 17:18:51 +0100
parents 2351a7690e8a
children 82024d32e333
files +scheme/Staggered1DAcoustics.m
diffstat 1 files changed, 29 insertions(+), 8 deletions(-) [+]
line wrap: on
line diff
--- a/+scheme/Staggered1DAcoustics.m	Mon Dec 04 11:18:13 2017 -0800
+++ b/+scheme/Staggered1DAcoustics.m	Wed Dec 27 17:18:51 2017 +0100
@@ -113,30 +113,30 @@
 
             % type = 'p' => boundary condition for p
             % type = 'v' => boundary condition for v
+            % type = 'characteristic' => incoming characteristc = data
             % No other types implemented yet
 
             % BC on the form Lu - g = 0;
 
+            A = obj.A;
+            B = obj.B;
+            C = inv(A)*B;
+
             % Diagonalize B
-            B = obj.B;
             [T, Lambda] = eig(B);
             lambda = diag(Lambda);
 
-            % Identify in- and outgoing characteristic variables
+            % Identify in- and outgoing variables
             Iplus = lambda > 0;
             Iminus = lambda < 0;
 
             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'
@@ -144,15 +144,36 @@
             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
 
             % Penalty parameters
-            A = obj.A;
             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);