Mercurial > repos > public > sbplib
changeset 885:18e10217dca9 feature/d1_staggered
Generalize Staggered1DAcoustics to handle variable coefficients defined by both function handles and vectors.
| author | Martin Almquist <malmquist@stanford.edu> |
|---|---|
| date | Sun, 04 Nov 2018 12:36:30 -0800 |
| parents | 675e571e6f4a |
| children | c70131daaa6e |
| files | +scheme/Staggered1DAcousticsVariable.m |
| diffstat | 1 files changed, 67 insertions(+), 35 deletions(-) [+] |
line wrap: on
line diff
--- a/+scheme/Staggered1DAcousticsVariable.m Tue Jun 19 14:25:43 2018 -0700 +++ b/+scheme/Staggered1DAcousticsVariable.m Sun Nov 04 12:36:30 2018 -0800 @@ -2,10 +2,10 @@ properties m % Number of points of primal grid in each direction, possibly a vector h % Grid spacing - + % Grids grid % Total grid object - grid_primal + grid_primal grid_dual order % Order accuracy for the approximation @@ -39,15 +39,15 @@ methods - % Scheme for A*u_t + B u_x = 0, + % Scheme for A*u_t + B u_x = 0, % u = [p, v]; % A: Diagonal and A > 0, % A = [a1, 0; % 0, a2] - % + % % B = B^T and with diagonal entries = 0. - % B = [0 b - % b 0] + % B = [0 b + % b 0] % Here we store p on the primal grid and v on the dual % A and B are cell matrices of function handles function obj = Staggered1DAcousticsVariable(g, order, A, B) @@ -69,16 +69,48 @@ x_primal = obj.grid_primal.points(); x_dual = obj.grid_dual.points(); + % If coefficients are function handles, evaluate them + if isa(A{1,1}, 'function_handle') + A{1,1} = A{1,1}(x_primal); + A{1,2} = A{1,2}(x_dual); + A{2,1} = A{2,1}(x_primal); + A{2,2} = A{2,2}(x_dual); + end + + if isa(B{1,1}, 'function_handle') + B{1,1} = B{1,1}(x_primal); + B{1,2} = B{1,2}(x_primal); + B{2,1} = B{2,1}(x_dual); + B{2,2} = B{2,2}(x_dual); + end + + % If coefficents are vectors, turn them into diagonal matrices + [m, n] = size(A{1,1}); + if m==1 || n == 1 + A{1,1} = spdiag(A{1,1}, 0); + A{2,1} = spdiag(A{2,1}, 0); + A{1,2} = spdiag(A{1,2}, 0); + A{2,2} = spdiag(A{2,2}, 0); + end + + [m, n] = size(B{1,1}); + if m==1 || n == 1 + B{1,1} = spdiag(B{1,1}, 0); + B{2,1} = spdiag(B{2,1}, 0); + B{1,2} = spdiag(B{1,2}, 0); + B{2,2} = spdiag(B{2,2}, 0); + end + % Boundary matrices - obj.A_l = [A{1,1}(xl), A{1,2}(xl);.... - A{2,1}(xl), A{2,2}(xl)]; - obj.A_r = [A{1,1}(xr), A{1,2}(xr);.... - A{2,1}(xr), A{2,2}(xr)]; - obj.B_l = [B{1,1}(xl), B{1,2}(xl);.... - B{2,1}(xl), B{2,2}(xl)]; - obj.B_r = [B{1,1}(xr), B{1,2}(xr);.... - B{2,1}(xr), B{2,2}(xr)]; - + obj.A_l = full([A{1,1}(1,1), A{1,2}(1,1);.... + A{2,1}(1,1), A{2,2}(1,1)]); + obj.A_r = full([A{1,1}(end,end), A{1,2}(end,end);.... + A{2,1}(end,end), A{2,2}(end,end)]); + obj.B_l = full([B{1,1}(1,1), B{1,2}(1,1);.... + B{2,1}(1,1), B{2,2}(1,1)]); + obj.B_r = full([B{1,1}(end,end), B{1,2}(end,end);.... + B{2,1}(end,end), B{2,2}(end,end)]); + % Get operators ops = sbp.D1StaggeredUpwind(obj.m, xlim, order); obj.h = ops.h; @@ -88,14 +120,14 @@ H_dual = ops.H_dual; obj.H = blockmatrix.toMatrix( {H_primal, []; [], H_dual } ); obj.Hi = inv(obj.H); - + D1_primal = ops.D1_primal; D1_dual = ops.D1_dual; - A11_B12 = spdiag(-1./A{1,1}(x_primal).*B{1,2}(x_primal), 0); - A22_B21 = spdiag(-1./A{2,2}(x_dual).*B{2,1}(x_dual), 0); + A11_B12 = -A{1,1}\B{1,2}; + A22_B21 = -A{2,2}\B{2,1}; D = {[], A11_B12*D1_primal;... A22_B21*D1_dual, []}; - obj.D = blockmatrix.toMatrix(D); + obj.D = blockmatrix.toMatrix(D); obj.D1_primal = D1_primal; obj.D1_dual = D1_dual; @@ -104,9 +136,9 @@ e_primal_r = ops.e_primal_r; e_dual_l = ops.e_dual_l; e_dual_r = ops.e_dual_r; - e_l = {e_primal_l, [];... + e_l = {e_primal_l, [];... [] , e_dual_l}; - e_r = {e_primal_r, [];... + e_r = {e_primal_r, [];... [] , e_dual_r}; obj.e_l = blockmatrix.toMatrix(e_l); obj.e_r = blockmatrix.toMatrix(e_r); @@ -126,7 +158,7 @@ % neighbour_scheme is an instance of Scheme that should be interfaced to. % neighbour_boundary is a string specifying which boundary to interface to. function [closure, penalty] = boundary_condition(obj, boundary, type) - + default_arg('type','p'); % type = 'p' => boundary condition for p @@ -202,25 +234,25 @@ sigma = sparse(sigma); L = sparse(L); Tin = sparse(Tin); - + switch boundary case 'l' - tau = -1*obj.e_l * inv(A) * inv(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) * inv(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 - + penalty = -obj.Hi*tau; - + end - + % Uses the hat variable method for the interface coupling, % see hypsyst_varcoeff.pdf in the hypsyst repository. - % Notation: u for left side of interface, v for right side. + % Notation: u for left side of interface, v for right side. function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) % Get boundary matrices @@ -249,7 +281,7 @@ C_u = inv(A_u)*B_u; C_v = inv(A_v)*B_v; - % Diagonalize C + % Diagonalize C [T_u, ~] = eig(C_u); Lambda_u = inv(T_u)*C_u*T_u; lambda_u = diag(Lambda_u); @@ -279,18 +311,18 @@ % Create S_tilde and T_tilde Stilde = [Sp_u; Sm_v]; Ttilde = inv(Stilde); - Ttilde_p = Ttilde(:,1); + Ttilde_p = Ttilde(:,1); Ttilde_m = Ttilde(:,2); % Solve for penalty parameters theta_1,2 LHS = Ttilde_m*Sm_v*Tm_u; RHS = B_u*Tm_u; - th_u = RHS./LHS; + th_u = RHS./LHS; TH_u = diag(th_u); LHS = Ttilde_p*Sp_u*Tp_v; RHS = -B_v*Tp_v; - th_v = RHS./LHS; + th_v = RHS./LHS; TH_v = diag(th_v); % Construct penalty matrices @@ -313,9 +345,9 @@ closure = closure_u; penalty = penalty_u; end - + end - + function N = size(obj) N = 2*obj.m+1; end