Mercurial > repos > public > sbplib

comparison +sbp/upwind6.m @ 52:4f5a65f49035

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Fixed the upwind operators.
author Jonatan Werpers <jonatan@werpers.com>
date Wed, 11 Nov 2015 15:43:15 -0800
parents ebd25e5a481a
children
comparison
equal deleted inserted replaced
51:0be702829bb9 52:4f5a65f49035
1 function [H, HI, Dp, Dm, e_1, e_m] = upwind6(m,h) 1 function [H, HI, Dp, Dm, e_1, e_m] = upwind6(m,h)
2 H=spdiag(ones(m,1),0); 2 Hv = ones(m,1);
3 H(1:6,1:6)=[0.13613e5 / 0.43200e5 0 0 0 0 0; 0 0.12049e5 / 0.8640e4 0 0 0 0; 0 0 0.535e3 / 0.864e3 0 0 0; 0 0 0 0.1079e4 / 0.864e3 0 0; 0 0 0 0 0.7841e4 / 0.8640e4 0; 0 0 0 0 0 0.43837e5 / 0.43200e5;]; 3 Hv(1:6)=[0.13613e5/0.43200e5; 0.12049e5/0.8640e4 ; 0.535e3/0.864e3 ; 0.1079e4/0.864e3 ; 0.7841e4/0.8640e4 ; 0.43837e5/0.43200e5];
4 H(m-5:m,m-5:m)=fliplr(flipud(H(1:6,1:6))); 4 Hv(m-5:m)=rot90(Hv(1:6),2);
5 H=H*h; 5 Hv = Hv*h;
6 HI=inv(H); 6 H = spdiag(Hv,0);
7 HI = spdiag(1./Hv,0);
7 8
9 q_diags = [-2 -1 0 1 2 3 4];
10 q_stencil = [1/30 -2/5 -7/12 4/3 -1/2 2/15 -1/60];
11 Qp = stripeMatrix(q_stencil, q_diags,m);
8 12
9 Qp=(1/30*spdiag(ones(m-2,1),-2)-2/5*spdiag(ones(m-1,1),-1)-7/12*spdiag(ones(m,1),0)+4/3*spdiag(ones(m-1,1),+1)-1/2*spdiag(ones(m-2,1),+2)+2/15*spdiag(ones(m-3,1),+3)-1/60*spdiag(ones(m-4,1),+4)); 13 Q_U =[
10 %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;]; 14 -0.265e3/0.128688e6 0.1146190567e10/0.1737288000e10 -0.1596619e7/0.18384000e8 -0.55265831e8/0.579096000e9 0.26269819e8/0.3474576000e10 0.2464501e7/0.144774000e9;
11 Q_U =[-0.265e3 / 0.128688e6 0.1146190567e10 / 0.1737288000e10 -0.1596619e7 / 0.18384000e8 -0.55265831e8 / 0.579096000e9 0.26269819e8 / 0.3474576000e10 0.2464501e7 / 0.144774000e9; -0.1116490567e10 / 0.1737288000e10 -0.8839e4 / 0.214480e6 0.190538869e9 / 0.347457600e9 0.102705469e9 / 0.694915200e9 0.413741e6 / 0.9651600e7 -0.191689861e9 / 0.3474576000e10; 0.1096619e7 / 0.18384000e8 -0.135385429e9 / 0.347457600e9 -0.61067e5 / 0.321720e6 0.45137333e8 / 0.57909600e8 -0.253641811e9 / 0.694915200e9 0.70665929e8 / 0.579096000e9; 0.66965831e8 / 0.579096000e9 -0.208765789e9 / 0.694915200e9 -0.17623253e8 / 0.57909600e8 -0.18269e5 / 0.45960e5 0.410905829e9 / 0.347457600e9 -0.477953317e9 / 0.1158192000e10; -0.49219819e8 / 0.3474576000e10 0.293299e6 / 0.9651600e7 0.26422771e8 / 0.694915200e9 -0.141938309e9 / 0.347457600e9 -0.346583e6 / 0.643440e6 0.2217185207e10 / 0.1737288000e10; -0.2374501e7 / 0.144774000e9 0.142906261e9 / 0.3474576000e10 -0.3137129e7 / 0.579096000e9 -0.29884283e8 / 0.1158192000e10 -0.630168407e9 / 0.1737288000e10 -0.3559e4 / 0.6128e4;]; 15 -0.1116490567e10/0.1737288000e10 -0.8839e4/0.214480e6 0.190538869e9/0.347457600e9 0.102705469e9/0.694915200e9 0.413741e6/0.9651600e7 -0.191689861e9/0.3474576000e10;
16 0.1096619e7/0.18384000e8 -0.135385429e9/0.347457600e9 -0.61067e5/0.321720e6 0.45137333e8/0.57909600e8 -0.253641811e9/0.694915200e9 0.70665929e8/0.579096000e9;
17 0.66965831e8/0.579096000e9 -0.208765789e9/0.694915200e9 -0.17623253e8/0.57909600e8 -0.18269e5/0.45960e5 0.410905829e9/0.347457600e9 -0.477953317e9/0.1158192000e10;
18 -0.49219819e8/0.3474576000e10 0.293299e6/0.9651600e7 0.26422771e8/0.694915200e9 -0.141938309e9/0.347457600e9 -0.346583e6/0.643440e6 0.2217185207e10/0.1737288000e10;
19 -0.2374501e7/0.144774000e9 0.142906261e9/0.3474576000e10 -0.3137129e7/0.579096000e9 -0.29884283e8/0.1158192000e10 -0.630168407e9/0.1737288000e10 -0.3559e4/0.6128e4;
20 ];
12 21
13 Qp(1:6,1:6)=Q_U; 22 Qp(1:6,1:6)=Q_U;
14 Qp(m-5:m,m-5:m)=flipud( fliplr(Q_U(1:6,1:6) ) )'; %%% This is different from standard SBP 23 Qp(m-5:m,m-5:m)=rot90(Q_U,2)'; %%% This is different from standard SBP
15 24
16 Qm=-Qp'; 25 Qm=-Qp';
17 26
18 e_1=zeros(m,1);e_1(1)=1; 27 e_1=sparse(m,1);e_1(1)=1;
19 e_m=zeros(m,1);e_m(m)=1; 28 e_m=sparse(m,1);e_m(m)=1;
20 29
21 Dp=HI*(Qp-1/2*e_1*e_1'+1/2*e_m*e_m') ; 30 Dp=HI*(Qp-1/2*e_1*e_1'+1/2*e_m*e_m') ;
22 31
23 Dm=HI*(Qm-1/2*e_1*e_1'+1/2*e_m*e_m') ; 32 Dm=HI*(Qm-1/2*e_1*e_1'+1/2*e_m*e_m') ;
24 33