Mercurial > repos > public > sbplib

annotate +sbp/upwind4.m @ 47:ebd25e5a481a

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Added upwind operators.
author Jonatan Werpers <jonatan@werpers.com>
date Thu, 05 Nov 2015 16:30:17 -0800
parents
children 4f5a65f49035
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
47
ebd25e5a481a Added upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
1 function [H, HI, Dp, Dm, e_1, e_m] = upwind4(m,h)
ebd25e5a481a Added upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
2
ebd25e5a481a Added upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
3 H=spdiag(ones(m,1),0);
ebd25e5a481a Added upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
4 H(1:4,1:4)=spdiag([49/144 61/48 41/48 149/144]);
ebd25e5a481a Added upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
5 H(m-3:m,m-3:m)=fliplr(flipud(spdiag([49/144 61/48 41/48 149/144])));
ebd25e5a481a Added upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
6 H=H*h;
ebd25e5a481a Added upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
7 HI=inv(H);
ebd25e5a481a Added upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
8
ebd25e5a481a Added upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
9
ebd25e5a481a Added upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
10 Qp=(-1/4*spdiag(ones(m-1,1),-1)-5/6*spdiag(ones(m,1),0)+3/2*spdiag(ones(m-1,1),+1)-1/2*spdiag(ones(m-2,1),+2)+1/12*spdiag(ones(m-3,1),+3));
ebd25e5a481a Added upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
11 %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;];
ebd25e5a481a Added upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
12 %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;];
ebd25e5a481a Added upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
13
ebd25e5a481a Added upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
14 Q_U = [-0.1e1 / 0.48e2 0.205e3 / 0.288e3 -0.29e2 / 0.144e3 0.1e1 / 0.96e2; -0.169e3 / 0.288e3 -0.11e2 / 0.48e2 0.33e2 / 0.32e2 -0.43e2 / 0.144e3; 0.11e2 / 0.144e3 -0.13e2 / 0.32e2 -0.29e2 / 0.48e2 0.389e3 / 0.288e3; 0.1e1 / 0.32e2 -0.11e2 / 0.144e3 -0.65e2 / 0.288e3 -0.13e2 / 0.16e2;];
ebd25e5a481a Added upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
15
ebd25e5a481a Added upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
16 Qp(1:4,1:4)=Q_U;
ebd25e5a481a Added upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
17 Qp(m-3:m,m-3:m)=flipud( fliplr(Q_U(1:4,1:4) ) )'; %%% This is different from standard SBP
ebd25e5a481a Added upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
18
ebd25e5a481a Added upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
19 Qm=-Qp';
ebd25e5a481a Added upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
20
ebd25e5a481a Added upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
21 e_1=zeros(m,1);e_1(1)=1;
ebd25e5a481a Added upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
22 e_m=zeros(m,1);e_m(m)=1;
ebd25e5a481a Added upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
23
ebd25e5a481a Added upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
24 Dp=HI*(Qp-1/2*e_1*e_1'+1/2*e_m*e_m') ;
ebd25e5a481a Added upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
25
ebd25e5a481a Added upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
26 Dm=HI*(Qm-1/2*e_1*e_1'+1/2*e_m*e_m') ;
ebd25e5a481a Added upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
27 end