Mercurial > repos > public > sbplib

comparison +sbp/+implementations/d2_noneq_variable_12.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
322 %%%% Difference operators %% 322 %%%% Difference operators %%
323 D1 = H \ Q; 323 D1 = H \ Q;
324 324
325 % Helper functions for constructing D2(c) 325 % Helper functions for constructing D2(c)
326 % TODO: Consider changing sparse(diag(...)) to spdiags(....) 326 % TODO: Consider changing sparse(diag(...)) to spdiags(....)
327 min_inds = sbp.implementations.d2_sparsity_pattern_inds(m, order, BP, 0, 6);
328 nonmin_inds = sbp.implementations.d2_sparsity_pattern_inds(m, order, BP, 1, 6);
327 329
328 % Minimal 13 point stencil width 330 % Minimal 13 point stencil width
329 function D2 = D2_fun_minimal(c) 331 function D2 = D2_fun_minimal(c)
330 % Here we add variable diffusion 332 % Here we add variable diffusion
331 C1 = sparse(diag(c)); 333 C1 = sparse(diag(c));
343 C6 = sparse(diag(C6 * c)); 345 C6 = sparse(diag(C6 * c));
344 346
345 % Remainder term added to wide second derivative operator 347 % Remainder term added to wide second derivative operator
346 R = (1/30735936 / h) * transpose(DD_12) * C1 * DD_12 + (1/6403320 / h) * transpose(DD_11) * C2 * DD_11 + (1/1293600 / h) * transpose(DD_10) * C3 * DD_10 + (1/249480 / h) * transpose(DD_9) * C4 * DD_9 + (1/44352 / h) * transpose(DD_8) * C5 * DD_8 + (1/6468 / h) * transpose(DD_7) * C6 * DD_7; 348 R = (1/30735936 / h) * transpose(DD_12) * C1 * DD_12 + (1/6403320 / h) * transpose(DD_11) * C2 * DD_11 + (1/1293600 / h) * transpose(DD_10) * C3 * DD_10 + (1/249480 / h) * transpose(DD_9) * C4 * DD_9 + (1/44352 / h) * transpose(DD_8) * C5 * DD_8 + (1/6468 / h) * transpose(DD_7) * C6 * DD_7;
347 D2 = D1 * C1 * D1 - H \ R; 349 D2 = D1 * C1 * D1 - H \ R;
350
351 % Remove potential round off zeros
352 D2_tmp = sparse(m,m);
353 D2_tmp(min_inds) = D2(min_inds);
354 D2 = D2_tmp;
348 end 355 end
349 356
350 % Non-minimal 15 point stencil width 357 % Non-minimal 15 point stencil width
351 function D2 = D2_fun_nonminimal(c) 358 function D2 = D2_fun_nonminimal(c)
352 % Here we add variable diffusion 359 % Here we add variable diffusion
362 C5 = sparse(diag(C5 * c)); 369 C5 = sparse(diag(C5 * c));
363 370
364 % Remainder term added to wide second derivative operator 371 % Remainder term added to wide second derivative operator
365 R = (1/30735936 / h) * transpose(DD_12) * C1 * DD_12 + (1/6403320 / h) * transpose(DD_11) * C2 * DD_11 + (1/1293600 / h) * transpose(DD_10) * C3 * DD_10 + (1/249480 / h) * transpose(DD_9) * C4 * DD_9 + (1/44352 / h) * transpose(DD_8) * C5 * DD_8; 372 R = (1/30735936 / h) * transpose(DD_12) * C1 * DD_12 + (1/6403320 / h) * transpose(DD_11) * C2 * DD_11 + (1/1293600 / h) * transpose(DD_10) * C3 * DD_10 + (1/249480 / h) * transpose(DD_9) * C4 * DD_9 + (1/44352 / h) * transpose(DD_8) * C5 * DD_8;
366 D2 = D1 * C1 * D1 - H \ R; 373 D2 = D1 * C1 * D1 - H \ R;
374
375 % Remove potential round off zeros
376 D2_tmp = sparse(m,m);
377 D2_tmp(nonmin_inds) = D2(nonmin_inds);
378 D2 = D2_tmp;
367 end 379 end
368 380
369 % Wide stencil 381 % Wide stencil
370 function D2 = D2_fun_wide(c) 382 function D2 = D2_fun_wide(c)
371 % Here we add variable diffusion 383 % Here we add variable diffusion