Mercurial > repos > public > sbplib

comparison +sbp/+implementations/d1_upwind_2.m @ 284:dae8c3a56f5e

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merged in operator_remake (pull request #2) Operator remake
author Jonatan Werpers <jonatan.werpers@it.uu.se>
date Mon, 12 Sep 2016 12:53:02 +0200
parents f7ac3cd6eeaa
children
comparison
equal deleted inserted replaced
282:18c023aaf3f7 284:dae8c3a56f5e
1 function [H, HI, Dp, Dm, e_1, e_m] = d1_upwind_2(m,h)
2
3 if(m<4)
4 error('Operator requires at least 4 grid points');
5 end
6
7 Hv=ones(m,1);
8 Hv(1:2)=[0.1e1/0.4e1; 0.5e1/0.4e1;];
9 Hv(m-1:m)=rot90(Hv(1:2),2);
10 Hv = Hv*h;
11 H = spdiag(Hv,0);
12 HI = spdiag(1./Hv,0);
13
14 q_diags = [0 1 2];
15 q_stencil = [-3/2 +2 -1/2];
16 Qp = stripeMatrix(q_stencil, q_diags,m);
17
18 Q_U = [
19 -0.1e1/0.4e1 0.5e1/0.4e1;
20 -0.1e1/0.4e1 -0.5e1/0.4e1;
21 ];
22
23 Qp(1:2,1:2)=Q_U;
24 Qp(m-1:m,m-1:m)=rot90(Q_U,2)'; %%% This is different from standard SBP
25
26 Qm=-Qp';
27
28 e_1=sparse(m,1);e_1(1)=1;
29 e_m=sparse(m,1);e_m(m)=1;
30
31 Dp=HI*(Qp-1/2*(e_1*e_1')+1/2*(e_m*e_m')) ;
32
33 Dm=HI*(Qm-1/2*(e_1*e_1')+1/2*(e_m*e_m')) ;
34 end
35