Mercurial > repos > public > sbplib

changeset 686:5ccf6aaf6d6b feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add D2VariablePeriodic orders 4 and 6.
author Martin Almquist <malmquist@stanford.edu>
date Thu, 08 Feb 2018 16:44:46 -0800
parents b035902869a8
children e8fc3aa1faf6
files +sbp/+implementations/d2_variable_periodic_2.m +sbp/+implementations/d2_variable_periodic_4.m +sbp/+implementations/d2_variable_periodic_6.m +sbp/D2VariablePeriodic.m
diffstat 4 files changed, 121 insertions(+), 11 deletions(-) [+]
line wrap: on
line diff
--- a/+sbp/+implementations/d2_variable_periodic_2.m	Thu Feb 08 16:43:43 2018 -0800
+++ b/+sbp/+implementations/d2_variable_periodic_2.m	Thu Feb 08 16:44:46 2018 -0800
@@ -1,9 +1,8 @@
 function [H, HI, D1, D2, e_l, e_r, d1_l, d1_r] = d2_variable_periodic_2(m,h)
     % m = number of unique grid points, i.e. h = L/m;
 
-    BP = 1;
-    if(m<2*BP)
-        error(['Operator requires at least ' num2str(2*BP) ' grid points']);
+    if(m<3)
+        error(['Operator requires at least ' num2str(3) ' grid points']);
     end
 
     % Norm
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+sbp/+implementations/d2_variable_periodic_4.m	Thu Feb 08 16:44:46 2018 -0800
@@ -0,0 +1,57 @@
+function [H, HI, D1, D2, e_l, e_r, d1_l, d1_r] = d2_variable_periodic_4(m,h)
+    % m = number of unique grid points, i.e. h = L/m;
+
+    if(m<5)
+        error(['Operator requires at least ' num2str(5) ' 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);
+
+    S = d1_l*d1_l' + d1_r*d1_r';
+
+    % D1 operator
+    stencil = [1/12 -2/3 0 2/3 -1/12];
+    diags = -2:2;
+    Q = stripeMatrixPeriodic(stencil, diags, m);
+    D1 = HI*(Q - 1/2*e_l*e_l' + 1/2*e_r*e_r');
+
+
+    scheme_width = 5;
+    scheme_radius = (scheme_width-1)/2;
+    
+    r = 1:m;
+    offset = scheme_width;
+    r = r + offset;
+
+    function D2 = D2_fun(c)
+        c = [c(end-scheme_width+1:end); c; c(1:scheme_width) ];
+
+        % Note: these coefficients are for -M.
+        Mm2 = -1/8*c(r-2) + 1/6*c(r-1) - 1/8*c(r);
+        Mm1 = 1/6 *c(r-2) + 1/2*c(r-1) + 1/2*c(r) + 1/6*c(r+1);
+        M0  = -1/24*c(r-2)- 5/6*c(r-1) - 3/4*c(r) - 5/6*c(r+1) - 1/24*c(r+2);
+        Mp1  = 0 * c(r-2) + 1/6*c(r-1) + 1/2*c(r) + 1/2*c(r+1) + 1/6 *c(r+2);
+        Mp2  = 0 * c(r-2) + 0 * c(r-1) - 1/8*c(r) + 1/6*c(r+1) - 1/8 *c(r+2);
+
+        vals = -[Mm2,Mm1,M0,Mp1,Mp2];
+        diags = -scheme_radius : scheme_radius;
+        M = spdiagsVariablePeriodic(vals,diags); 
+
+        M=M/h;
+        D2=HI*(-M );
+
+    end
+    D2 = @D2_fun;
+end
\ No newline at end of file
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+sbp/+implementations/d2_variable_periodic_6.m	Thu Feb 08 16:44:46 2018 -0800
@@ -0,0 +1,58 @@
+function [H, HI, D1, D2, e_l, e_r, d1_l, d1_r] = d2_variable_periodic_6(m,h)
+    % m = number of unique grid points, i.e. h = L/m;
+
+    if(m<7)
+        error(['Operator requires at least ' num2str(7) ' 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   = -3:3;
+    stencil = [-1/60 9/60 -45/60 0 45/60 -9/60 1/60];
+    D1 = stripeMatrixPeriodic(stencil, diags, m);
+    D1 = D1/h;
+
+    % D2 operator
+    scheme_width = 7;
+    scheme_radius = (scheme_width-1)/2;
+
+    r = 1:m;
+    offset = scheme_width;
+    r = r + offset;
+
+    function D2 = D2_fun(c)
+        c = [c(end-scheme_width+1:end); c; c(1:scheme_width) ];
+
+        Mm3 =  c(r-2)/0.40e2 + c(r-1)/0.40e2 - 0.11e2/0.360e3 * c(r-3) - 0.11e2/0.360e3 * c(r);
+        Mm2 =  c(r-3)/0.20e2 - 0.3e1/0.10e2 * c(r-1) + c(r+1)/0.20e2 + 0.7e1/0.40e2 * c(r) + 0.7e1/0.40e2 * c(r-2);
+        Mm1 = -c(r-3)/0.40e2 - 0.3e1/0.10e2 * c(r-2) - 0.3e1/0.10e2 * c(r+1) - c(r+2)/0.40e2 - 0.17e2/0.40e2 * c(r) - 0.17e2/0.40e2 * c(r-1);
+        M0 =  c(r-3)/0.180e3 + c(r-2)/0.8e1 + 0.19e2/0.20e2 * c(r-1) + 0.19e2/0.20e2 * c(r+1) + c(r+2)/0.8e1 + c(r+3)/0.180e3 + 0.101e3/0.180e3 * c(r);
+        Mp1 = -c(r-2)/0.40e2 - 0.3e1/0.10e2 * c(r-1) - 0.3e1/0.10e2 * c(r+2) - c(r+3)/0.40e2 - 0.17e2/0.40e2 * c(r) - 0.17e2/0.40e2 * c(r+1);
+        Mp2 =  c(r-1)/0.20e2 - 0.3e1/0.10e2 * c(r+1) + c(r+3)/0.20e2 + 0.7e1/0.40e2 * c(r) + 0.7e1/0.40e2 * c(r+2);
+        Mp3 =  c(r+1)/0.40e2 + c(r+2)/0.40e2 - 0.11e2/0.360e3 * c(r) - 0.11e2/0.360e3 * c(r+3);
+
+        vals = [Mm3,Mm2,Mm1,M0,Mp1,Mp2,Mp3];
+        diags = -scheme_radius : scheme_radius;
+        M = spdiagsVariablePeriodic(vals,diags); 
+
+        M=M/h;
+        D2=HI*(-M );
+    end
+    D2 = @D2_fun;
+
+    
+end
--- a/+sbp/D2VariablePeriodic.m	Thu Feb 08 16:43:43 2018 -0800
+++ b/+sbp/D2VariablePeriodic.m	Thu Feb 08 16:44:46 2018 -0800
@@ -29,19 +29,14 @@
             switch order
 
                 case 6
-                     error('Not impl')
-
-                    [obj.H, obj.HI, obj.D1, obj.D2, ...
-                    ~, obj.e_l, obj.e_r, ~, ~, ~, ~, ~,...
-                     obj.d1_l, obj.d1_r] = ...
-                        sbp.implementations.d4_variable_periodic_6(m, obj.h);
+                    [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_6(m,obj.h);
                     obj.borrowing.M.d1 = 0.1878;
                     obj.borrowing.R.delta_D = 0.3696;
                     % Borrowing e^T*D1 - d1 from R
 
                 case 4
-                    error('Not impl')
-
                     [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_4(m,obj.h);
@@ -63,6 +58,7 @@
                     error('Invalid operator order %d.',order);
             end
             obj.borrowing.H11 = obj.H(1,1)/obj.h; % First element in H/h,
+
             obj.m = m;
             obj.M = [];
         end