Mercurial > repos > public > sbplib

comparison +sbp/+implementations/d2_variable_periodic_4.m @ 686:5ccf6aaf6d6b feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add D2VariablePeriodic orders 4 and 6.
author Martin Almquist <malmquist@stanford.edu>
date Thu, 08 Feb 2018 16:44:46 -0800
parents
children bbf303c1f0cf
comparison
equal deleted inserted replaced
685:b035902869a8 686:5ccf6aaf6d6b
1 function [H, HI, D1, D2, e_l, e_r, d1_l, d1_r] = d2_variable_periodic_4(m,h)
2 % m = number of unique grid points, i.e. h = L/m;
3
4 if(m<5)
5 error(['Operator requires at least ' num2str(5) ' grid points']);
6 end
7
8 % Norm
9 Hv = ones(m,1);
10 Hv = h*Hv;
11 H = spdiag(Hv, 0);
12 HI = spdiag(1./Hv, 0);
13
14
15 % Dummy boundary operators
16 e_l = sparse(m,1);
17 e_r = rot90(e_l, 2);
18
19 d1_l = sparse(m,1);
20 d1_r = -rot90(d1_l, 2);
21
22 S = d1_l*d1_l' + d1_r*d1_r';
23
24 % D1 operator
25 stencil = [1/12 -2/3 0 2/3 -1/12];
26 diags = -2:2;
27 Q = stripeMatrixPeriodic(stencil, diags, m);
28 D1 = HI*(Q - 1/2*e_l*e_l' + 1/2*e_r*e_r');
29
30
31 scheme_width = 5;
32 scheme_radius = (scheme_width-1)/2;
33
34 r = 1:m;
35 offset = scheme_width;
36 r = r + offset;
37
38 function D2 = D2_fun(c)
39 c = [c(end-scheme_width+1:end); c; c(1:scheme_width) ];
40
41 % Note: these coefficients are for -M.
42 Mm2 = -1/8*c(r-2) + 1/6*c(r-1) - 1/8*c(r);
43 Mm1 = 1/6 *c(r-2) + 1/2*c(r-1) + 1/2*c(r) + 1/6*c(r+1);
44 M0 = -1/24*c(r-2)- 5/6*c(r-1) - 3/4*c(r) - 5/6*c(r+1) - 1/24*c(r+2);
45 Mp1 = 0 * c(r-2) + 1/6*c(r-1) + 1/2*c(r) + 1/2*c(r+1) + 1/6 *c(r+2);
46 Mp2 = 0 * c(r-2) + 0 * c(r-1) - 1/8*c(r) + 1/6*c(r+1) - 1/8 *c(r+2);
47
48 vals = -[Mm2,Mm1,M0,Mp1,Mp2];
49 diags = -scheme_radius : scheme_radius;
50 M = spdiagsVariablePeriodic(vals,diags);
51
52 M=M/h;
53 D2=HI*(-M );
54
55 end
56 D2 = @D2_fun;
57 end