Mercurial > repos > public > sbplib
diff +sbp/upwind3_3bp.m @ 52:4f5a65f49035
Fixed the upwind operators.
| author | Jonatan Werpers <jonatan@werpers.com> |
|---|---|
| date | Wed, 11 Nov 2015 15:43:15 -0800 |
| parents | ebd25e5a481a |
| children |
line wrap: on
line diff
--- a/+sbp/upwind3_3bp.m Wed Nov 11 15:15:26 2015 -0800 +++ b/+sbp/upwind3_3bp.m Wed Nov 11 15:43:15 2015 -0800 @@ -1,22 +1,29 @@ function [H, HI, Dp, Dm, e_1, e_m] = upwind3(m,h) - H(1:3,1:3)=[0.3e1 / 0.8e1 0 0; 0 0.7e1 / 0.6e1 0; 0 0 0.23e2 / 0.24e2;]; - H(m-2:m,m-2:m)=fliplr(flipud(H(1:3,1:3))); - H=H*h; - H=spdiag(ones(m,1),0); - HI=inv(H); + Hv = ones(m,1); + Hv(1:3) = [3/8; 7/6; 23/24]; + Hv(m-2:m) = rot90(Hv(1:3),2); + Hv = Hv*h; + H = spdiag(Hv,0); + HI = spdiag(1./Hv,0); - Qp=(-1/3*spdiag(ones(m-1,1),-1)-1/2*spdiag(ones(m,1),0)+1*spdiag(ones(m-1,1),+1)-1/6*spdiag(ones(m-2,1),+2)); - %Q_U = [-0.49e2 / 0.2640e4 0.11167e5 / 0.15840e5 -0.1541e4 / 0.7920e4 0.43e2 / 0.5280e4; -0.9403e4 / 0.15840e5 -0.551e3 / 0.2640e4 0.1779e4 / 0.1760e4 -0.2311e4 / 0.7920e4; 0.659e3 / 0.7920e4 -0.751e3 / 0.1760e4 -0.1541e4 / 0.2640e4 0.21287e5 / 0.15840e5; 0.51e2 / 0.1760e4 -0.551e3 / 0.7920e4 -0.3683e4 / 0.15840e5 -0.713e3 / 0.880e3;]; - %Q_U = [-0.29e2 / 0.1680e4 0.7067e4 / 0.10080e5 -0.961e3 / 0.5040e4 0.23e2 / 0.3360e4; -0.6023e4 / 0.10080e5 -0.331e3 / 0.1680e4 0.1119e4 / 0.1120e4 -0.1451e4 / 0.5040e4; 0.439e3 / 0.5040e4 -0.491e3 / 0.1120e4 -0.961e3 / 0.1680e4 0.13507e5 / 0.10080e5; 0.31e2 / 0.1120e4 -0.331e3 / 0.5040e4 -0.2383e4 / 0.10080e5 -0.453e3 / 0.560e3;]; + q_diags = [-1, 0, 1, 2]; + q_stencil = [-1/3 -1/2 1 -1/6]; + Qp = stripeMatrix(q_stencil, q_diags,m); - Q_U = [-0.1e1 / 0.24e2 0.17e2 / 0.24e2 -0.1e1 / 0.6e1; -0.13e2 / 0.24e2 -0.1e1 / 0.4e1 0.23e2 / 0.24e2; 0.1e1 / 0.12e2 -0.11e2 / 0.24e2 -0.11e2 / 0.24e2;]; + Q_U = [ + -1/24 17/24 -1/6; + -13/24 -1/4 23/24; + 1/12 -11/24 -11/24; + ]; Qp(1:3,1:3)=Q_U; - Qp(m-2:m,m-2:m)=flipud( fliplr(Q_U(1:3,1:3) ) )'; %%% This is different from standard SBP + Qp(m-2:m,m-2:m)=rot90(Q_U,2)'; %%% This is different from standard SBP Qm=-Qp'; - e_1=sparse(m,1);e_1(1)=1; - e_m=sparse(m,1);e_m(m)=1; + e_1=sparse(m,1); + e_1(1)=1; + e_m=sparse(m,1); + e_m(m)=1; Dp=HI*(Qp-1/2*e_1*e_1'+1/2*e_m*e_m') ;