Mercurial > repos > public > sbplib

changeset 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 09a5783a3d37
children
files +sbp/+implementations/d2_noneq_variable_10.m +sbp/+implementations/d2_noneq_variable_12.m +sbp/+implementations/d2_noneq_variable_4.m +sbp/+implementations/d2_noneq_variable_6.m +sbp/+implementations/d2_noneq_variable_8.m +sbp/+implementations/d2_sparsity_pattern_inds.m
diffstat 6 files changed, 74 insertions(+), 1 deletions(-) [+]
line wrap: on
line diff
--- a/+sbp/+implementations/d2_noneq_variable_10.m	Fri Oct 14 14:42:32 2022 +0200
+++ b/+sbp/+implementations/d2_noneq_variable_10.m	Sat Oct 15 15:48:20 2022 +0200
@@ -249,6 +249,8 @@
 
     % Helper functions for constructing D2(c)
     % TODO: Consider changing sparse(diag(...)) to spdiags(....)
+    min_inds = sbp.implementations.d2_sparsity_pattern_inds(m, order, BP, 0, 5);
+    nonmin_inds = sbp.implementations.d2_sparsity_pattern_inds(m, order, BP, 1, 5);
 
     % Minimal 11 point stencil width
     function D2 = D2_fun_minimal(c)
@@ -267,6 +269,11 @@
         % Remainder term added to wide second derivative operator
         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;
         D2 = D1 * C1 * D1 - H \ R;
+
+        % Remove potential round off zeros
+        D2_tmp = sparse(m,m);
+        D2_tmp(min_inds) = D2(min_inds);
+        D2 = D2_tmp;
     end
 
     %  Non-minimal 13 point stencil width
@@ -283,6 +290,11 @@
         % Remainder term added to wide second derivative operator
         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;
         D2 = D1 * C1 * D1 - H \ R;
+
+        % Remove potential round off zeros
+        D2_tmp = sparse(m,m);
+        D2_tmp(nonmin_inds) = D2(nonmin_inds);
+        D2 = D2_tmp;
     end
 
     % Wide stencil
--- a/+sbp/+implementations/d2_noneq_variable_12.m	Fri Oct 14 14:42:32 2022 +0200
+++ b/+sbp/+implementations/d2_noneq_variable_12.m	Sat Oct 15 15:48:20 2022 +0200
@@ -324,6 +324,8 @@
 
     % Helper functions for constructing D2(c)
     % TODO: Consider changing sparse(diag(...)) to spdiags(....)
+    min_inds = sbp.implementations.d2_sparsity_pattern_inds(m, order, BP, 0, 6);
+    nonmin_inds = sbp.implementations.d2_sparsity_pattern_inds(m, order, BP, 1, 6);
 
     % Minimal 13 point stencil width
     function D2 = D2_fun_minimal(c)
@@ -345,6 +347,11 @@
         % Remainder term added to wide second derivative operator
         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;
         D2 = D1 * C1 * D1 - H \ R;
+
+        % Remove potential round off zeros
+        D2_tmp = sparse(m,m);
+        D2_tmp(min_inds) = D2(min_inds);
+        D2 = D2_tmp;
     end
 
     %  Non-minimal 15 point stencil width
@@ -364,6 +371,11 @@
         % Remainder term added to wide second derivative operator
         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;
         D2 = D1 * C1 * D1 - H \ R;
+
+        % Remove potential round off zeros
+        D2_tmp = sparse(m,m);
+        D2_tmp(nonmin_inds) = D2(nonmin_inds);
+        D2 = D2_tmp;
     end
 
     % Wide stencil
--- a/+sbp/+implementations/d2_noneq_variable_4.m	Fri Oct 14 14:42:32 2022 +0200
+++ b/+sbp/+implementations/d2_noneq_variable_4.m	Sat Oct 15 15:48:20 2022 +0200
@@ -100,6 +100,8 @@
 
     % Helper functions for constructing D2(c)
     % TODO: Consider changing sparse(diag(...)) to spdiags(....)
+    min_inds = sbp.implementations.d2_sparsity_pattern_inds(m, order, BP, 0, 2);
+    nonmin_inds = sbp.implementations.d2_sparsity_pattern_inds(m, order, BP, 1, 2);
 
     % Minimal 5 point stencil width
     function D2 = D2_fun_minimal(c)
@@ -113,6 +115,11 @@
         % point narrow stencil.
         R = (1/144 / h) * transpose(DD_4) * C1 * DD_4 + (1/18 / h) * transpose(DD_3) * C2 * DD_3;
         D2 = D1 * C1 * D1 - H \ R;
+
+        % Remove potential round off zeros
+        D2tmp = sparse(m,m);
+        D2tmp(min_inds) = D2(min_inds);
+        D2 = D2tmp;
     end
 
     %  Non-minimal 7 point stencil width
@@ -123,6 +130,11 @@
         % Remainder term added to wide second derivative operator
         R = (1/144 / h) * transpose(DD_4) * C1 * DD_4;
         D2 = D1 * C1 * D1 - H \ R;
+        
+        % Remove potential round off zeros
+        D2tmp = sparse(m,m);
+        D2tmp(nonmin_inds) = D2(nonmin_inds);
+        D2 = D2tmp;
     end
 
     % Wide stencil
--- a/+sbp/+implementations/d2_noneq_variable_6.m	Fri Oct 14 14:42:32 2022 +0200
+++ b/+sbp/+implementations/d2_noneq_variable_6.m	Sat Oct 15 15:48:20 2022 +0200
@@ -138,6 +138,8 @@
 
     % Helper functions for constructing D2(c)
     % TODO: Consider changing sparse(diag(...)) to spdiags(....)
+    min_inds = sbp.implementations.d2_sparsity_pattern_inds(m, order, BP, 0, 3);
+    nonmin_inds = sbp.implementations.d2_sparsity_pattern_inds(m, order, BP, 1, 3);
 
     % Minimal 7 point stencil width
     function D2 = D2_fun_minimal(c)
@@ -152,6 +154,11 @@
         % Remainder term added to wide second derivative operator
         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;
         D2 = D1 * C1 * D1 - H \ R;
+
+        % Remove potential round off zeros
+        D2_tmp = sparse(m,m);
+        D2_tmp(min_inds) = D2(min_inds);
+        D2 = D2_tmp;
     end
 
     %  Non-minimal 9 point stencil width
@@ -165,6 +172,11 @@
         % Remainder term added to wide second derivative operator
         R = (1/3600 / h) * transpose(DD_6) * C1 * DD_6 + (1/600 / h) * transpose(DD_5) * C2 * DD_5;
         D2 = D1 * C1 * D1 - H \ R;
+
+        % Remove potential round off zeros
+        D2_tmp = sparse(m,m);
+        D2_tmp(nonmin_inds) = D2(nonmin_inds);
+        D2 = D2_tmp;
     end
 
     % Wide stencil
--- a/+sbp/+implementations/d2_noneq_variable_8.m	Fri Oct 14 14:42:32 2022 +0200
+++ b/+sbp/+implementations/d2_noneq_variable_8.m	Sat Oct 15 15:48:20 2022 +0200
@@ -184,12 +184,14 @@
 
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
+
     %%%% Difference operators %%
     D1 = H \ Q;
 
     % Helper functions for constructing D2(c)
     % TODO: Consider changing sparse(diag(...)) to spdiags(....)
-
+    min_inds =  sbp.implementations.d2_sparsity_pattern_inds(m,order,BP,0,4);
+    nonmin_inds = sbp.implementations.d2_sparsity_pattern_inds(m,order,BP,1,4);
     % Minimal 9 point stencil width
     function D2 = D2_fun_minimal(c)
         % Here we add variable diffusion
@@ -206,6 +208,11 @@
         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;
 
         D2 = D1 * C1 * D1 - H \ R;
+
+        % Remove potential round off zeros
+        D2_tmp = sparse(m,m);
+        D2_tmp(min_inds) = D2(min_inds);
+        D2 = D2_tmp;
     end
 
     %  Non-minimal 11 point stencil width
@@ -221,6 +228,11 @@
         % Remainder term added to wide second derivative operator
         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;
         D2 = D1 * C1 * D1 - H \ R;
+        
+        % Remove potential round off zeros
+        D2_tmp = sparse(m,m);
+        D2_tmp(nonmin_inds) = D2(nonmin_inds);
+        D2 = D2_tmp;
     end
 
     % Wide stencil
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+sbp/+implementations/d2_sparsity_pattern_inds.m	Sat Oct 15 15:48:20 2022 +0200
@@ -0,0 +1,13 @@
+function inds = d2_sparsity_pattern_inds(m, order, BP, interior_offset, closure_offset)
+    sparsity_pattern = sparse(m,m);
+    inner_stencil_inds = -order/2-interior_offset:order/2+interior_offset;
+    inner_stencil_pattern = ones(m,length(inner_stencil_inds));
+    sparsity_pattern = spdiags(inner_stencil_pattern,inner_stencil_inds,sparsity_pattern);
+    sparsity_pattern(1:BP,1:BP+closure_offset) = 1;
+    sparsity_pattern(end-BP+1:end,end-(BP+closure_offset)+1:end) = 1;
+    for k = 1:closure_offset
+        sparsity_pattern(BP+k,1:BP+closure_offset+k) = 1;
+        sparsity_pattern(end-(BP+k)+1,end-(BP+k+closure_offset)+1:end) = 1;
+    end
+    inds  = find(sparsity_pattern);
+end
\ No newline at end of file