Mercurial > repos > public > sbplib

annotate +sbp/upwind4.m @ 105:1df4f3704b76 feature/arclen-param

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Added error checking for all methods that use the derivative gp. Removed that from the to do list. Updated some comments.
author Martin Almquist <martin.almquist@it.uu.se>
date Mon, 07 Dec 2015 18:51:41 +0100
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