Mercurial > repos > public > sbplib

annotate +sbp/+implementations/d2_variable_periodic_6.m @ 1337:bf2554f1825d feature/D2_boundary_opt

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add periodic D1 and D2 operators for orders 8,10,12
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Fri, 13 May 2022 13:28:10 +0200
parents bbf303c1f0cf
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
686
5ccf6aaf6d6b Add D2VariablePeriodic orders 4 and 6.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
1 function [H, HI, D1, D2, e_l, e_r, d1_l, d1_r] = d2_variable_periodic_6(m,h)
5ccf6aaf6d6b Add D2VariablePeriodic orders 4 and 6.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
2 % m = number of unique grid points, i.e. h = L/m;
5ccf6aaf6d6b Add D2VariablePeriodic orders 4 and 6.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
3
5ccf6aaf6d6b Add D2VariablePeriodic orders 4 and 6.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
4 if(m<7)
5ccf6aaf6d6b Add D2VariablePeriodic orders 4 and 6.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
5 error(['Operator requires at least ' num2str(7) ' grid points']);
5ccf6aaf6d6b Add D2VariablePeriodic orders 4 and 6.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
6 end
5ccf6aaf6d6b Add D2VariablePeriodic orders 4 and 6.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
7
5ccf6aaf6d6b Add D2VariablePeriodic orders 4 and 6.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
8 % Norm
5ccf6aaf6d6b Add D2VariablePeriodic orders 4 and 6.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
9 Hv = ones(m,1);
5ccf6aaf6d6b Add D2VariablePeriodic orders 4 and 6.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
10 Hv = h*Hv;
5ccf6aaf6d6b Add D2VariablePeriodic orders 4 and 6.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
11 H = spdiag(Hv, 0);
5ccf6aaf6d6b Add D2VariablePeriodic orders 4 and 6.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
12 HI = spdiag(1./Hv, 0);
5ccf6aaf6d6b Add D2VariablePeriodic orders 4 and 6.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
13
5ccf6aaf6d6b Add D2VariablePeriodic orders 4 and 6.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
14
5ccf6aaf6d6b Add D2VariablePeriodic orders 4 and 6.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
15 % Dummy boundary operators
5ccf6aaf6d6b Add D2VariablePeriodic orders 4 and 6.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
16 e_l = sparse(m,1);
5ccf6aaf6d6b Add D2VariablePeriodic orders 4 and 6.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
17 e_r = rot90(e_l, 2);
5ccf6aaf6d6b Add D2VariablePeriodic orders 4 and 6.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
18
5ccf6aaf6d6b Add D2VariablePeriodic orders 4 and 6.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
19 d1_l = sparse(m,1);
5ccf6aaf6d6b Add D2VariablePeriodic orders 4 and 6.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
20 d1_r = -rot90(d1_l, 2);
5ccf6aaf6d6b Add D2VariablePeriodic orders 4 and 6.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
21
5ccf6aaf6d6b Add D2VariablePeriodic orders 4 and 6.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
22
5ccf6aaf6d6b Add D2VariablePeriodic orders 4 and 6.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
23 % D1 operator
5ccf6aaf6d6b Add D2VariablePeriodic orders 4 and 6.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
24 diags = -3:3;
5ccf6aaf6d6b Add D2VariablePeriodic orders 4 and 6.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
25 stencil = [-1/60 9/60 -45/60 0 45/60 -9/60 1/60];
5ccf6aaf6d6b Add D2VariablePeriodic orders 4 and 6.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
26 D1 = stripeMatrixPeriodic(stencil, diags, m);
5ccf6aaf6d6b Add D2VariablePeriodic orders 4 and 6.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
27 D1 = D1/h;
5ccf6aaf6d6b Add D2VariablePeriodic orders 4 and 6.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
28
5ccf6aaf6d6b Add D2VariablePeriodic orders 4 and 6.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
29 % D2 operator
5ccf6aaf6d6b Add D2VariablePeriodic orders 4 and 6.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
30 scheme_width = 7;
5ccf6aaf6d6b Add D2VariablePeriodic orders 4 and 6.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
31 scheme_radius = (scheme_width-1)/2;
5ccf6aaf6d6b Add D2VariablePeriodic orders 4 and 6.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
32
5ccf6aaf6d6b Add D2VariablePeriodic orders 4 and 6.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
33 r = 1:m;
5ccf6aaf6d6b Add D2VariablePeriodic orders 4 and 6.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
34 offset = scheme_width;
5ccf6aaf6d6b Add D2VariablePeriodic orders 4 and 6.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
35 r = r + offset;
5ccf6aaf6d6b Add D2VariablePeriodic orders 4 and 6.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
36
5ccf6aaf6d6b Add D2VariablePeriodic orders 4 and 6.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
37 function D2 = D2_fun(c)
5ccf6aaf6d6b Add D2VariablePeriodic orders 4 and 6.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
38 c = [c(end-scheme_width+1:end); c; c(1:scheme_width) ];
5ccf6aaf6d6b Add D2VariablePeriodic orders 4 and 6.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
39
5ccf6aaf6d6b Add D2VariablePeriodic orders 4 and 6.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
40 Mm3 = c(r-2)/0.40e2 + c(r-1)/0.40e2 - 0.11e2/0.360e3 * c(r-3) - 0.11e2/0.360e3 * c(r);
5ccf6aaf6d6b Add D2VariablePeriodic orders 4 and 6.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
41 Mm2 = c(r-3)/0.20e2 - 0.3e1/0.10e2 * c(r-1) + c(r+1)/0.20e2 + 0.7e1/0.40e2 * c(r) + 0.7e1/0.40e2 * c(r-2);
5ccf6aaf6d6b Add D2VariablePeriodic orders 4 and 6.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
42 Mm1 = -c(r-3)/0.40e2 - 0.3e1/0.10e2 * c(r-2) - 0.3e1/0.10e2 * c(r+1) - c(r+2)/0.40e2 - 0.17e2/0.40e2 * c(r) - 0.17e2/0.40e2 * c(r-1);
5ccf6aaf6d6b Add D2VariablePeriodic orders 4 and 6.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
43 M0 = c(r-3)/0.180e3 + c(r-2)/0.8e1 + 0.19e2/0.20e2 * c(r-1) + 0.19e2/0.20e2 * c(r+1) + c(r+2)/0.8e1 + c(r+3)/0.180e3 + 0.101e3/0.180e3 * c(r);
5ccf6aaf6d6b Add D2VariablePeriodic orders 4 and 6.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
44 Mp1 = -c(r-2)/0.40e2 - 0.3e1/0.10e2 * c(r-1) - 0.3e1/0.10e2 * c(r+2) - c(r+3)/0.40e2 - 0.17e2/0.40e2 * c(r) - 0.17e2/0.40e2 * c(r+1);
5ccf6aaf6d6b Add D2VariablePeriodic orders 4 and 6.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
45 Mp2 = c(r-1)/0.20e2 - 0.3e1/0.10e2 * c(r+1) + c(r+3)/0.20e2 + 0.7e1/0.40e2 * c(r) + 0.7e1/0.40e2 * c(r+2);
5ccf6aaf6d6b Add D2VariablePeriodic orders 4 and 6.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
46 Mp3 = c(r+1)/0.40e2 + c(r+2)/0.40e2 - 0.11e2/0.360e3 * c(r) - 0.11e2/0.360e3 * c(r+3);
5ccf6aaf6d6b Add D2VariablePeriodic orders 4 and 6.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
47
5ccf6aaf6d6b Add D2VariablePeriodic orders 4 and 6.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
48 vals = [Mm3,Mm2,Mm1,M0,Mp1,Mp2,Mp3];
5ccf6aaf6d6b Add D2VariablePeriodic orders 4 and 6.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
49 diags = -scheme_radius : scheme_radius;
801
bbf303c1f0cf Rename spdaigsVariablePeriodic spdiagsPeriodic
Jonatan Werpers <jonatan@werpers.com>
parents: 686
diff changeset
50 M = spdiagsPeriodic(vals,diags);
686
5ccf6aaf6d6b Add D2VariablePeriodic orders 4 and 6.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
51
5ccf6aaf6d6b Add D2VariablePeriodic orders 4 and 6.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
52 M=M/h;
5ccf6aaf6d6b Add D2VariablePeriodic orders 4 and 6.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
53 D2=HI*(-M );
5ccf6aaf6d6b Add D2VariablePeriodic orders 4 and 6.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
54 end
5ccf6aaf6d6b Add D2VariablePeriodic orders 4 and 6.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
55 D2 = @D2_fun;
5ccf6aaf6d6b Add D2VariablePeriodic orders 4 and 6.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
56
801
bbf303c1f0cf Rename spdaigsVariablePeriodic spdiagsPeriodic
Jonatan Werpers <jonatan@werpers.com>
parents: 686
diff changeset
57
686
5ccf6aaf6d6b Add D2VariablePeriodic orders 4 and 6.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
58 end