Mercurial > repos > public > sbplib

annotate +sbp/upwind3_3bp.m @ 87:0a29a60e0b21

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
In Curve: Rearranged for speed. arc_length_fun is now a property of Curve. If it is not supplied, it is computed via the derivative and spline fitting. Switching to the arc length parameterization is much faster now. The new stuff can be tested with testArcLength.m (which should be deleted after that).
author Martin Almquist <martin.almquist@it.uu.se>
date Sun, 29 Nov 2015 22:23:09 +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] = upwind3(m,h)
52
4f5a65f49035 Fixed the upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents: 47
diff changeset
2 Hv = ones(m,1);
4f5a65f49035 Fixed the upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents: 47
diff changeset
3 Hv(1:3) = [3/8; 7/6; 23/24];
4f5a65f49035 Fixed the upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents: 47
diff changeset
4 Hv(m-2:m) = rot90(Hv(1:3),2);
4f5a65f49035 Fixed the upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents: 47
diff changeset
5 Hv = Hv*h;
4f5a65f49035 Fixed the upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents: 47
diff changeset
6 H = spdiag(Hv,0);
4f5a65f49035 Fixed the upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents: 47
diff changeset
7 HI = spdiag(1./Hv,0);
47
ebd25e5a481a Added upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
8
52
4f5a65f49035 Fixed the upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents: 47
diff changeset
9 q_diags = [-1, 0, 1, 2];
4f5a65f49035 Fixed the upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents: 47
diff changeset
10 q_stencil = [-1/3 -1/2 1 -1/6];
4f5a65f49035 Fixed the upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents: 47
diff changeset
11 Qp = stripeMatrix(q_stencil, q_diags,m);
47
ebd25e5a481a Added upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
12
52
4f5a65f49035 Fixed the upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents: 47
diff changeset
13 Q_U = [
4f5a65f49035 Fixed the upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents: 47
diff changeset
14 -1/24 17/24 -1/6;
4f5a65f49035 Fixed the upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents: 47
diff changeset
15 -13/24 -1/4 23/24;
4f5a65f49035 Fixed the upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents: 47
diff changeset
16 1/12 -11/24 -11/24;
4f5a65f49035 Fixed the upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents: 47
diff changeset
17 ];
47
ebd25e5a481a Added upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
18 Qp(1:3,1:3)=Q_U;
52
4f5a65f49035 Fixed the upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents: 47
diff changeset
19 Qp(m-2:m,m-2:m)=rot90(Q_U,2)'; %%% This is different from standard SBP
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 Qm=-Qp';
ebd25e5a481a Added upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
22
52
4f5a65f49035 Fixed the upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents: 47
diff changeset
23 e_1=sparse(m,1);
4f5a65f49035 Fixed the upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents: 47
diff changeset
24 e_1(1)=1;
4f5a65f49035 Fixed the upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents: 47
diff changeset
25 e_m=sparse(m,1);
4f5a65f49035 Fixed the upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents: 47
diff changeset
26 e_m(m)=1;
47
ebd25e5a481a Added upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
27
ebd25e5a481a Added upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
28 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
29
ebd25e5a481a Added upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
30 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
31 end
ebd25e5a481a Added upwind operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
32