diff +sbp/+implementations/d1_upwind_6.m @ 307:fefb2f9884f7 feature/beams

Merge with default.
author Jonatan Werpers <jonatan@werpers.com>
date Fri, 23 Sep 2016 10:40:12 +0200
parents f7ac3cd6eeaa
children
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+sbp/+implementations/d1_upwind_6.m	Fri Sep 23 10:40:12 2016 +0200
@@ -0,0 +1,39 @@
+function [H, HI, Dp, Dm, e_1, e_m] = d1_upwind_6(m,h)
+    
+    if(m<12)
+        error('Operator requires at least 12 grid points');
+    end
+
+    Hv = ones(m,1);
+    Hv(1:6)=[0.13613e5/0.43200e5; 0.12049e5/0.8640e4 ; 0.535e3/0.864e3 ; 0.1079e4/0.864e3 ;  0.7841e4/0.8640e4 ; 0.43837e5/0.43200e5];
+    Hv(m-5:m)=rot90(Hv(1:6),2);
+    Hv = Hv*h;
+    H = spdiag(Hv,0);
+    HI = spdiag(1./Hv,0);
+
+    q_diags   = [-2 -1 0 1 2 3 4];
+    q_stencil = [1/30 -2/5 -7/12 4/3 -1/2 2/15 -1/60];
+    Qp = stripeMatrix(q_stencil, q_diags,m);
+
+    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)=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