Mercurial > repos > public > sbplib

comparison +sbp/+implementations/d2_noneq_variable_4.m @ 1325:1b0f2415237f feature/D2_boundary_opt

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add variable coefficient boundary-optimized second derivatives.
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Sun, 13 Feb 2022 19:32:34 +0100
parents
children c2d716c4f1ed
comparison
equal deleted inserted replaced
1301:8978521b0f06 1325:1b0f2415237f
1 function [H, HI, D1, D2, DI] = d2_noneq_variable_4(N, h, options)
2 % N: Number of grid points
3 % h: grid spacing
4 % options: struct containing options for constructing the operator
5 % current options are:
6 % options.stencil_type ('minimal','nonminimal','wide')
7 % options.AD ('upwind', 'op')
8
9 % BP: Number of boundary points
10 % order: Accuracy of interior stencil
11 BP = 4;
12 order = 4;
13
14 %%%% Norm matrix %%%%%%%%
15 P = zeros(BP, 1);
16 P0 = 2.1259737557798e-01;
17 P1 = 1.0260290400758e+00;
18 P2 = 1.0775123588954e+00;
19 P3 = 9.8607273802835e-01;
20
21 for i = 0:BP - 1
22 P(i + 1) = eval(['P' num2str(i)]);
23 end
24
25 Hv = ones(N, 1);
26 Hv(1:BP) = P;
27 Hv(end - BP + 1:end) = flip(P);
28 Hv = h * Hv;
29 H = spdiags(Hv, 0, N, N);
30 HI = spdiags(1 ./ Hv, 0, N, N);
31 %%%%%%%%%%%%%%%%%%%%%%%%%
32
33 %%%% Q matrix %%%%%%%%%%%
34 d = [1/12, -2/3, 0, 2/3, -1/12];
35 d = repmat(d, N, 1);
36 Q = spdiags(d, -order / 2:order / 2, N, N);
37
38 % Boundaries
39 Q0_0 = -5.0000000000000e-01;
40 Q0_1 = 6.5605279837843e-01;
41 Q0_2 = -1.9875859409017e-01;
42 Q0_3 = 4.2705795711740e-02;
43 Q0_4 = 0.0000000000000e+00;
44 Q0_5 = 0.0000000000000e+00;
45 Q1_0 = -6.5605279837843e-01;
46 Q1_1 = 0.0000000000000e+00;
47 Q1_2 = 8.1236966439895e-01;
48 Q1_3 = -1.5631686602052e-01;
49 Q1_4 = 0.0000000000000e+00;
50 Q1_5 = 0.0000000000000e+00;
51 Q2_0 = 1.9875859409017e-01;
52 Q2_1 = -8.1236966439895e-01;
53 Q2_2 = 0.0000000000000e+00;
54 Q2_3 = 6.9694440364211e-01;
55 Q2_4 = -8.3333333333333e-02;
56 Q2_5 = 0.0000000000000e+00;
57 Q3_0 = -4.2705795711740e-02;
58 Q3_1 = 1.5631686602052e-01;
59 Q3_2 = -6.9694440364211e-01;
60 Q3_3 = 0.0000000000000e+00;
61 Q3_4 = 6.6666666666667e-01;
62 Q3_5 = -8.3333333333333e-02;
63
64 for i = 1:BP
65
66 for j = 1:BP
67 Q(i, j) = eval(['Q' num2str(i - 1) '_' num2str(j - 1)]);
68 Q(N + 1 - i, N + 1 - j) = -eval(['Q' num2str(i - 1) '_' num2str(j - 1)]);
69 end
70
71 end
72
73 %%%%%%%%%%%%%%%%%%%%%%%%%%%
74
75 %%% Undivided difference operators %%%%
76 % Closed with zeros at the first boundary nodes.
77 m = N;
78
79 DD_2 = (diag(ones(m - 1, 1), -1) - 2 * diag(ones(m, 1), 0) + diag(ones(m - 1, 1), 1));
80 DD_2(1:3, 1:4) = [0 0 0 0; 0.16138369498429727170e1 -0.26095138364100825853e1 0.99567688656710986834e0 0; 0 0.84859980956172494512e0 -0.17944203477786665350e1 0.94582053821694158989e0; ];
81 DD_2(m - 2:m, m - 3:m) = [0.94582053821694158989e0 -0.17944203477786665350e1 0.84859980956172494512e0 0; 0 0.99567688656710986834e0 -0.26095138364100825853e1 0.16138369498429727170e1; 0 0 0 0; ];
82 DD_2 = sparse(DD_2);
83
84 DD_3 = (-diag(ones(m - 2, 1), -2) + 3 * diag(ones(m - 1, 1), -1) - 3 * diag(ones(m, 1), 0) + diag(ones(m - 1, 1), 1));
85 DD_3(1:4, 1:5) = [0 0 0 0 0; 0 0 0 0 0; -0.17277463987989539852e1 0.37021976718569105700e1 -0.29870306597013296050e1 0.10125793866433730203e1 0; 0 -0.81738495424057284493e0 0.26916305216679998025e1 -0.28374616146508247697e1 0.96321604722339781208e0; ];
86 DD_3(m - 2:m, m - 4:m) = [-0.96321604722339781208e0 0.28374616146508247697e1 -0.26916305216679998025e1 0.81738495424057284493e0 0; 0 -0.10125793866433730203e1 0.29870306597013296050e1 -0.37021976718569105700e1 0.17277463987989539852e1; 0 0 0 0 0; ];
87 DD_3 = sparse(DD_3);
88
89 DD_4 = (diag(ones(m - 2, 1), 2) - 4 * diag(ones(m - 1, 1), 1) + 6 * diag(ones(m, 1), 0) - 4 * diag(ones(m - 1, 1), -1) + diag(ones(m - 2, 1), -2));
90 DD_4(1:4, 1:6) = [0 0 0 0 0 0; 0 0 0 0 0 0; 0.18176226052481525189e1 -0.47546882767009058782e1 0.59740613194026592100e1 -0.40503175465734920811e1 0.10133218986235862303e1 0; 0 0.79462567299107735362e0 -0.35888406955573330700e1 0.56749232293016495393e1 -0.38528641888935912483e1 0.97215598215819742539e0; ];
91 DD_4(m - 3:m, m - 5:m) = [0.97215598215819742539e0 -0.38528641888935912483e1 0.56749232293016495393e1 -0.35888406955573330700e1 0.79462567299107735362e0 0; 0 0.10133218986235862303e1 -0.40503175465734920811e1 0.59740613194026592100e1 -0.47546882767009058782e1 0.18176226052481525189e1; 0 0 0 0 0 0; 0 0 0 0 0 0; ];
92 DD_4 = sparse(DD_4);
93 %%%%%%%%%%%%%%%%%%%%%%%%%%
94
95 %%%% Difference operators %%%
96 D1 = H \ Q;
97
98 % Helper functions for constructing D2(c)
99 % TODO: Consider changing sparse(diag(...)) to spdiags(....)
100
101 % Minimal 5 point stencil width
102 function D2 = D2_fun_minimal(c)
103 % Here we add variable diffusion
104 C1 = sparse(diag(c));
105 C2 = 1/2 * diag(ones(m - 1, 1), -1) + 1/2 * diag(ones(m, 1), 0); C2(1, 2) = 1/2;
106
107 C2 = sparse(diag(C2 * c));
108
109 % Remainder term added to wide second drivative opereator, to obtain a 5
110 % point narrow stencil.
111 R = (1/144 / h) * transpose(DD_4) * C1 * DD_4 + (1/18 / h) * transpose(DD_3) * C2 * DD_3;
112 D2 = D1 * C1 * D1 - H \ R;
113 end
114
115 % Few additional grid point in interior stencil cmp. to minimal
116 function D2 = D2_fun_nonminimal(c)
117 % Here we add variable diffusion
118 C1 = sparse(diag(c));
119
120 % Remainder term added to wide second derivative operator
121 R = (1/144 / h) * transpose(DD_4) * C1 * DD_4;
122 D2 = D1 * C1 * D1 - H \ R;
123 end
124
125 % Wide stencil
126 function D2 = D2_wide(c)
127 % Here we add variable diffusion
128 C1 = sparse(diag(c));
129 D2 = D1 * C1 * D1;
130 end
131
132 switch options.stencil_width
133 case 'minimal'
134 D2 = @D2_fun_minimal;
135 case 'nonminimal'
136 D2 = @D2_fun_nonminimal;
137 case 'wide'
138 D2 = @D2_fun_wide;
139 otherwise
140 error('No option %s for stencil width', options.stencil_width)
141 end
142
143 %%%%%%%%%%%%%%%%%%%%%%%%%%%
144
145 %%%% Artificial dissipation operator %%%
146 switch options.AD
147 case 'upwind'
148 % This is the choice that yield 3rd order Upwind
149 DI = H \ (transpose(DD_2) * DD_2) * (-1/12);
150 case 'op'
151 % This choice will preserve the order of the underlying
152 % Non-dissipative D1 SBP operator
153 DI = H \ (transpose(DD_3) * DD_3) * (-1 / (5 * 12));
154 otherwise
155 error("Artificial dissipation options '%s' not implemented.", option.AD)
156 end
157
158 %%%%%%%%%%%%%%%%%%%%%%%%%%%
159 end