Mercurial > repos > public > sbplib

comparison +sbp/+implementations/d2_noneq_variable_8.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
182 DD_8(m - 7:m, m - 11:m) = [0.99112262457261614959e0 -0.79189382269199820415e1 0.27669555809755867922e2 -0.55208801421796203457e2 0.68768103448343019521e2 -0.54693640344683071895e2 0.27031526461705782491e2 -0.74650882590612358780e1 0.82615990808320718845e0 0 0 0; 0 0.97411941507587258409e0 -0.77621755876022479980e1 0.27021942543523272783e2 -0.53627357074178360373e2 0.66225803906075245320e2 -0.51828716216050995413e2 0.24540872831466741457e2 -0.63048462739199410204e1 0.76035645561041265933e0 0 0; 0 0 0.97783801558437253538e0 -0.77894362989766952447e1 0.27093191644704986344e2 -0.53644769362408908249e2 0.65804221723591524304e2 -0.50115808052716920383e2 0.23888038475126051744e2 -0.76427065234692260122e1 0.14294303785648149613e1 0; 0 0 0 0.10642345748642342173e1 -0.85747647355583382784e1 0.30254875259437643639e2 -0.60939919211778894557e2 0.75491869612115314813e2 -0.64530162797168947524e2 0.45789440151310044599e2 -0.29767555907566203748e2 0.11211983054345146839e2; 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0; ]; 182 DD_8(m - 7:m, m - 11:m) = [0.99112262457261614959e0 -0.79189382269199820415e1 0.27669555809755867922e2 -0.55208801421796203457e2 0.68768103448343019521e2 -0.54693640344683071895e2 0.27031526461705782491e2 -0.74650882590612358780e1 0.82615990808320718845e0 0 0 0; 0 0.97411941507587258409e0 -0.77621755876022479980e1 0.27021942543523272783e2 -0.53627357074178360373e2 0.66225803906075245320e2 -0.51828716216050995413e2 0.24540872831466741457e2 -0.63048462739199410204e1 0.76035645561041265933e0 0 0; 0 0 0.97783801558437253538e0 -0.77894362989766952447e1 0.27093191644704986344e2 -0.53644769362408908249e2 0.65804221723591524304e2 -0.50115808052716920383e2 0.23888038475126051744e2 -0.76427065234692260122e1 0.14294303785648149613e1 0; 0 0 0 0.10642345748642342173e1 -0.85747647355583382784e1 0.30254875259437643639e2 -0.60939919211778894557e2 0.75491869612115314813e2 -0.64530162797168947524e2 0.45789440151310044599e2 -0.29767555907566203748e2 0.11211983054345146839e2; 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0; ];
183 DD_8 = sparse(DD_8); 183 DD_8 = sparse(DD_8);
184 184
185 %%%%%%%%%%%%%%%%%%%%%%%%%%%%% 185 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
186 186
187
187 %%%% Difference operators %% 188 %%%% Difference operators %%
188 D1 = H \ Q; 189 D1 = H \ Q;
189 190
190 % Helper functions for constructing D2(c) 191 % Helper functions for constructing D2(c)
191 % TODO: Consider changing sparse(diag(...)) to spdiags(....) 192 % TODO: Consider changing sparse(diag(...)) to spdiags(....)
192 193 min_inds = sbp.implementations.d2_sparsity_pattern_inds(m,order,BP,0,4);
194 nonmin_inds = sbp.implementations.d2_sparsity_pattern_inds(m,order,BP,1,4);
193 % Minimal 9 point stencil width 195 % Minimal 9 point stencil width
194 function D2 = D2_fun_minimal(c) 196 function D2 = D2_fun_minimal(c)
195 % Here we add variable diffusion 197 % Here we add variable diffusion
196 C1 = sparse(diag(c)); 198 C1 = sparse(diag(c));
197 C2 = 1/2 * diag(ones(m - 1, 1), -1) + 1/2 * diag(ones(m, 1), 0); C2(1, 2) = 1/2; 199 C2 = 1/2 * diag(ones(m - 1, 1), -1) + 1/2 * diag(ones(m, 1), 0); C2(1, 2) = 1/2;
204 206
205 % Remainder term added to wide second derivative operator 207 % Remainder term added to wide second derivative operator
206 R = (1/78400 / h) * transpose(DD_8) * C1 * DD_8 + (1/14700 / h) * transpose(DD_7) * C2 * DD_7 + (1/2520 / h) * transpose(DD_6) * C3 * DD_6 + (1/350 / h) * transpose(DD_5) * C4 * DD_5; 208 R = (1/78400 / h) * transpose(DD_8) * C1 * DD_8 + (1/14700 / h) * transpose(DD_7) * C2 * DD_7 + (1/2520 / h) * transpose(DD_6) * C3 * DD_6 + (1/350 / h) * transpose(DD_5) * C4 * DD_5;
207 209
208 D2 = D1 * C1 * D1 - H \ R; 210 D2 = D1 * C1 * D1 - H \ R;
211
212 % Remove potential round off zeros
213 D2_tmp = sparse(m,m);
214 D2_tmp(min_inds) = D2(min_inds);
215 D2 = D2_tmp;
209 end 216 end
210 217
211 % Non-minimal 11 point stencil width 218 % Non-minimal 11 point stencil width
212 function D2 = D2_fun_nonminimal(c) 219 function D2 = D2_fun_nonminimal(c)
213 % Here we add variable diffusion 220 % Here we add variable diffusion
219 C3 = sparse(diag(C3 * c)); 226 C3 = sparse(diag(C3 * c));
220 227
221 % Remainder term added to wide second derivative operator 228 % Remainder term added to wide second derivative operator
222 R = (1/78400 / h) * transpose(DD_8) * C1 * DD_8 + (1/14700 / h) * transpose(DD_7) * C2 * DD_7 + (1/2520 / h) * transpose(DD_6) * C3 * DD_6; 229 R = (1/78400 / h) * transpose(DD_8) * C1 * DD_8 + (1/14700 / h) * transpose(DD_7) * C2 * DD_7 + (1/2520 / h) * transpose(DD_6) * C3 * DD_6;
223 D2 = D1 * C1 * D1 - H \ R; 230 D2 = D1 * C1 * D1 - H \ R;
231
232 % Remove potential round off zeros
233 D2_tmp = sparse(m,m);
234 D2_tmp(nonmin_inds) = D2(nonmin_inds);
235 D2 = D2_tmp;
224 end 236 end
225 237
226 % Wide stencil 238 % Wide stencil
227 function D2 = D2_fun_wide(c) 239 function D2 = D2_fun_wide(c)
228 % Here we add variable diffusion 240 % Here we add variable diffusion