Mercurial > repos > public > sbplib

comparison +sbp/+implementations/d2_noneq_variable_4.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
98 %%%% Difference operators %%% 98 %%%% Difference operators %%%
99 D1 = H \ Q; 99 D1 = H \ Q;
100 100
101 % Helper functions for constructing D2(c) 101 % Helper functions for constructing D2(c)
102 % TODO: Consider changing sparse(diag(...)) to spdiags(....) 102 % TODO: Consider changing sparse(diag(...)) to spdiags(....)
103 min_inds = sbp.implementations.d2_sparsity_pattern_inds(m, order, BP, 0, 2);
104 nonmin_inds = sbp.implementations.d2_sparsity_pattern_inds(m, order, BP, 1, 2);
103 105
104 % Minimal 5 point stencil width 106 % Minimal 5 point stencil width
105 function D2 = D2_fun_minimal(c) 107 function D2 = D2_fun_minimal(c)
106 % Here we add variable diffusion 108 % Here we add variable diffusion
107 C1 = sparse(diag(c)); 109 C1 = sparse(diag(c));
111 113
112 % Remainder term added to wide second drivative opereator, to obtain a 5 114 % Remainder term added to wide second drivative opereator, to obtain a 5
113 % point narrow stencil. 115 % point narrow stencil.
114 R = (1/144 / h) * transpose(DD_4) * C1 * DD_4 + (1/18 / h) * transpose(DD_3) * C2 * DD_3; 116 R = (1/144 / h) * transpose(DD_4) * C1 * DD_4 + (1/18 / h) * transpose(DD_3) * C2 * DD_3;
115 D2 = D1 * C1 * D1 - H \ R; 117 D2 = D1 * C1 * D1 - H \ R;
118
119 % Remove potential round off zeros
120 D2tmp = sparse(m,m);
121 D2tmp(min_inds) = D2(min_inds);
122 D2 = D2tmp;
116 end 123 end
117 124
118 % Non-minimal 7 point stencil width 125 % Non-minimal 7 point stencil width
119 function D2 = D2_fun_nonminimal(c) 126 function D2 = D2_fun_nonminimal(c)
120 % Here we add variable diffusion 127 % Here we add variable diffusion
121 C1 = sparse(diag(c)); 128 C1 = sparse(diag(c));
122 129
123 % Remainder term added to wide second derivative operator 130 % Remainder term added to wide second derivative operator
124 R = (1/144 / h) * transpose(DD_4) * C1 * DD_4; 131 R = (1/144 / h) * transpose(DD_4) * C1 * DD_4;
125 D2 = D1 * C1 * D1 - H \ R; 132 D2 = D1 * C1 * D1 - H \ R;
133
134 % Remove potential round off zeros
135 D2tmp = sparse(m,m);
136 D2tmp(nonmin_inds) = D2(nonmin_inds);
137 D2 = D2tmp;
126 end 138 end
127 139
128 % Wide stencil 140 % Wide stencil
129 function D2 = D2_fun_wide(c) 141 function D2 = D2_fun_wide(c)
130 % Here we add variable diffusion 142 % Here we add variable diffusion