Mercurial > repos > public > sbplib

changeset 47:ebd25e5a481a

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Added upwind operators.
author Jonatan Werpers <jonatan@werpers.com>
date Thu, 05 Nov 2015 16:30:17 -0800
parents dd4db0c97a02
children b21c53ff61d4
files +sbp/Upwind.m +sbp/upwind3_3bp.m +sbp/upwind4.m +sbp/upwind6.m +sbp/upwind8.m
diffstat 5 files changed, 152 insertions(+), 0 deletions(-) [+]
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+sbp/Upwind.m	Thu Nov 05 16:30:17 2015 -0800
@@ -0,0 +1,51 @@
+classdef Ordinary < sbp.OpSet
+    properties
+        norms % Struct containing norm matrices such as H,Q, M
+        boundary  % Struct contanging vectors for boundry point approximations
+        derivatives % Struct containging differentiation operators
+        borrowing % Struct with borrowing limits for different norm matrices
+        m % Number of grid points.
+        h % Step size
+    end
+
+    methods
+        function obj = Ordinary(m,h,order)
+
+            if order == 3
+                [H, HI, Dp, Dm, e_1, e_m] = sbp.upwind3(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);
+            end
+
+            obj.h = h;
+            obj.m = m;
+
+            obj.norms.H = H;
+            obj.norms.HI = HI;
+            obj.norms.Q = Q;
+
+            obj.boundary.e_1 = e_1;
+            obj.boundary.e_m = e_m;
+
+            obj.derivatives.Dp = Dp;
+            obj.derivatives.Dm = Dm;
+        end
+    end
+
+    methods (Static)
+        function lambda = smallestGrid(obj)
+            error('Not implmented')
+        end
+    end
+end
+
+
+
+
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+sbp/upwind3_3bp.m	Thu Nov 05 16:30:17 2015 -0800
@@ -0,0 +1,25 @@
+function [H, HI, Dp, Dm, e_1, e_m] = upwind3(m,h)
+    H(1:3,1:3)=[0.3e1 / 0.8e1 0 0; 0 0.7e1 / 0.6e1 0; 0 0 0.23e2 / 0.24e2;];
+    H(m-2:m,m-2:m)=fliplr(flipud(H(1:3,1:3)));
+    H=H*h;
+    H=spdiag(ones(m,1),0);
+    HI=inv(H);
+
+    Qp=(-1/3*spdiag(ones(m-1,1),-1)-1/2*spdiag(ones(m,1),0)+1*spdiag(ones(m-1,1),+1)-1/6*spdiag(ones(m-2,1),+2));
+    %Q_U = [-0.49e2 / 0.2640e4 0.11167e5 / 0.15840e5 -0.1541e4 / 0.7920e4 0.43e2 / 0.5280e4; -0.9403e4 / 0.15840e5 -0.551e3 / 0.2640e4 0.1779e4 / 0.1760e4 -0.2311e4 / 0.7920e4; 0.659e3 / 0.7920e4 -0.751e3 / 0.1760e4 -0.1541e4 / 0.2640e4 0.21287e5 / 0.15840e5; 0.51e2 / 0.1760e4 -0.551e3 / 0.7920e4 -0.3683e4 / 0.15840e5 -0.713e3 / 0.880e3;];
+    %Q_U = [-0.29e2 / 0.1680e4 0.7067e4 / 0.10080e5 -0.961e3 / 0.5040e4 0.23e2 / 0.3360e4; -0.6023e4 / 0.10080e5 -0.331e3 / 0.1680e4 0.1119e4 / 0.1120e4 -0.1451e4 / 0.5040e4; 0.439e3 / 0.5040e4 -0.491e3 / 0.1120e4 -0.961e3 / 0.1680e4 0.13507e5 / 0.10080e5; 0.31e2 / 0.1120e4 -0.331e3 / 0.5040e4 -0.2383e4 / 0.10080e5 -0.453e3 / 0.560e3;];
+
+    Q_U = [-0.1e1 / 0.24e2 0.17e2 / 0.24e2 -0.1e1 / 0.6e1; -0.13e2 / 0.24e2 -0.1e1 / 0.4e1 0.23e2 / 0.24e2; 0.1e1 / 0.12e2 -0.11e2 / 0.24e2 -0.11e2 / 0.24e2;];
+    Qp(1:3,1:3)=Q_U;
+    Qp(m-2:m,m-2:m)=flipud( fliplr(Q_U(1:3,1:3) ) )'; %%% 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/upwind4.m	Thu Nov 05 16:30:17 2015 -0800
@@ -0,0 +1,27 @@
+function [H, HI, Dp, Dm, e_1, e_m] = upwind4(m,h)
+
+    H=spdiag(ones(m,1),0);
+    H(1:4,1:4)=spdiag([49/144 61/48 41/48 149/144]);
+    H(m-3:m,m-3:m)=fliplr(flipud(spdiag([49/144 61/48 41/48 149/144])));
+    H=H*h;
+    HI=inv(H);
+
+
+    Qp=(-1/4*spdiag(ones(m-1,1),-1)-5/6*spdiag(ones(m,1),0)+3/2*spdiag(ones(m-1,1),+1)-1/2*spdiag(ones(m-2,1),+2)+1/12*spdiag(ones(m-3,1),+3));
+    %Q_U = [-0.49e2 / 0.2640e4 0.11167e5 / 0.15840e5 -0.1541e4 / 0.7920e4 0.43e2 / 0.5280e4; -0.9403e4 / 0.15840e5 -0.551e3 / 0.2640e4 0.1779e4 / 0.1760e4 -0.2311e4 / 0.7920e4; 0.659e3 / 0.7920e4 -0.751e3 / 0.1760e4 -0.1541e4 / 0.2640e4 0.21287e5 / 0.15840e5; 0.51e2 / 0.1760e4 -0.551e3 / 0.7920e4 -0.3683e4 / 0.15840e5 -0.713e3 / 0.880e3;];
+    %Q_U = [-0.29e2 / 0.1680e4 0.7067e4 / 0.10080e5 -0.961e3 / 0.5040e4 0.23e2 / 0.3360e4; -0.6023e4 / 0.10080e5 -0.331e3 / 0.1680e4 0.1119e4 / 0.1120e4 -0.1451e4 / 0.5040e4; 0.439e3 / 0.5040e4 -0.491e3 / 0.1120e4 -0.961e3 / 0.1680e4 0.13507e5 / 0.10080e5; 0.31e2 / 0.1120e4 -0.331e3 / 0.5040e4 -0.2383e4 / 0.10080e5 -0.453e3 / 0.560e3;];
+
+    Q_U = [-0.1e1 / 0.48e2 0.205e3 / 0.288e3 -0.29e2 / 0.144e3 0.1e1 / 0.96e2; -0.169e3 / 0.288e3 -0.11e2 / 0.48e2 0.33e2 / 0.32e2 -0.43e2 / 0.144e3; 0.11e2 / 0.144e3 -0.13e2 / 0.32e2 -0.29e2 / 0.48e2 0.389e3 / 0.288e3; 0.1e1 / 0.32e2 -0.11e2 / 0.144e3 -0.65e2 / 0.288e3 -0.13e2 / 0.16e2;];
+
+    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/upwind6.m	Thu Nov 05 16:30:17 2015 -0800
@@ -0,0 +1,25 @@
+function [H, HI, Dp, Dm, e_1, e_m] = upwind6(m,h)
+    H=spdiag(ones(m,1),0);
+    H(1:6,1:6)=[0.13613e5 / 0.43200e5 0 0 0 0 0; 0 0.12049e5 / 0.8640e4 0 0 0 0; 0 0 0.535e3 / 0.864e3 0 0 0; 0 0 0 0.1079e4 / 0.864e3 0 0; 0 0 0 0 0.7841e4 / 0.8640e4 0; 0 0 0 0 0 0.43837e5 / 0.43200e5;];
+    H(m-5:m,m-5:m)=fliplr(flipud(H(1:6,1:6)));
+    H=H*h;
+    HI=inv(H);
+
+
+    Qp=(1/30*spdiag(ones(m-2,1),-2)-2/5*spdiag(ones(m-1,1),-1)-7/12*spdiag(ones(m,1),0)+4/3*spdiag(ones(m-1,1),+1)-1/2*spdiag(ones(m-2,1),+2)+2/15*spdiag(ones(m-3,1),+3)-1/60*spdiag(ones(m-4,1),+4));
+    %Q_U = [-0.375523e6 / 0.237369744e9 0.73393569793e11 / 0.111599208000e12 -0.1244182786433e13 / 0.14954293872000e14 -0.245555527093e12 / 0.2492382312000e13 0.131952309637e12 / 0.14954293872000e14 0.62907543023e11 / 0.3738573468000e13; -0.4818422888131e13 / 0.7477146936000e13 -0.10098911e8 / 0.307701520e9 0.796148045807e12 / 0.1495429387200e13 0.487507785187e12 / 0.2990858774400e13 0.8887918703e10 / 0.249238231200e12 -0.804641655323e12 / 0.14954293872000e14; 0.925771462433e12 / 0.14954293872000e14 -0.601394962127e12 / 0.1495429387200e13 -0.674340047e9 / 0.4153970520e10 0.560418114917e12 / 0.747714693600e12 -0.1043208352973e13 / 0.2990858774400e13 0.42207619661e11 / 0.356054616000e12; 0.285728443093e12 / 0.2492382312000e13 -0.873874460227e12 / 0.2990858774400e13 -0.34960363091e11 / 0.106816384800e12 -0.1528791703e10 / 0.4153970520e10 0.1741615892927e13 / 0.1495429387200e13 -0.6107723941553e13 / 0.14954293872000e14; -0.213784143637e12 / 0.14954293872000e14 0.1011411511e10 / 0.35605461600e11 0.141331778093e12 / 0.2990858774400e13 -0.632571739487e12 / 0.1495429387200e13 -0.1464298201e10 / 0.2769313680e10 0.9523869708211e13 / 0.7477146936000e13; -0.60750297023e11 / 0.3738573468000e13 0.614597810123e12 / 0.14954293872000e14 -0.16980260027e11 / 0.2492382312000e13 -0.343200536047e12 / 0.14954293872000e14 -0.389757912373e12 / 0.1068163848000e13 -0.964053659e9 / 0.1661588208e10;];
+    Q_U =[-0.265e3 / 0.128688e6 0.1146190567e10 / 0.1737288000e10 -0.1596619e7 / 0.18384000e8 -0.55265831e8 / 0.579096000e9 0.26269819e8 / 0.3474576000e10 0.2464501e7 / 0.144774000e9; -0.1116490567e10 / 0.1737288000e10 -0.8839e4 / 0.214480e6 0.190538869e9 / 0.347457600e9 0.102705469e9 / 0.694915200e9 0.413741e6 / 0.9651600e7 -0.191689861e9 / 0.3474576000e10; 0.1096619e7 / 0.18384000e8 -0.135385429e9 / 0.347457600e9 -0.61067e5 / 0.321720e6 0.45137333e8 / 0.57909600e8 -0.253641811e9 / 0.694915200e9 0.70665929e8 / 0.579096000e9; 0.66965831e8 / 0.579096000e9 -0.208765789e9 / 0.694915200e9 -0.17623253e8 / 0.57909600e8 -0.18269e5 / 0.45960e5 0.410905829e9 / 0.347457600e9 -0.477953317e9 / 0.1158192000e10; -0.49219819e8 / 0.3474576000e10 0.293299e6 / 0.9651600e7 0.26422771e8 / 0.694915200e9 -0.141938309e9 / 0.347457600e9 -0.346583e6 / 0.643440e6 0.2217185207e10 / 0.1737288000e10; -0.2374501e7 / 0.144774000e9 0.142906261e9 / 0.3474576000e10 -0.3137129e7 / 0.579096000e9 -0.29884283e8 / 0.1158192000e10 -0.630168407e9 / 0.1737288000e10 -0.3559e4 / 0.6128e4;];
+
+    Qp(1:6,1:6)=Q_U;
+    Qp(m-5:m,m-5:m)=flipud( fliplr(Q_U(1:6,1:6) ) )'; %%% 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/upwind8.m	Thu Nov 05 16:30:17 2015 -0800
@@ -0,0 +1,24 @@
+function [H, HI, Dp, Dm, e_1, e_m] = upwind8(m,h)
+    H=spdiag(ones(m,1),0);
+    H(1:8,1:8)=[0.7489399e7 / 0.25401600e8 0 0 0 0 0 0 0; 0 0.5537831e7 / 0.3628800e7 0 0 0 0 0 0; 0 0 0.103373e6 / 0.403200e6 0 0 0 0 0; 0 0 0 0.261259e6 / 0.145152e6 0 0 0 0; 0 0 0 0 0.298231e6 / 0.725760e6 0 0 0; 0 0 0 0 0 0.515917e6 / 0.403200e6 0 0; 0 0 0 0 0 0 0.3349159e7 / 0.3628800e7 0; 0 0 0 0 0 0 0 0.25639991e8 / 0.25401600e8;];
+
+    H(m-7:m,m-7:m)=fliplr(flipud(H(1:8,1:8)));
+    H=H*h;
+    HI=inv(H);
+
+    Qp=(-1/168*spdiag(ones(m-3,1),-3)+1/14*spdiag(ones(m-2,1),-2)-1/2*spdiag(ones(m-1,1),-1)-9/20*spdiag(ones(m,1),0)+5/4*spdiag(ones(m-1,1),+1)-1/2*spdiag(ones(m-2,1),+2)+1/6*spdiag(ones(m-3,1),+3)-1/28*spdiag(ones(m-4,1),+4)+1/280*spdiag(ones(m-5,1),+5));
+    %Q_U = [-0.375523e6 / 0.237369744e9 0.73393569793e11 / 0.111599208000e12 -0.1244182786433e13 / 0.14954293872000e14 -0.245555527093e12 / 0.2492382312000e13 0.131952309637e12 / 0.14954293872000e14 0.62907543023e11 / 0.3738573468000e13; -0.4818422888131e13 / 0.7477146936000e13 -0.10098911e8 / 0.307701520e9 0.796148045807e12 / 0.1495429387200e13 0.487507785187e12 / 0.2990858774400e13 0.8887918703e10 / 0.249238231200e12 -0.804641655323e12 / 0.14954293872000e14; 0.925771462433e12 / 0.14954293872000e14 -0.601394962127e12 / 0.1495429387200e13 -0.674340047e9 / 0.4153970520e10 0.560418114917e12 / 0.747714693600e12 -0.1043208352973e13 / 0.2990858774400e13 0.42207619661e11 / 0.356054616000e12; 0.285728443093e12 / 0.2492382312000e13 -0.873874460227e12 / 0.2990858774400e13 -0.34960363091e11 / 0.106816384800e12 -0.1528791703e10 / 0.4153970520e10 0.1741615892927e13 / 0.1495429387200e13 -0.6107723941553e13 / 0.14954293872000e14; -0.213784143637e12 / 0.14954293872000e14 0.1011411511e10 / 0.35605461600e11 0.141331778093e12 / 0.2990858774400e13 -0.632571739487e12 / 0.1495429387200e13 -0.1464298201e10 / 0.2769313680e10 0.9523869708211e13 / 0.7477146936000e13; -0.60750297023e11 / 0.3738573468000e13 0.614597810123e12 / 0.14954293872000e14 -0.16980260027e11 / 0.2492382312000e13 -0.343200536047e12 / 0.14954293872000e14 -0.389757912373e12 / 0.1068163848000e13 -0.964053659e9 / 0.1661588208e10;];
+    Q_U =[-0.16683e5 / 0.63018560e8 0.29345822969987e14 / 0.43354248537600e14 -0.2734625426591e13 / 0.40644608004000e14 -0.113480208109603e15 / 0.780376473676800e15 -0.830250230261e12 / 0.26012549122560e14 0.2500519492033e13 / 0.32515686403200e14 0.2235718279643e13 / 0.390188236838400e15 -0.388481888477e12 / 0.26543417472000e14; -0.29227665839987e14 / 0.43354248537600e14 -0.493793e6 / 0.63018560e8 0.8302717120817e13 / 0.26543417472000e14 0.3739408501537e13 / 0.9290196115200e13 0.2684481534461e13 / 0.13935294172800e14 -0.4450185662513e13 / 0.18580392230400e14 -0.1221838279381e13 / 0.37160784460800e14 0.90595000956023e14 / 0.1950941184192000e16; 0.2505689537591e13 / 0.40644608004000e14 -0.7312922392817e13 / 0.26543417472000e14 -0.69332623e8 / 0.1323389760e10 0.10994933811709e14 / 0.18580392230400e14 -0.9270952411151e13 / 0.18580392230400e14 0.3191238635141e13 / 0.20644880256000e14 0.4442211176987e13 / 0.92901961152000e14 -0.940661365031e12 / 0.32515686403200e14; 0.118016946570403e15 / 0.780376473676800e15 -0.4173878828737e13 / 0.9290196115200e13 -0.7990503962509e13 / 0.18580392230400e14 -0.207799621e9 / 0.1323389760e10 0.2044021560341e13 / 0.2477385630720e13 0.511197701761e12 / 0.18580392230400e14 0.1237681717213e13 / 0.13935294172800e14 -0.7784834666617e13 / 0.130062745612800e15; 0.68609076271e11 / 0.2364777192960e13 -0.2235651762161e13 / 0.13935294172800e14 0.6527681584751e13 / 0.18580392230400e14 -0.1115980068821e13 / 0.2477385630720e13 -0.55386253e8 / 0.189055680e9 0.3208334350649e13 / 0.3716078446080e13 -0.407569013461e12 / 0.844563283200e12 0.136474842626653e15 / 0.780376473676800e15; -0.2487637785013e13 / 0.32515686403200e14 0.4244231077313e13 / 0.18580392230400e14 -0.1550378843141e13 / 0.20644880256000e14 -0.5726967564961e13 / 0.18580392230400e14 -0.1017898941929e13 / 0.3716078446080e13 -0.526653889e9 / 0.1323389760e10 0.45241297077547e14 / 0.37160784460800e14 -0.2279608411897e13 / 0.5080576000500e13; -0.2164019088443e13 / 0.390188236838400e15 0.1263196075861e13 / 0.37160784460800e14 -0.6600697610987e13 / 0.92901961152000e14 0.556610591687e12 / 0.13935294172800e14 0.926842346471e12 / 0.9290196115200e13 -0.18757693936747e14 / 0.37160784460800e14 -0.584765899e9 / 0.1323389760e10 0.204646287449e12 / 0.168431424000e12; 0.387091928477e12 / 0.26543417472000e14 -0.90231551688023e14 / 0.1950941184192000e16 0.1032404418251e13 / 0.32515686403200e14 0.3502353445417e13 / 0.130062745612800e15 -0.15385068876253e14 / 0.780376473676800e15 0.262499068919e12 / 0.10161152001000e14 -0.867004691939e12 / 0.1852745664000e13 -0.85017967e8 / 0.189055680e9;];
+
+    Qp(1:8,1:8)=Q_U;
+    Qp(m-7:m,m-7:m)=flipud( fliplr(Q_U(1:8,1:8) ) )'; %%% 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