Mercurial > repos > public > sbplib

comparison +sbp/+implementations/d2_noneq_variable_6.m @ 1344:b4e5e45bd239 feature/D2_boundary_opt

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Remove round off zeros from D2Nonequidistant operators
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Sat, 15 Oct 2022 15:48:20 +0200
parents bcdb14b05d03
children
comparison
equal deleted inserted replaced
1343:09a5783a3d37 1344:b4e5e45bd239
136 %%%% Difference operators %%% 136 %%%% Difference operators %%%
137 D1 = H \ Q; 137 D1 = H \ Q;
138 138
139 % Helper functions for constructing D2(c) 139 % Helper functions for constructing D2(c)
140 % TODO: Consider changing sparse(diag(...)) to spdiags(....) 140 % TODO: Consider changing sparse(diag(...)) to spdiags(....)
141 min_inds = sbp.implementations.d2_sparsity_pattern_inds(m, order, BP, 0, 3);
142 nonmin_inds = sbp.implementations.d2_sparsity_pattern_inds(m, order, BP, 1, 3);
141 143
142 % Minimal 7 point stencil width 144 % Minimal 7 point stencil width
143 function D2 = D2_fun_minimal(c) 145 function D2 = D2_fun_minimal(c)
144 % Here we add variable diffusion 146 % Here we add variable diffusion
145 C1 = sparse(diag(c)); 147 C1 = sparse(diag(c));
150 C3 = sparse(diag(C3 * c)); 152 C3 = sparse(diag(C3 * c));
151 153
152 % Remainder term added to wide second derivative operator 154 % Remainder term added to wide second derivative operator
153 R = (1/3600 / h) * transpose(DD_6) * C1 * DD_6 + (1/600 / h) * transpose(DD_5) * C2 * DD_5 + (1/80 / h) * transpose(DD_4) * C3 * DD_4; 155 R = (1/3600 / h) * transpose(DD_6) * C1 * DD_6 + (1/600 / h) * transpose(DD_5) * C2 * DD_5 + (1/80 / h) * transpose(DD_4) * C3 * DD_4;
154 D2 = D1 * C1 * D1 - H \ R; 156 D2 = D1 * C1 * D1 - H \ R;
157
158 % Remove potential round off zeros
159 D2_tmp = sparse(m,m);
160 D2_tmp(min_inds) = D2(min_inds);
161 D2 = D2_tmp;
155 end 162 end
156 163
157 % Non-minimal 9 point stencil width 164 % Non-minimal 9 point stencil width
158 function D2 = D2_fun_nonminimal(c) 165 function D2 = D2_fun_nonminimal(c)
159 % Here we add variable diffusion 166 % Here we add variable diffusion
163 C2 = sparse(diag(C2 * c)); 170 C2 = sparse(diag(C2 * c));
164 171
165 % Remainder term added to wide second derivative operator 172 % Remainder term added to wide second derivative operator
166 R = (1/3600 / h) * transpose(DD_6) * C1 * DD_6 + (1/600 / h) * transpose(DD_5) * C2 * DD_5; 173 R = (1/3600 / h) * transpose(DD_6) * C1 * DD_6 + (1/600 / h) * transpose(DD_5) * C2 * DD_5;
167 D2 = D1 * C1 * D1 - H \ R; 174 D2 = D1 * C1 * D1 - H \ R;
175
176 % Remove potential round off zeros
177 D2_tmp = sparse(m,m);
178 D2_tmp(nonmin_inds) = D2(nonmin_inds);
179 D2 = D2_tmp;
168 end 180 end
169 181
170 % Wide stencil 182 % Wide stencil
171 function D2 = D2_fun_wide(c) 183 function D2 = D2_fun_wide(c)
172 % Here we add variable diffusion 184 % Here we add variable diffusion