Mercurial > repos > public > sbplib

comparison +scheme/Staggered1DAcoustics.m @ 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 8e55298657b9
children ec0ac76006e2
comparison
equal deleted inserted replaced
653:2351a7690e8a 671:993aac771efd
111 111
112 default_arg('type','p'); 112 default_arg('type','p');
113 113
114 % type = 'p' => boundary condition for p 114 % type = 'p' => boundary condition for p
115 % type = 'v' => boundary condition for v 115 % type = 'v' => boundary condition for v
116 % type = 'characteristic' => incoming characteristc = data
116 % No other types implemented yet 117 % No other types implemented yet
117 118
118 % BC on the form Lu - g = 0; 119 % BC on the form Lu - g = 0;
119 120
121 A = obj.A;
122 B = obj.B;
123 C = inv(A)*B;
124
120 % Diagonalize B 125 % Diagonalize B
121 B = obj.B;
122 [T, Lambda] = eig(B); 126 [T, Lambda] = eig(B);
123 lambda = diag(Lambda); 127 lambda = diag(Lambda);
124 128
125 % Identify in- and outgoing characteristic variables 129 % Identify in- and outgoing variables
126 Iplus = lambda > 0; 130 Iplus = lambda > 0;
127 Iminus = lambda < 0; 131 Iminus = lambda < 0;
128 132
129 switch boundary 133 switch boundary
130 case 'l' 134 case 'l'
131 Iout = Iminus;
132 Iin = Iplus; 135 Iin = Iplus;
133 case 'r' 136 case 'r'
134 Iout = Iplus;
135 Iin = Iminus; 137 Iin = Iminus;
136 end 138 end
137
138 Tin = T(:,Iin); 139 Tin = T(:,Iin);
139 Tout = T(:,Iout);
140 140
141 switch type 141 switch type
142 case 'p' 142 case 'p'
143 L = [1, 0]; 143 L = [1, 0];
144 case 'v' 144 case 'v'
145 L = [0, 1]; 145 L = [0, 1];
146 case 'characteristic' 146 case 'characteristic'
147 L = Tin'; 147 % Diagonalize C
148 [T_C, Lambda_C] = eig(C);
149 lambda_C = diag(Lambda_C);
150 % Identify in- and outgoing characteristic variables
151 Iplus = lambda_C > 0;
152 Iminus = lambda_C < 0;
153
154 switch boundary
155 case 'l'
156 Iin_C = Iplus;
157 case 'r'
158 Iin_C = Iminus;
159 end
160 T_C_inv = inv(T_C);
161 L = T_C_inv(Iin_C,:);
148 otherwise 162 otherwise
149 error('Boundary condition not implemented.'); 163 error('Boundary condition not implemented.');
150 end 164 end
151 165
152 % Penalty parameters 166 % Penalty parameters
153 A = obj.A;
154 sigma = [0; 0]; 167 sigma = [0; 0];
155 sigma(Iin) = lambda(Iin); 168 sigma(Iin) = lambda(Iin);
169
170 % Sparsify
171 A = sparse(A);
172 T = sparse(T);
173 sigma = sparse(sigma);
174 L = sparse(L);
175 Tin = sparse(Tin);
176
156 switch boundary 177 switch boundary
157 case 'l' 178 case 'l'
158 tau = -1*obj.e_l * inv(A) * T * sigma * inv(L*Tin); 179 tau = -1*obj.e_l * inv(A) * T * sigma * inv(L*Tin);
159 closure = obj.Hi*tau*L*obj.e_l'; 180 closure = obj.Hi*tau*L*obj.e_l';
160 181