Mercurial > repos > public > sbplib

annotate +sbp/upwind4.m @ 238:5f75d7094983 feature/beams

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
grid: Fixed returns of funcToXXX
author Jonatan Werpers <jonatan@werpers.com>
date Thu, 28 Jul 2016 11:45:02 +0200
parents 4f5a65f49035
children
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
52
4f5a65f49035 Fixed the upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents: 47
diff changeset
3 Hv=ones(m,1);
4f5a65f49035 Fixed the upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents: 47
diff changeset
4 Hv(1:4)=[49/144 61/48 41/48 149/144];
4f5a65f49035 Fixed the upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents: 47
diff changeset
5 Hv(m-3:m)=rot90(Hv(1:4),2);
4f5a65f49035 Fixed the upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents: 47
diff changeset
6 Hv = Hv*h;
4f5a65f49035 Fixed the upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents: 47
diff changeset
7 H = spdiag(Hv,0);
4f5a65f49035 Fixed the upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents: 47
diff changeset
8 HI = spdiag(1./Hv,0);
47
ebd25e5a481a Added upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
9
52
4f5a65f49035 Fixed the upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents: 47
diff changeset
10 q_diags = [-1, 0, 1, 2 3];
4f5a65f49035 Fixed the upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents: 47
diff changeset
11 q_stencil = [-1/4 -5/6 3/2 -1/2 1/12];
4f5a65f49035 Fixed the upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents: 47
diff changeset
12 Qp = stripeMatrix(q_stencil, q_diags,m);
47
ebd25e5a481a Added upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
13
52
4f5a65f49035 Fixed the upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents: 47
diff changeset
14 Q_U = [
4f5a65f49035 Fixed the upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents: 47
diff changeset
15 -0.1e1/0.48e2 0.205e3/0.288e3 -0.29e2/0.144e3 0.1e1/0.96e2;
4f5a65f49035 Fixed the upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents: 47
diff changeset
16 -0.169e3/0.288e3 -0.11e2/0.48e2 0.33e2/0.32e2 -0.43e2/0.144e3;
4f5a65f49035 Fixed the upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents: 47
diff changeset
17 0.11e2/0.144e3 -0.13e2/0.32e2 -0.29e2/0.48e2 0.389e3/0.288e3;
4f5a65f49035 Fixed the upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents: 47
diff changeset
18 0.1e1/0.32e2 -0.11e2/0.144e3 -0.65e2/0.288e3 -0.13e2/0.16e2;
4f5a65f49035 Fixed the upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents: 47
diff changeset
19 ];
47
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 Qp(1:4,1:4)=Q_U;
52
4f5a65f49035 Fixed the upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents: 47
diff changeset
22 Qp(m-3:m,m-3:m)=rot90(Q_U,2)'; %%% This is different from standard SBP
47
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 Qm=-Qp';
ebd25e5a481a Added upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
25
52
4f5a65f49035 Fixed the upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents: 47
diff changeset
26 e_1=sparse(m,1);e_1(1)=1;
4f5a65f49035 Fixed the upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents: 47
diff changeset
27 e_m=sparse(m,1);e_m(m)=1;
47
ebd25e5a481a Added upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
28
ebd25e5a481a Added upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
29 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
30
ebd25e5a481a Added upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
31 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
32 end