Mercurial > repos > public > sbplib
diff +sbp/upwind8.m @ 47:ebd25e5a481a
Added upwind operators.
author | Jonatan Werpers <jonatan@werpers.com> |
---|---|
date | Thu, 05 Nov 2015 16:30:17 -0800 |
parents | |
children | 4f5a65f49035 |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/upwind8.m Thu Nov 05 16:30:17 2015 -0800 @@ -0,0 +1,24 @@ +function [H, HI, Dp, Dm, e_1, e_m] = upwind8(m,h) + H=spdiag(ones(m,1),0); + H(1:8,1:8)=[0.7489399e7 / 0.25401600e8 0 0 0 0 0 0 0; 0 0.5537831e7 / 0.3628800e7 0 0 0 0 0 0; 0 0 0.103373e6 / 0.403200e6 0 0 0 0 0; 0 0 0 0.261259e6 / 0.145152e6 0 0 0 0; 0 0 0 0 0.298231e6 / 0.725760e6 0 0 0; 0 0 0 0 0 0.515917e6 / 0.403200e6 0 0; 0 0 0 0 0 0 0.3349159e7 / 0.3628800e7 0; 0 0 0 0 0 0 0 0.25639991e8 / 0.25401600e8;]; + + H(m-7:m,m-7:m)=fliplr(flipud(H(1:8,1:8))); + H=H*h; + HI=inv(H); + + Qp=(-1/168*spdiag(ones(m-3,1),-3)+1/14*spdiag(ones(m-2,1),-2)-1/2*spdiag(ones(m-1,1),-1)-9/20*spdiag(ones(m,1),0)+5/4*spdiag(ones(m-1,1),+1)-1/2*spdiag(ones(m-2,1),+2)+1/6*spdiag(ones(m-3,1),+3)-1/28*spdiag(ones(m-4,1),+4)+1/280*spdiag(ones(m-5,1),+5)); + %Q_U = [-0.375523e6 / 0.237369744e9 0.73393569793e11 / 0.111599208000e12 -0.1244182786433e13 / 0.14954293872000e14 -0.245555527093e12 / 0.2492382312000e13 0.131952309637e12 / 0.14954293872000e14 0.62907543023e11 / 0.3738573468000e13; -0.4818422888131e13 / 0.7477146936000e13 -0.10098911e8 / 0.307701520e9 0.796148045807e12 / 0.1495429387200e13 0.487507785187e12 / 0.2990858774400e13 0.8887918703e10 / 0.249238231200e12 -0.804641655323e12 / 0.14954293872000e14; 0.925771462433e12 / 0.14954293872000e14 -0.601394962127e12 / 0.1495429387200e13 -0.674340047e9 / 0.4153970520e10 0.560418114917e12 / 0.747714693600e12 -0.1043208352973e13 / 0.2990858774400e13 0.42207619661e11 / 0.356054616000e12; 0.285728443093e12 / 0.2492382312000e13 -0.873874460227e12 / 0.2990858774400e13 -0.34960363091e11 / 0.106816384800e12 -0.1528791703e10 / 0.4153970520e10 0.1741615892927e13 / 0.1495429387200e13 -0.6107723941553e13 / 0.14954293872000e14; -0.213784143637e12 / 0.14954293872000e14 0.1011411511e10 / 0.35605461600e11 0.141331778093e12 / 0.2990858774400e13 -0.632571739487e12 / 0.1495429387200e13 -0.1464298201e10 / 0.2769313680e10 0.9523869708211e13 / 0.7477146936000e13; -0.60750297023e11 / 0.3738573468000e13 0.614597810123e12 / 0.14954293872000e14 -0.16980260027e11 / 0.2492382312000e13 -0.343200536047e12 / 0.14954293872000e14 -0.389757912373e12 / 0.1068163848000e13 -0.964053659e9 / 0.1661588208e10;]; + Q_U =[-0.16683e5 / 0.63018560e8 0.29345822969987e14 / 0.43354248537600e14 -0.2734625426591e13 / 0.40644608004000e14 -0.113480208109603e15 / 0.780376473676800e15 -0.830250230261e12 / 0.26012549122560e14 0.2500519492033e13 / 0.32515686403200e14 0.2235718279643e13 / 0.390188236838400e15 -0.388481888477e12 / 0.26543417472000e14; -0.29227665839987e14 / 0.43354248537600e14 -0.493793e6 / 0.63018560e8 0.8302717120817e13 / 0.26543417472000e14 0.3739408501537e13 / 0.9290196115200e13 0.2684481534461e13 / 0.13935294172800e14 -0.4450185662513e13 / 0.18580392230400e14 -0.1221838279381e13 / 0.37160784460800e14 0.90595000956023e14 / 0.1950941184192000e16; 0.2505689537591e13 / 0.40644608004000e14 -0.7312922392817e13 / 0.26543417472000e14 -0.69332623e8 / 0.1323389760e10 0.10994933811709e14 / 0.18580392230400e14 -0.9270952411151e13 / 0.18580392230400e14 0.3191238635141e13 / 0.20644880256000e14 0.4442211176987e13 / 0.92901961152000e14 -0.940661365031e12 / 0.32515686403200e14; 0.118016946570403e15 / 0.780376473676800e15 -0.4173878828737e13 / 0.9290196115200e13 -0.7990503962509e13 / 0.18580392230400e14 -0.207799621e9 / 0.1323389760e10 0.2044021560341e13 / 0.2477385630720e13 0.511197701761e12 / 0.18580392230400e14 0.1237681717213e13 / 0.13935294172800e14 -0.7784834666617e13 / 0.130062745612800e15; 0.68609076271e11 / 0.2364777192960e13 -0.2235651762161e13 / 0.13935294172800e14 0.6527681584751e13 / 0.18580392230400e14 -0.1115980068821e13 / 0.2477385630720e13 -0.55386253e8 / 0.189055680e9 0.3208334350649e13 / 0.3716078446080e13 -0.407569013461e12 / 0.844563283200e12 0.136474842626653e15 / 0.780376473676800e15; -0.2487637785013e13 / 0.32515686403200e14 0.4244231077313e13 / 0.18580392230400e14 -0.1550378843141e13 / 0.20644880256000e14 -0.5726967564961e13 / 0.18580392230400e14 -0.1017898941929e13 / 0.3716078446080e13 -0.526653889e9 / 0.1323389760e10 0.45241297077547e14 / 0.37160784460800e14 -0.2279608411897e13 / 0.5080576000500e13; -0.2164019088443e13 / 0.390188236838400e15 0.1263196075861e13 / 0.37160784460800e14 -0.6600697610987e13 / 0.92901961152000e14 0.556610591687e12 / 0.13935294172800e14 0.926842346471e12 / 0.9290196115200e13 -0.18757693936747e14 / 0.37160784460800e14 -0.584765899e9 / 0.1323389760e10 0.204646287449e12 / 0.168431424000e12; 0.387091928477e12 / 0.26543417472000e14 -0.90231551688023e14 / 0.1950941184192000e16 0.1032404418251e13 / 0.32515686403200e14 0.3502353445417e13 / 0.130062745612800e15 -0.15385068876253e14 / 0.780376473676800e15 0.262499068919e12 / 0.10161152001000e14 -0.867004691939e12 / 0.1852745664000e13 -0.85017967e8 / 0.189055680e9;]; + + Qp(1:8,1:8)=Q_U; + Qp(m-7:m,m-7:m)=flipud( fliplr(Q_U(1:8,1:8) ) )'; %%% This is different from standard SBP + + Qm=-Qp'; + + e_1=zeros(m,1);e_1(1)=1; + e_m=zeros(m,1);e_m(m)=1; + + Dp=HI*(Qp-1/2*e_1*e_1'+1/2*e_m*e_m') ; + + Dm=HI*(Qm-1/2*e_1*e_1'+1/2*e_m*e_m') ; +end