Mercurial > repos > public > sbplib

comparison +sbp/+implementations/d1_upwind_4.m @ 423:a2cb0d4f4a02 feature/grids

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge in default.
author Jonatan Werpers <jonatan@werpers.com>
date Tue, 07 Feb 2017 15:47:51 +0100
parents f7ac3cd6eeaa
children 88584b0cfba1
comparison
equal deleted inserted replaced
218:da058ce66876 423:a2cb0d4f4a02
1 function [H, HI, Dp, Dm, e_1, e_m] = d1_upwind_4(m,h)
2
3 if(m<8)
4 error('Operator requires at least 8 grid points');
5 end
6
7 Hv=ones(m,1);
8 Hv(1:4)=[49/144 61/48 41/48 149/144];
9 Hv(m-3:m)=rot90(Hv(1:4),2);
10 Hv = Hv*h;
11 H = spdiag(Hv,0);
12 HI = spdiag(1./Hv,0);
13
14 q_diags = [-1, 0, 1, 2 3];
15 q_stencil = [-1/4 -5/6 3/2 -1/2 1/12];
16 Qp = stripeMatrix(q_stencil, q_diags,m);
17
18 Q_U = [
19 -0.1e1/0.48e2 0.205e3/0.288e3 -0.29e2/0.144e3 0.1e1/0.96e2;
20 -0.169e3/0.288e3 -0.11e2/0.48e2 0.33e2/0.32e2 -0.43e2/0.144e3;
21 0.11e2/0.144e3 -0.13e2/0.32e2 -0.29e2/0.48e2 0.389e3/0.288e3;
22 0.1e1/0.32e2 -0.11e2/0.144e3 -0.65e2/0.288e3 -0.13e2/0.16e2;
23 ];
24
25 Qp(1:4,1:4)=Q_U;
26 Qp(m-3:m,m-3:m)=rot90(Q_U,2)'; %%% This is different from standard SBP
27
28 Qm=-Qp';
29
30 e_1=sparse(m,1);e_1(1)=1;
31 e_m=sparse(m,1);e_m(m)=1;
32
33 Dp=HI*(Qp-1/2*(e_1*e_1')+1/2*(e_m*e_m')) ;
34
35 Dm=HI*(Qm-1/2*(e_1*e_1')+1/2*(e_m*e_m')) ;
36 end