Mercurial > repos > public > sbplib

changeset 636:476b7e411ce3 feature/d1_staggered

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add staggered operators
author Martin Almquist <malmquist@stanford.edu>
date Tue, 07 Nov 2017 20:44:13 -0800
parents bb1180becc30
children eead18a8964d
files +sbp/+implementations/d1_staggered_2.m +sbp/+implementations/d1_staggered_4.m +sbp/+implementations/d1_staggered_6.m +sbp/D1Staggered.m
diffstat 4 files changed, 315 insertions(+), 0 deletions(-) [+]
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+sbp/+implementations/d1_staggered_2.m	Tue Nov 07 20:44:13 2017 -0800
@@ -0,0 +1,71 @@
+function [xp, xm, Pp, Pm, Qp, Qm] = d1_staggered_2(n, L)
+% [xp, xm, Pp, Pm, Qp, Qm] = sbp_staggered_2(n, L)
+% Input arguments:
+% n : number of grid points: n on xp grid and n+1 on xm grid
+% L : interval length
+
+n = n - 1;
+assert(n >= 3, 'Not enough grid points');  
+
+h = L/(n-1);
+
+qm20 = -1/4;
+qm01 = 1/2;
+qp11 = -1/4;
+qm11 = 1/4;
+qp02 = 1/4;
+pp0 = 1/2;
+pm1 = 1/4;
+qm00 = -1/2;
+qp10 = -1/2;
+qp12 = 3/4;
+qp00 = -1/2;
+qm21 = -3/4;
+pp1 = 1;
+qp01 = 1/4;
+pm0 = 1/2;
+pm2 = 5/4;
+qm10 = -1/4;
+
+% Q+ and Q-, top-left corner
+QpL = [...
+qp00, qp01, qp02;
+ qp10, qp11, qp12
+];
+QmL = [...
+qm00, qm01;
+ qm10, qm11;
+ qm20, qm21
+];
+
+% Q+ and Q-
+b = 2;
+w = b; 
+s = [-1, 1]';  
+Qp = spdiags(repmat(-s(end:-1:1)',[n+2 1]), -(w/2-1):w/2, n+2, n+2); 
+Qm = spdiags(repmat(s(:)',[n+2 1]), -(w/2-1)-1:w/2-1, n+2, n+2);
+Qp(end,:) = [];
+Qm(:,end) = [];
+
+% Add SBP boundary closures
+Qp(1:b,1:b+1) = QpL;
+Qp(end-b+1:end,end-b:end) = -fliplr(flipud(QpL));
+Qm(1:b+1,1:b) = QmL;
+Qm(end-b:end,end-b+1:end) = -fliplr(flipud(QmL));
+
+% P+ and P-
+Pp = ones(n+1,1);
+Pm = ones(n+2,1);
+
+Pp(1:b) = [pp0,  pp1]; 
+Pp(end-b+1:end) = Pp(b:-1:1);
+Pm(1:b+1) = [pm0,  pm1,  pm2];
+Pm(end-b:end) = Pm(b+1:-1:1);
+Pp = spdiags(Pp,0,n+1,n+1);
+Pm = spdiags(Pm,0,n+2,n+2);
+
+Pp = h*Pp;
+Pm = h*Pm;
+
+xp = h*[0:n]';
+xm = h*[0 1/2+0:n n]';  
\ No newline at end of file
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+sbp/+implementations/d1_staggered_4.m	Tue Nov 07 20:44:13 2017 -0800
@@ -0,0 +1,58 @@
+function [xp, xm, Pp, Pm, Qp, Qm] = d1_staggered_4(n, L)
+% [xp, xm, Pp, Pm, Qp, Qm] = sbp_staggered_4(n, L)
+% Input arguments:
+% n : number of grid points: n on xp grid and n+1 on xm grid
+% L : interval length
+
+n = n - 1;
+assert(n >= 7, 'Not enough grid points');  
+
+h = L/(n-1);
+
+% Q+ and Q-, top-left corner
+QpL = [...
+-323/378, 2783/3072, -1085/27648, -17/9216, -667/64512;
+ -5/21, -847/1024, 1085/1024, -37/1024, 899/21504;
+ 5/42, -121/1024, -1085/1024, 3575/3072, -753/7168;
+ -5/189, 121/3072, 1085/27648, -10759/9216, 74635/64512
+];
+QmL = [...
+-55/378, 5/21, -5/42, 5/189;
+ -2783/3072, 847/1024, 121/1024, -121/3072;
+ 1085/27648, -1085/1024, 1085/1024, -1085/27648;
+ 17/9216, 37/1024, -3575/3072, 10759/9216;
+ 667/64512, -899/21504, 753/7168, -74635/64512
+];
+
+% Q+ and Q-
+w = 4; 
+s = [ 1/24,  -9/8,  9/8,  -1/24];  
+Qp = spdiags(repmat(-s(end:-1:1),[n+2 1]), -(w/2-1):w/2, n+2, n+2); 
+Qm = spdiags(repmat(s(:)',[n+2 1]), -(w/2-1)-1:w/2-1, n+2, n+2);
+Qp(end,:) = [];
+Qm(:,end) = [];
+
+% Add SBP boundary closures
+bp = 4; 
+bm = 5;
+Qp(1:bp,1:bm) = double(QpL);
+Qp(end-bp+1:end,end-bm+1:end) = -fliplr(flipud(QpL));
+Qm(1:bm,1:bp) = double(QmL);
+Qm(end-bm+1:end,end-bp+1:end) = -fliplr(flipud(QmL));
+
+% P+ and P-
+Pp = ones(n+1,1);
+Pm = ones(n+2,1);
+
+Pp(1:bp) = [407/1152,  473/384,  343/384,  1177/1152]; 
+Pp(end-bp+1:end) = Pp(bp:-1:1);
+Pm(1:bm) = [5/63,  121/128,  1085/1152,  401/384,  2659/2688];
+Pm(end-bm+1:end) = Pm(bm:-1:1);
+Pp = spdiags(Pp,0,n+1,n+1);
+Pm = spdiags(Pm,0,n+2,n+2);
+
+Pp = h*Pp;
+Pm = h*Pm;
+
+xp = h*[0:n]';
+xm = h*[0 1/2+0:n n]';  
\ No newline at end of file
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+sbp/+implementations/d1_staggered_6.m	Tue Nov 07 20:44:13 2017 -0800
@@ -0,0 +1,103 @@
+function [xp, xm, Pp, Pm, Qp, Qm] = d1_staggered_6(n, L)
+% [xp, xm, Pp, Pm, Qp, Qm] = sbp_staggered_6(n, L)
+% Input arguments:
+% n : number of grid points: n on xp grid and n+1 on xm grid
+% L : interval length
+
+n = n - 1;
+assert(n >= 9, 'Not enough grid points');  
+
+h = L/(n-1);
+
+% Coefficients determined such that the SBP property is satisfied
+qp44 = -144662434693018337/131093333881896960;
+pp3 = 793/720;
+qp22 = -729260819980499/780317463582720;
+pm5 = 1262401/1244160;
+pp4 = 157/160;
+qp03 = -1056618403248889/17835827739033600;
+qp32 = 17518644727427/2340952390748160;
+pm4 = 224201/241920;
+pm2 = 90101/103680;
+qp40 = -572605772036807/27311111225395200;
+qp30 = 3572646778087103/61450000257139200;
+pm3 = 26369/23040;
+qp41 = 935926732952051/29964190601576448;
+qp02 = -67953322736279/2340952390748160;
+pp0 = 95/288;
+qp01 = 143647626126052207/149820953007882240;
+pp1 = 317/240;
+qp12 = 444281115722813/468190478149632;
+qp25 = 312883915256963/24970158834647040;
+qp11 = -36400588465707331/37455238251970560;
+qp20 = -137012145866261/13655555612697600;
+qp10 = -1008199398207353/6827777806348800;
+qp15 = -1088192686015471/37455238251970560;
+qp04 = -239061311631881/131093333881896960;
+qp35 = -2480203967217637/112365714755911680;
+qp21 = 1658590350384343/24970158834647040;
+qp05 = 4843362440877649/449462859023646720;
+qp24 = -372310079307409/4369777796063232;
+qp33 = -4809914536549211/4458956934758400;
+qp23 = 2826535075774271/2972637956505600;
+qp42 = 4202889834071/585238097687040;
+qp14 = 492703999899311/32773333470474240;
+pm1 = 133801/138240;
+qp34 = 36544940032227659/32773333470474240;
+qp00 = -216175739231516543/245800001028556800;
+qp43 = 128340370711031/17835827739033600;
+qp13 = 823082791653949/4458956934758400;
+pm0 = 5251/68040;
+pp2 = 23/30;
+qp45 = 170687698457070187/149820953007882240;
+qp31 = -3169112007572299/37455238251970560;
+
+% Q+ and Q-, top-left corner
+QpL = [...
+-216175739231516543/245800001028556800, 143647626126052207/149820953007882240, -67953322736279/2340952390748160, -1056618403248889/17835827739033600, -239061311631881/131093333881896960, 4843362440877649/449462859023646720;
+ -1008199398207353/6827777806348800, -36400588465707331/37455238251970560, 444281115722813/468190478149632, 823082791653949/4458956934758400, 492703999899311/32773333470474240, -1088192686015471/37455238251970560;
+ -137012145866261/13655555612697600, 1658590350384343/24970158834647040, -729260819980499/780317463582720, 2826535075774271/2972637956505600, -372310079307409/4369777796063232, 312883915256963/24970158834647040;
+ 3572646778087103/61450000257139200, -3169112007572299/37455238251970560, 17518644727427/2340952390748160, -4809914536549211/4458956934758400, 36544940032227659/32773333470474240, -2480203967217637/112365714755911680;
+ -572605772036807/27311111225395200, 935926732952051/29964190601576448, 4202889834071/585238097687040, 128340370711031/17835827739033600, -144662434693018337/131093333881896960, 170687698457070187/149820953007882240
+];
+QmL = [...
+-29624261797040257/245800001028556800, 1008199398207353/6827777806348800, 137012145866261/13655555612697600, -3572646778087103/61450000257139200, 572605772036807/27311111225395200;
+ -143647626126052207/149820953007882240, 36400588465707331/37455238251970560, -1658590350384343/24970158834647040, 3169112007572299/37455238251970560, -935926732952051/29964190601576448;
+ 67953322736279/2340952390748160, -444281115722813/468190478149632, 729260819980499/780317463582720, -17518644727427/2340952390748160, -4202889834071/585238097687040;
+ 1056618403248889/17835827739033600, -823082791653949/4458956934758400, -2826535075774271/2972637956505600, 4809914536549211/4458956934758400, -128340370711031/17835827739033600;
+ 239061311631881/131093333881896960, -492703999899311/32773333470474240, 372310079307409/4369777796063232, -36544940032227659/32773333470474240, 144662434693018337/131093333881896960;
+ -4843362440877649/449462859023646720, 1088192686015471/37455238251970560, -312883915256963/24970158834647040, 2480203967217637/112365714755911680, -170687698457070187/149820953007882240
+];
+
+% Q+ and Q-
+w = 6; 
+s = [ -3/640,  25/384,  -75/64,  75/64,  -25/384,  3/640];  
+Qp = spdiags(repmat(-s(end:-1:1),[n+2 1]), -(w/2-1):w/2, n+2, n+2); 
+Qm = spdiags(repmat(s(:)',[n+2 1]), -(w/2-1)-1:w/2-1, n+2, n+2);
+Qp(end,:) = [];
+Qm(:,end) = [];
+
+% Add SBP boundary closures
+bp = 5;
+bm = 6;
+Qp(1:bp,1:bm) = double(QpL);
+Qp(end-bp+1:end,end-bm+1:end) = -fliplr(flipud(QpL));
+Qm(1:bm,1:bp) = double(QmL);
+Qm(end-bm+1:end,end-bp+1:end) = -fliplr(flipud(QmL));
+
+% P+ and P-
+Pp = ones(n+1,1);
+Pm = ones(n+2,1);
+
+Pp(1:bp) = [95/288,  317/240,  23/30,  793/720,  157/160]; 
+Pp(end-bp+1:end) = Pp(bp:-1:1);
+Pm(1:bm) = [5251/68040,  133801/138240,  90101/103680,  26369/23040,  224201/241920,  1262401/1244160];
+Pm(end-bm+1:end) = Pm(bm:-1:1);
+Pp = spdiags(Pp,0,n+1,n+1);
+Pm = spdiags(Pm,0,n+2,n+2);
+
+Pp = h*Pp;
+Pm = h*Pm;
+
+xp = h*[0:n]';
+xm = h*[0 1/2+0:n n]';
\ No newline at end of file
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+sbp/D1Staggered.m	Tue Nov 07 20:44:13 2017 -0800
@@ -0,0 +1,83 @@
+classdef D1Staggered < sbp.OpSet
+    properties
+        % xp: "plus" grid with m points
+        % xm: "minus" grid with m+1 points
+
+        D1p % SBP operator approximating first derivative
+        D1m % SBP operator approximating first derivative
+        Hp % Norm matrix
+        Hm % Norm matrix
+        HpI % H^-1
+        HmI % H^-1
+        ep_l % Left boundary operator
+        em_l % Left boundary operator
+        ep_r % Right boundary operator
+        em_r % Right boundary operator
+        m % Number of grid points.
+        mp % Number of grid points.
+        mm % Number of grid points.
+        h % Step size
+        xp % grid
+        xm % grid
+        x
+        borrowing % Struct with borrowing limits for different norm matrices
+    end
+
+    methods
+        function obj = D1Staggered(m,lim,order)
+
+          x_l = lim{1};
+          x_r = lim{2};
+          L = x_r-x_l;
+
+          mp = m;
+          mm = m+1;
+
+          switch order
+          case 2
+            [xp, xm, Pp, Pm, Qp, Qm] = sbp.implementations.d1_staggered_2(m, L);
+          case 4
+            [xp, xm, Pp, Pm, Qp, Qm] = sbp.implementations.d1_staggered_4(m, L);
+          case 6
+            [xp, xm, Pp, Pm, Qp, Qm] = sbp.implementations.d1_staggered_6(m, L);
+          otherwise
+           error('Invalid operator order %d.',order);
+          end
+
+          obj.m = m;
+          obj.mp = mp;
+          obj.mm = mm;
+          obj.xp = x_l + xp';
+          obj.xm = x_l + xm';
+
+          D1p = Pp\Qp;
+          D1m = Pm\Qm;
+
+          obj.D1p = D1p;
+          obj.D1m = D1m;
+          pbj.Hp = Pp;
+          obj.Hm = Pm; 
+
+          obj.ep_l = sparse(mp,1);
+          obj.ep_r = sparse(mp,1);
+          obj.ep_l(1) = 1;
+          obj.ep_r(mp) = 1;
+
+          obj.em_l = sparse(mm,1);
+          obj.em_r = sparse(mm,1);
+          obj.em_l(1) = 1;
+          obj.em_r(mm) = 1;
+
+          obj.HpI = inv(obj.Hp);
+          obj.HmI = inv(obj.Hm);
+
+          obj.borrowing = [];
+          obj.x = [];
+
+        end
+
+        function str = string(obj)
+            str = [class(obj) '_' num2str(obj.order)];
+        end
+    end
+end