Mercurial > repos > public > sbplib

changeset 1337:bf2554f1825d feature/D2_boundary_opt

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add periodic D1 and D2 operators for orders 8,10,12
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Fri, 13 May 2022 13:28:10 +0200
parents 0666629aa183
children da61892884a4
files +sbp/+implementations/d2_variable_periodic_10.m +sbp/+implementations/d2_variable_periodic_12.m +sbp/+implementations/d2_variable_periodic_8.m +sbp/D2VariablePeriodic.m
diffstat 4 files changed, 265 insertions(+), 0 deletions(-) [+]
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+sbp/+implementations/d2_variable_periodic_10.m	Fri May 13 13:28:10 2022 +0200
@@ -0,0 +1,80 @@
+function [H, HI, D1, D2, e_l, e_r, d1_l, d1_r] = d2_variable_periodic_10(m,h)
+    % m = number of unique grid points, i.e. h = L/m;
+
+    if(m<11)
+        error(['Operator requires at least ' num2str(8) ' grid points']);
+    end
+
+    % Norm
+    Hv = ones(m,1);
+    Hv = h*Hv;
+    H = spdiag(Hv, 0);
+    HI = spdiag(1./Hv, 0);
+
+
+    % Dummy boundary operators
+    e_l = sparse(m,1);
+    e_r = rot90(e_l, 2);
+
+    d1_l = sparse(m,1);
+    d1_r = -rot90(d1_l, 2);
+
+
+    % D1 operator
+    diags   = -5:5;
+    stencil = [-1/1260, 5/504, -5/84, 5/21, -5/6, 0, 5/6, -5/21, 5/84, -5/504, 1/1260];
+    D1 = stripeMatrixPeriodic(stencil, diags, m);
+    D1 = D1/h;
+
+    % Undivided differences    
+    diags   = -3:3;
+    stencil = [1 -6 15 -20 15 -6 1];
+    DD_6 = stripeMatrixPeriodic(stencil, diags, m);    
+
+    diags   = -4:3;
+    stencil = [-1 7 -21 35 -35 21 -7 1];
+    DD_7 = stripeMatrixPeriodic(stencil, diags, m);
+
+    diags   = -4:4;
+    stencil = [1 -8 28 -56 70 -56 28 -8 1];
+    DD_8 = stripeMatrixPeriodic(stencil, diags, m);
+
+    diags   = -5:4;
+    stencil = [-1 9 -36 84 -126 126 -84 36 -9 1];
+    DD_9 = stripeMatrixPeriodic(stencil, diags, m);
+    
+    diags   = -5:5;
+    stencil = [1 -10 45 -120 210 -252 210 -120 45 -10 1];
+    DD_10 = stripeMatrixPeriodic(stencil, diags, m);
+    % D2 operator
+    function D2 = D2_fun(c)
+
+        diags   = -1:0;
+        stencil = [1/2, 1/2];
+        C2 = stripeMatrixPeriodic(stencil, diags, m);
+
+        diags   = -1:1;
+        stencil = [2/7, 3/7, 2/7];
+        C3 = stripeMatrixPeriodic(stencil, diags, m);
+
+        diags   = -2:1;
+        stencil = [1/5, 3/10, 1/5, 3/10];
+        C4 = stripeMatrixPeriodic(stencil, diags, m);
+
+        diags   = -2:2;
+        stencil = [1/5, 1/5, 1/5, 1/5, 1/5];
+        C5 = stripeMatrixPeriodic(stencil, diags, m);
+
+        C1 = sparse(diag(c));
+        C2 = sparse(diag(C2 * c));
+        C3 = sparse(diag(C3 * c));
+        C4 = sparse(diag(C4 * c));
+
+        % 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;
+    end
+    D2 = @D2_fun;
+
+
+end
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+sbp/+implementations/d2_variable_periodic_12.m	Fri May 13 13:28:10 2022 +0200
@@ -0,0 +1,89 @@
+function [H, HI, D1, D2, e_l, e_r, d1_l, d1_r] = d2_variable_periodic_12(m,h)
+    % m = number of unique grid points, i.e. h = L/m;
+
+    if(m<13)
+        error(['Operator requires at least ' num2str(8) ' grid points']);
+    end
+
+    % Norm
+    Hv = ones(m,1);
+    Hv = h*Hv;
+    H = spdiag(Hv, 0);
+    HI = spdiag(1./Hv, 0);
+
+
+    % Dummy boundary operators
+    e_l = sparse(m,1);
+    e_r = rot90(e_l, 2);
+
+    d1_l = sparse(m,1);
+    d1_r = -rot90(d1_l, 2);
+
+
+    % D1 operator
+    diags   = -6:6;
+    stencil = [1/5544, -1/385, 1/56, -5/63, 15/56, -6/7, 0, 6/7, -15/56, 5/63, -1/56, 1/385, -1/5544];
+    D1 = stripeMatrixPeriodic(stencil, diags, m);
+    D1 = D1/h;
+
+    % Undivided differences      
+    diags   = -4:3;
+    stencil = [-1 7 -21 35 -35 21 -7 1];
+    DD_7 = stripeMatrixPeriodic(stencil, diags, m);
+
+    diags   = -4:4;
+    stencil = [1 -8 28 -56 70 -56 28 -8 1];
+    DD_8 = stripeMatrixPeriodic(stencil, diags, m);
+
+    diags   = -5:4;
+    stencil = [-1 9 -36 84 -126 126 -84 36 -9 1];
+    DD_9 = stripeMatrixPeriodic(stencil, diags, m);
+    
+    diags   = -5:5;
+    stencil = [1 -10 45 -120 210 -252 210 -120 45 -10 1];
+    DD_10 = stripeMatrixPeriodic(stencil, diags, m);
+
+    diags   = -6:5;
+    stencil = [-1, 11, -55, 165, -330, 462, -462, 330, -165, 55, -11, 1];
+    DD_11 = stripeMatrixPeriodic(stencil, diags, m);
+
+    diags   = -6:6;
+    stencil = [1, -12, 66, -220, 495, -792, 924, -792, 495, -220, 66,-12, 1];
+    DD_12 = stripeMatrixPeriodic(stencil, diags, m);
+
+    % D2 operator
+    function D2 = D2_fun(c)
+
+        diags   = -1:0;
+        stencil = [1/2, 1/2];
+        C2 = stripeMatrixPeriodic(stencil, diags, m);
+
+        diags   = -1:1;
+        stencil = [5/18, 4/9, 5/18];
+        C3 = stripeMatrixPeriodic(stencil, diags, m);
+
+        diags   = -2:1;
+        stencil = [5/28, 9/28, 9/28, 5/28];
+        C4 = stripeMatrixPeriodic(stencil, diags, m);
+
+        diags   = -2:2;
+        stencil = [1/7, 8/35, 9/35, 8/35, 1/7];
+        C5 = stripeMatrixPeriodic(stencil, diags, m);
+
+        diags   = -3:2;
+        stencil = [1/6, 1/6, 1/6, 1/6, 1/6, 1/6];
+        C6 = stripeMatrixPeriodic(stencil, diags, m);
+
+        C1 = sparse(diag(c));
+        C2 = sparse(diag(C2 * c));
+        C3 = sparse(diag(C3 * c));
+        C4 = sparse(diag(C4 * c));
+
+        % 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;
+    end
+    D2 = @D2_fun;
+
+
+end
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+sbp/+implementations/d2_variable_periodic_8.m	Fri May 13 13:28:10 2022 +0200
@@ -0,0 +1,75 @@
+function [H, HI, D1, D2, e_l, e_r, d1_l, d1_r] = d2_variable_periodic_8(m,h)
+    % m = number of unique grid points, i.e. h = L/m;
+
+    if(m<9)
+        error(['Operator requires at least ' num2str(8) ' grid points']);
+    end
+
+    % Norm
+    Hv = ones(m,1);
+    Hv = h*Hv;
+    H = spdiag(Hv, 0);
+    HI = spdiag(1./Hv, 0);
+
+
+    % Dummy boundary operators
+    e_l = sparse(m,1);
+    e_r = rot90(e_l, 2);
+
+    d1_l = sparse(m,1);
+    d1_r = -rot90(d1_l, 2);
+
+
+    % D1 operator
+    diags   = -4:4;
+    stencil = [1/280, -4/105, 1/5, -4/5, 0, 4/5, -1/5, 4/105, -1/280];
+    D1 = stripeMatrixPeriodic(stencil, diags, m);
+    D1 = D1/h;
+
+    % Undivided differences
+    diags   = -3:2;
+    stencil = [-1 5 -10 10 -5 1];
+    DD_5 = stripeMatrixPeriodic(stencil, diags, m);
+
+    diags   = -3:3;
+    stencil = [1 -6 15 -20 15 -6 1];
+    DD_6 = stripeMatrixPeriodic(stencil, diags, m);
+
+    diags   = -4:3;
+    stencil = [-1 7 -21 35 -35 21 -7 1];
+    DD_7 = stripeMatrixPeriodic(stencil, diags, m);
+
+    diags   = -4:4;
+    stencil = [1 -8 28 -56 70 -56 28 -8 1];
+    DD_8 = stripeMatrixPeriodic(stencil, diags, m);
+
+    % D2 operator
+    function D2 = D2_fun(c)
+
+        C1 = sparse(diag(c));
+        
+        diags   = -1:0;
+        stencil = [1/2, 1/2];
+        C2 = stripeMatrixPeriodic(stencil, diags, m);
+
+        diags   = -1:1;
+        stencil = [3/10, 2/5, 3/10];
+        C3 = stripeMatrixPeriodic(stencil, diags, m);
+
+        diags   = -2:1;
+        stencil = [1/4, 1/4, 1/4, 1/4];
+        C4 = stripeMatrixPeriodic(stencil, diags, m);
+
+        C2 = sparse(diag(C2 * c));
+        C3 = sparse(diag(C3 * c));
+        C4 = sparse(diag(C4 * c));
+
+        % 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 + (1/350 / h) * transpose(DD_5) * C4 * DD_5;
+
+        D2 = D1 * C1 * D1 - H \ R;
+    end
+    D2 = @D2_fun;
+
+
+end
--- a/+sbp/D2VariablePeriodic.m	Fri May 13 13:26:16 2022 +0200
+++ b/+sbp/D2VariablePeriodic.m	Fri May 13 13:28:10 2022 +0200
@@ -28,6 +28,27 @@
 
             switch order
 
+            case 12
+                [obj.H, obj.HI, obj.D1, obj.D2, obj.e_l,...
+                    obj.e_r, obj.d1_l, obj.d1_r] = ...
+                    sbp.implementations.d2_variable_periodic_12(m,obj.h);
+                obj.borrowing.M.d1 = obj.H(1, 1) / obj.h; % First element in H/h is borrowing also for M
+                obj.borrowing.R.delta_D = inf;
+                
+            case 10
+                [obj.H, obj.HI, obj.D1, obj.D2, obj.e_l,...
+                    obj.e_r, obj.d1_l, obj.d1_r] = ...
+                    sbp.implementations.d2_variable_periodic_10(m,obj.h);
+                obj.borrowing.M.d1 = obj.H(1, 1) / obj.h; % First element in H/h is borrowing also for M
+                obj.borrowing.R.delta_D = inf;
+
+            case 8
+                [obj.H, obj.HI, obj.D1, obj.D2, obj.e_l,...
+                    obj.e_r, obj.d1_l, obj.d1_r] = ...
+                    sbp.implementations.d2_variable_periodic_8(m,obj.h);
+                obj.borrowing.M.d1 = obj.H(1, 1) / obj.h; % First element in H/h is borrowing also for M
+                obj.borrowing.R.delta_D = inf;
+
                 case 6
                     [obj.H, obj.HI, obj.D1, obj.D2, obj.e_l,...
                         obj.e_r, obj.d1_l, obj.d1_r] = ...