Mercurial > repos > public > sbplib

comparison +sbp/+implementations/d2_noneq_variable_10.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 8e9df030a0a5
children
comparison
equal deleted inserted replaced
1343:09a5783a3d37 1344:b4e5e45bd239
247 %%%% Difference operator %% 247 %%%% Difference operator %%
248 D1 = H \ Q; 248 D1 = H \ Q;
249 249
250 % Helper functions for constructing D2(c) 250 % Helper functions for constructing D2(c)
251 % TODO: Consider changing sparse(diag(...)) to spdiags(....) 251 % TODO: Consider changing sparse(diag(...)) to spdiags(....)
252 min_inds = sbp.implementations.d2_sparsity_pattern_inds(m, order, BP, 0, 5);
253 nonmin_inds = sbp.implementations.d2_sparsity_pattern_inds(m, order, BP, 1, 5);
252 254
253 % Minimal 11 point stencil width 255 % Minimal 11 point stencil width
254 function D2 = D2_fun_minimal(c) 256 function D2 = D2_fun_minimal(c)
255 % Here we add variable diffusion 257 % Here we add variable diffusion
256 C1 = sparse(diag(c)); 258 C1 = sparse(diag(c));
265 C5 = sparse(diag(C5 * c)); 267 C5 = sparse(diag(C5 * c));
266 268
267 % Remainder term added to wide second derivative operator 269 % Remainder term added to wide second derivative operator
268 R = (1/1587600 / h) * transpose(DD_10) * C1 * DD_10 + (1/317520 / h) * transpose(DD_9) * C2 * DD_9 + (1/60480 / h) * transpose(DD_8) * C3 * DD_8 + (1/10584 / h) * transpose(DD_7) * C4 * DD_7 + (1/1512 / h) * transpose(DD_6) * C5 * DD_6; 270 R = (1/1587600 / h) * transpose(DD_10) * C1 * DD_10 + (1/317520 / h) * transpose(DD_9) * C2 * DD_9 + (1/60480 / h) * transpose(DD_8) * C3 * DD_8 + (1/10584 / h) * transpose(DD_7) * C4 * DD_7 + (1/1512 / h) * transpose(DD_6) * C5 * DD_6;
269 D2 = D1 * C1 * D1 - H \ R; 271 D2 = D1 * C1 * D1 - H \ R;
272
273 % Remove potential round off zeros
274 D2_tmp = sparse(m,m);
275 D2_tmp(min_inds) = D2(min_inds);
276 D2 = D2_tmp;
270 end 277 end
271 278
272 % Non-minimal 13 point stencil width 279 % Non-minimal 13 point stencil width
273 function D2 = D2_fun_nonminimal(c) 280 function D2 = D2_fun_nonminimal(c)
274 C1 = sparse(diag(c)); 281 C1 = sparse(diag(c));
281 C4 = sparse(diag(C4 * c)); 288 C4 = sparse(diag(C4 * c));
282 289
283 % Remainder term added to wide second derivative operator 290 % Remainder term added to wide second derivative operator
284 R = (1/1587600 / h) * transpose(DD_10) * C1 * DD_10 + (1/317520 / h) * transpose(DD_9) * C2 * DD_9 + (1/60480 / h) * transpose(DD_8) * C3 * DD_8 + (1/10584 / h) * transpose(DD_7) * C4 * DD_7; 291 R = (1/1587600 / h) * transpose(DD_10) * C1 * DD_10 + (1/317520 / h) * transpose(DD_9) * C2 * DD_9 + (1/60480 / h) * transpose(DD_8) * C3 * DD_8 + (1/10584 / h) * transpose(DD_7) * C4 * DD_7;
285 D2 = D1 * C1 * D1 - H \ R; 292 D2 = D1 * C1 * D1 - H \ R;
293
294 % Remove potential round off zeros
295 D2_tmp = sparse(m,m);
296 D2_tmp(nonmin_inds) = D2(nonmin_inds);
297 D2 = D2_tmp;
286 end 298 end
287 299
288 % Wide stencil 300 % Wide stencil
289 function D2 = D2_fun_wide(c) 301 function D2 = D2_fun_wide(c)
290 % Here we add variable diffusion 302 % Here we add variable diffusion