Mercurial > repos > public > sbplib

changeset 119:c56437d097de

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
I added remaining upwind operators.
author Jonatan Werpers <jonatan@werpers.com>
date Tue, 15 Dec 2015 16:38:17 +0100
parents 5046ff7d13b8
children 8eb4e39df8a5
files +sbp/Upwind.m +sbp/upwind2.m +sbp/upwind3.m +sbp/upwind3_3bp.m +sbp/upwind5.m +sbp/upwind7.m +sbp/upwind9.m
diffstat 7 files changed, 189 insertions(+), 42 deletions(-) [+]
line wrap: on
line diff
--- a/+sbp/Upwind.m	Tue Dec 15 16:38:02 2015 +0100
+++ b/+sbp/Upwind.m	Tue Dec 15 16:38:17 2015 +0100
@@ -11,16 +11,25 @@
     methods
         function obj = Upwind(m,h,order)
 
-            if order == 3
-                [H, HI, Dp, Dm, e_1, e_m] = sbp.upwind3_3bp(m,h);
-            elseif order == 4
-                [H, HI, Dp, Dm, e_1, e_m] = sbp.upwind4(m,h);
-            elseif order == 6
-                [H, HI, Dp, Dm, e_1, e_m] = sbp.upwind6(m,h);
-            elseif order == 8
-                [H, HI, Dp, Dm, e_1, e_m] = sbp.upwind8(m,h);
-            else
-                error('Invalid operator order %d.',order);
+            switch order
+                case 2
+                    [H, HI, Dp, Dm, e_1, e_m] = sbp.upwind2(m,h);
+                case 3
+                    [H, HI, Dp, Dm, e_1, e_m] = sbp.upwind3(m,h);
+                case 4
+                    [H, HI, Dp, Dm, e_1, e_m] = sbp.upwind4(m,h);
+                case 5
+                    [H, HI, Dp, Dm, e_1, e_m] = sbp.upwind5(m,h);
+                case 6
+                    [H, HI, Dp, Dm, e_1, e_m] = sbp.upwind6(m,h);
+                case 7
+                    [H, HI, Dp, Dm, e_1, e_m] = sbp.upwind7(m,h);
+                case 8
+                    [H, HI, Dp, Dm, e_1, e_m] = sbp.upwind8(m,h);
+                case 9
+                    [H, HI, Dp, Dm, e_1, e_m] = sbp.upwind9(m,h);
+                otherwise
+                    error('Invalid operator order %d.',order);
             end
 
             obj.h = h;
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+sbp/upwind2.m	Tue Dec 15 16:38:17 2015 +0100
@@ -0,0 +1,34 @@
+function [H, HI, Dp, Dm, e_1, e_m] = upwind2(m,h)
+    Hv=ones(m,1);
+    Hv(1:2)=[0.1e1/0.4e1; 0.5e1/0.4e1;];
+    Hv(m-1:m)=rot90(Hv(1:2),2);
+    Hv = Hv*h;
+    H = spdiag(Hv,0);
+    HI = spdiag(1./Hv,0);
+
+    q_diags   = [0 1 2];
+    q_stencil = [-3/2 +2 -1/2];
+    Qp = stripeMatrix(q_stencil, q_diags,m);
+
+    Q_U = [
+        -0.1e1/0.4e1 0.5e1/0.4e1;
+        -0.1e1/0.4e1 -0.5e1/0.4e1;
+    ];
+
+    Qp(1:2,1:2)=Q_U;
+    Qp(m-1:m,m-1:m)=rot90(Q_U,2)'; %%% This is different from standard SBP
+
+    Qm=-Qp';
+
+    e_1=zeros(m,1);e_1(1)=1;
+    e_m=zeros(m,1);e_m(m)=1;
+
+    Dp=HI*(Qp-1/2*e_1*e_1'+1/2*e_m*e_m') ;
+
+    Dm=HI*(Qm-1/2*e_1*e_1'+1/2*e_m*e_m') ;
+end
+
+
+
+% KOLLA NAMN PÅ ALLA FUNKTIONER
+% TESTA ALLA FUNKTIONER
\ No newline at end of file
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+sbp/upwind3.m	Tue Dec 15 16:38:17 2015 +0100
@@ -0,0 +1,32 @@
+function [H, HI, Dp, Dm, e_1, e_m] = upwind3(m,h)
+    Hv = ones(m,1);
+    Hv(1:3) = [3/8; 7/6; 23/24];
+    Hv(m-2:m) = rot90(Hv(1:3),2);
+    Hv = Hv*h;
+    H = spdiag(Hv,0);
+    HI = spdiag(1./Hv,0);
+
+    q_diags   = [-1, 0, 1, 2];
+    q_stencil = [-1/3 -1/2 1 -1/6];
+    Qp = stripeMatrix(q_stencil, q_diags,m);
+
+    Q_U = [
+         -1/24  17/24   -1/6;
+        -13/24   -1/4  23/24;
+          1/12 -11/24 -11/24;
+    ];
+    Qp(1:3,1:3)=Q_U;
+    Qp(m-2:m,m-2:m)=rot90(Q_U,2)'; %%% This is different from standard SBP
+
+    Qm=-Qp';
+
+    e_1=sparse(m,1);
+    e_1(1)=1;
+    e_m=sparse(m,1);
+    e_m(m)=1;
+
+    Dp=HI*(Qp-1/2*e_1*e_1'+1/2*e_m*e_m') ;
+
+    Dm=HI*(Qm-1/2*e_1*e_1'+1/2*e_m*e_m') ;
+end
+
--- a/+sbp/upwind3_3bp.m	Tue Dec 15 16:38:02 2015 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,32 +0,0 @@
-function [H, HI, Dp, Dm, e_1, e_m] = upwind3(m,h)
-    Hv = ones(m,1);
-    Hv(1:3) = [3/8; 7/6; 23/24];
-    Hv(m-2:m) = rot90(Hv(1:3),2);
-    Hv = Hv*h;
-    H = spdiag(Hv,0);
-    HI = spdiag(1./Hv,0);
-
-    q_diags   = [-1, 0, 1, 2];
-    q_stencil = [-1/3 -1/2 1 -1/6];
-    Qp = stripeMatrix(q_stencil, q_diags,m);
-
-    Q_U = [
-         -1/24  17/24   -1/6;
-        -13/24   -1/4  23/24;
-          1/12 -11/24 -11/24;
-    ];
-    Qp(1:3,1:3)=Q_U;
-    Qp(m-2:m,m-2:m)=rot90(Q_U,2)'; %%% This is different from standard SBP
-
-    Qm=-Qp';
-
-    e_1=sparse(m,1);
-    e_1(1)=1;
-    e_m=sparse(m,1);
-    e_m(m)=1;
-
-    Dp=HI*(Qp-1/2*e_1*e_1'+1/2*e_m*e_m') ;
-
-    Dm=HI*(Qm-1/2*e_1*e_1'+1/2*e_m*e_m') ;
-end
-
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+sbp/upwind5.m	Tue Dec 15 16:38:17 2015 +0100
@@ -0,0 +1,31 @@
+function [H, HI, Dp, Dm, e_1, e_m] = upwind5(m,h)
+    Hv=ones(m,1);
+    Hv(1:4)=[0.251e3/0.720e3; 0.299e3/0.240e3; 0.211e3/0.240e3; 0.739e3/0.720e3;];
+    Hv(m-3:m)=rot90(Hv(1:4),2);
+    Hv = Hv*h;
+    H = spdiag(Hv,0);
+    HI = spdiag(1./Hv,0);
+
+    q_diags   = [-2 -1 0 1 2 3];
+    q_stencil = [1/20 -1/2 -1/3 +1 -1/4 +1/30];
+    Qp = stripeMatrix(q_stencil, q_diags,m);
+
+    Q_U = [
+        -0.1e1/0.120e3 0.941e3/0.1440e4 -0.47e2/0.360e3 -0.7e1/0.480e3;
+        -0.869e3/0.1440e4 -0.11e2/0.120e3 0.25e2/0.32e2 -0.43e2/0.360e3;
+        0.29e2/0.360e3 -0.17e2/0.32e2 -0.29e2/0.120e3 0.1309e4/0.1440e4;
+        0.1e1/0.32e2 -0.11e2/0.360e3 -0.661e3/0.1440e4 -0.13e2/0.40e2;
+    ];
+
+    Qp(1:4,1:4)=Q_U;
+    Qp(m-3:m,m-3:m)=flipud( fliplr(Q_U(1:4,1:4) ) )'; %%% This is different from standard SBP
+
+    Qm=-Qp';
+
+    e_1=zeros(m,1);e_1(1)=1;
+    e_m=zeros(m,1);e_m(m)=1;
+
+    Dp=HI*(Qp-1/2*e_1*e_1'+1/2*e_m*e_m') ;
+
+    Dm=HI*(Qm-1/2*e_1*e_1'+1/2*e_m*e_m') ;
+end
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+sbp/upwind7.m	Tue Dec 15 16:38:17 2015 +0100
@@ -0,0 +1,35 @@
+function [H, HI, Dp, Dm, e_1, e_m] = upwind7(m,h)
+
+    Hv=ones(m,1);
+    Hv(1:6) = [0.19087e5/0.60480e5; 0.84199e5/0.60480e5; 0.18869e5/0.30240e5; 0.37621e5/0.30240e5; 0.55031e5/0.60480e5; 0.61343e5/0.60480e5; ];
+    Hv(m-5:m)=rot90(Hv(1:6),2);
+    Hv = Hv*h;
+    H = spdiag(Hv,0);
+    HI = spdiag(1./Hv,0);
+
+    q_diags = [-3 -2 -1 0 1 2 3 4];
+    q_stencil = [-1/105 +1/10 -3/5 -1/4 +1 -3/10 +1/15 -1/140];
+    Qp = stripeMatrix(q_stencil, q_diags,m);
+
+    Q_U =[
+        -0.265e3/0.300272e6 0.1587945773e10/0.2432203200e10 -0.1926361e7/0.25737600e8 -0.84398989e8/0.810734400e9 0.48781961e8/0.4864406400e10 0.3429119e7/0.202683600e9;
+        -0.1570125773e10/0.2432203200e10 -0.26517e5/0.1501360e7 0.240029831e9/0.486440640e9 0.202934303e9/0.972881280e9 0.118207e6/0.13512240e8 -0.231357719e9/0.4864406400e10;
+        0.1626361e7/0.25737600e8 -0.206937767e9/0.486440640e9 -0.61067e5/0.750680e6 0.49602727e8/0.81073440e8 -0.43783933e8/0.194576256e9 0.51815011e8/0.810734400e9;
+        0.91418989e8/0.810734400e9 -0.53314099e8/0.194576256e9 -0.33094279e8/0.81073440e8 -0.18269e5/0.107240e6 0.440626231e9/0.486440640e9 -0.365711063e9/0.1621468800e10;
+        -0.62551961e8/0.4864406400e10 0.799e3/0.35280e5 0.82588241e8/0.972881280e9 -0.279245719e9/0.486440640e9 -0.346583e6/0.1501360e7 0.2312302333e10/0.2432203200e10;
+        -0.3375119e7/0.202683600e9 0.202087559e9/0.4864406400e10 -0.11297731e8/0.810734400e9 0.61008503e8/0.1621468800e10 -0.1360092253e10/0.2432203200e10 -0.10677e5/0.42896e5;
+    ];
+
+
+    Qp(1:6,1:6)=Q_U;
+    Qp(m-5:m,m-5:m)=rot90(Q_U,2)';
+
+    Qm=-Qp';
+
+    e_1=zeros(m,1);e_1(1)=1;
+    e_m=zeros(m,1);e_m(m)=1;
+
+    Dp=HI*(Qp-1/2*e_1*e_1'+1/2*e_m*e_m') ;
+
+    Dm=HI*(Qm-1/2*e_1*e_1'+1/2*e_m*e_m') ;
+end
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+sbp/upwind9.m	Tue Dec 15 16:38:17 2015 +0100
@@ -0,0 +1,38 @@
+function [H, HI, Dp, Dm, e_1, e_m] = upwind9(m,h)
+
+    Hv=ones(m,1);
+    Hv(1:8) = [0.1070017e7/0.3628800e7; 0.5537111e7/0.3628800e7; 0.103613e6/0.403200e6; 0.261115e6/0.145152e6; 0.298951e6/0.725760e6; 0.515677e6/0.403200e6; 0.3349879e7/0.3628800e7; 0.3662753e7/0.3628800e7;];
+    Hv(m-7:m)=rot90(Hv(1:8),2);
+    Hv=Hv*h;
+
+    H = spdiag(Hv,0);
+    HI = spdiag(1./Hv,0);
+
+    q_diags   = [-4 -3 -2 -1 0 1 2 3 4 5];
+    q_stencil = [1/504 -1/42 +1/7 -2/3 -1/5 +1 -1/3 +2/21 -1/56 +1/630];
+    Qp = stripeMatrix(q_stencil, q_diags,m);
+
+    Q_U =[
+        -0.5561e4/0.47263920e8 0.4186300102421e13/0.6193464076800e13 -0.377895002003e12/0.5806372572000e13 -0.16485548951749e14/0.111482353382400e15 -0.113245973003e12/0.3716078446080e13 0.355360297339e12/0.4645098057600e13 0.321012170669e12/0.55741176691200e14 -0.388397049437e12/0.26543417472000e14;
+        -0.4178798062421e13/0.6193464076800e13 -0.493793e6/0.141791760e9 0.725405227507e12/0.2413037952000e13 0.3904159533697e13/0.9290196115200e13 0.2483046570341e13/0.13935294172800e14 -0.4336328670953e13/0.18580392230400e14 -0.1258688487061e13/0.37160784460800e14 0.12931584852209e14/0.278705883456000e15;
+        0.363359390003e12/0.5806372572000e13 -0.7539548734577e13/0.26543417472000e14 -0.69332623e8/0.2977626960e10 0.9994352248429e13/0.18580392230400e14 -0.8195655811631e13/0.18580392230400e14 0.7361486640463e13/0.61934640768000e14 0.5539855071347e13/0.92901961152000e14 -0.12898722943e11/0.422281641600e12;
+        0.16773595838149e14/0.111482353382400e15 -0.372477950627e12/0.844563283200e12 -0.8659050093229e13/0.18580392230400e14 -0.207799621e9/0.2977626960e10 0.1734921317461e13/0.2477385630720e13 0.2530020015841e13/0.18580392230400e14 0.441856623253e12/0.13935294172800e14 -0.115132773073e12/0.2654341747200e13;
+        0.108449122763e12/0.3716078446080e13 -0.2283566671541e13/0.13935294172800e14 0.6976424333231e13/0.18580392230400e14 -0.440819477447e12/0.825795210240e12 -0.55386253e8/0.425375280e9 0.2479572560009e13/0.3716078446080e13 -0.40258468963e11/0.120651897600e12 0.11808221047099e14/0.111482353382400e15;
+        -0.32231128289e11/0.422281641600e12 0.4244793299753e13/0.18580392230400e14 -0.5173673584463e13/0.61934640768000e14 -0.4848139955041e13/0.18580392230400e14 -0.1506045711689e13/0.3716078446080e13 -0.526653889e9/0.2977626960e10 0.36411368691307e14/0.37160784460800e14 -0.825434105779e12/0.2903186286000e13;
+        -0.316459841069e12/0.55741176691200e14 0.1277069729941e13/0.37160784460800e14 -0.6499182375347e13/0.92901961152000e14 0.355606625147e12/0.13935294172800e14 0.1519272420551e13/0.9290196115200e13 -0.2240079855137e13/0.3378253132800e13 -0.584765899e9/0.2977626960e10 0.2301241355533e13/0.2382101568000e13;
+        0.387779289437e12/0.26543417472000e14 -0.12908508708209e14/0.278705883456000e15 0.147710908133e12/0.4645098057600e13 0.534025841911e12/0.18580392230400e14 -0.4119981443899e13/0.111482353382400e15 0.279819152779e12/0.2903186286000e13 -0.1510324515533e13/0.2382101568000e13 -0.85017967e8/0.425375280e9;
+    ];
+
+
+    Qp(1:8,1:8)=Q_U;
+    Qp(m-7:m,m-7:m)=rot90(Q_U,2)'; %%% This is different from standard SBP
+
+    Qm=-Qp';
+
+    e_1=zeros(m,1);e_1(1)=1;
+    e_m=zeros(m,1);e_m(m)=1;
+
+    Dp=HI*(Qp-1/2*e_1*e_1'+1/2*e_m*e_m') ;
+
+    Dm=HI*(Qm-1/2*e_1*e_1'+1/2*e_m*e_m') ;
+end