Mercurial > repos > public > sbplib

annotate +sbp/+implementations/d2_variable_periodic_2.m @ 681:7368affc8f78 feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add D2 variable periodic for second order.
author Martin Almquist <malmquist@stanford.edu>
date Wed, 07 Feb 2018 15:42:50 -0800
parents
children 5ccf6aaf6d6b
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
681
7368affc8f78 Add D2 variable periodic for second order.
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_2(m,h)
7368affc8f78 Add D2 variable periodic for second order.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
2 % m = number of unique grid points, i.e. h = L/m;
7368affc8f78 Add D2 variable periodic for second order.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
3
7368affc8f78 Add D2 variable periodic for second order.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
4 BP = 1;
7368affc8f78 Add D2 variable periodic for second order.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
5 if(m<2*BP)
7368affc8f78 Add D2 variable periodic for second order.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
6 error(['Operator requires at least ' num2str(2*BP) ' grid points']);
7368affc8f78 Add D2 variable periodic for second order.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
7 end
7368affc8f78 Add D2 variable periodic for second order.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
8
7368affc8f78 Add D2 variable periodic for second order.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
9 % Norm
7368affc8f78 Add D2 variable periodic for second order.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
10 Hv = ones(m,1);
7368affc8f78 Add D2 variable periodic for second order.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
11 Hv = h*Hv;
7368affc8f78 Add D2 variable periodic for second order.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
12 H = spdiag(Hv, 0);
7368affc8f78 Add D2 variable periodic for second order.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
13 HI = spdiag(1./Hv, 0);
7368affc8f78 Add D2 variable periodic for second order.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
14
7368affc8f78 Add D2 variable periodic for second order.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
15
7368affc8f78 Add D2 variable periodic for second order.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
16 % Dummy boundary operators
7368affc8f78 Add D2 variable periodic for second order.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
17 e_l = sparse(m,1);
7368affc8f78 Add D2 variable periodic for second order.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
18 e_r = rot90(e_l, 2);
7368affc8f78 Add D2 variable periodic for second order.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
19
7368affc8f78 Add D2 variable periodic for second order.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
20 d1_l = sparse(m,1);
7368affc8f78 Add D2 variable periodic for second order.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
21 d1_r = -rot90(d1_l, 2);
7368affc8f78 Add D2 variable periodic for second order.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
22
7368affc8f78 Add D2 variable periodic for second order.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
23 % D1 operator
7368affc8f78 Add D2 variable periodic for second order.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
24 diags = -1:1;
7368affc8f78 Add D2 variable periodic for second order.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
25 stencil = [-1/2 0 1/2];
7368affc8f78 Add D2 variable periodic for second order.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
26 D1 = stripeMatrixPeriodic(stencil, diags, m);
7368affc8f78 Add D2 variable periodic for second order.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
27 D1 = D1/h;
7368affc8f78 Add D2 variable periodic for second order.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
28
7368affc8f78 Add D2 variable periodic for second order.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
29 scheme_width = 3;
7368affc8f78 Add D2 variable periodic for second order.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
30 scheme_radius = (scheme_width-1)/2;
7368affc8f78 Add D2 variable periodic for second order.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
31
7368affc8f78 Add D2 variable periodic for second order.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
32 r = 1:m;
7368affc8f78 Add D2 variable periodic for second order.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
33 offset = scheme_width;
7368affc8f78 Add D2 variable periodic for second order.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
34 r = r + offset;
7368affc8f78 Add D2 variable periodic for second order.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
35
7368affc8f78 Add D2 variable periodic for second order.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
36 function D2 = D2_fun(c)
7368affc8f78 Add D2 variable periodic for second order.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
37 c = [c(end-scheme_width+1:end); c; c(1:scheme_width) ];
7368affc8f78 Add D2 variable periodic for second order.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
38
7368affc8f78 Add D2 variable periodic for second order.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
39 Mm1 = -c(r-1)/2 - c(r)/2;
7368affc8f78 Add D2 variable periodic for second order.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
40 M0 = c(r-1)/2 + c(r) + c(r+1)/2;
7368affc8f78 Add D2 variable periodic for second order.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
41 Mp1 = -c(r)/2 - c(r+1)/2;
7368affc8f78 Add D2 variable periodic for second order.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
42
7368affc8f78 Add D2 variable periodic for second order.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
43 vals = [Mm1,M0,Mp1];
7368affc8f78 Add D2 variable periodic for second order.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
44 diags = -scheme_radius : scheme_radius;
7368affc8f78 Add D2 variable periodic for second order.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
45 M = spdiagsVariablePeriodic(vals,diags);
7368affc8f78 Add D2 variable periodic for second order.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
46
7368affc8f78 Add D2 variable periodic for second order.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
47 M=M/h;
7368affc8f78 Add D2 variable periodic for second order.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
48 D2=HI*(-M );
7368affc8f78 Add D2 variable periodic for second order.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
49 end
7368affc8f78 Add D2 variable periodic for second order.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
50 D2 = @D2_fun;
7368affc8f78 Add D2 variable periodic for second order.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
51 end