Mercurial > repos > public > sbplib

changeset 323:7579c2abbf9f feature/beams

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Cleaning d4_variable_8_higher_boundary_order.m
author Jonatan Werpers <jonatan@werpers.com>
date Mon, 26 Sep 2016 09:19:08 +0200
parents def409c10800
children c0cbffcf6513
files +sbp/+implementations/d4_variable_8_higher_boundary_order.m
diffstat 1 files changed, 31 insertions(+), 47 deletions(-) [+]
line wrap: on
line diff
--- a/+sbp/+implementations/d4_variable_8_higher_boundary_order.m	Mon Sep 26 09:12:08 2016 +0200
+++ b/+sbp/+implementations/d4_variable_8_higher_boundary_order.m	Mon Sep 26 09:19:08 2016 +0200
@@ -17,59 +17,44 @@
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     %This is 3rd order accurate at the boundary. Not same norm as D1 operator
 
-    H = diag(ones(m,1),0);
-    H(1:8,1:8) = [
-        0.7488203e7/0.25401600e8 0 0 0 0 0 0 0;
-        0 0.5539027e7/0.3628800e7 0 0 0 0 0 0;
-        0 0 0.308923e6/0.1209600e7 0 0 0 0 0;
-        0 0 0 0.1307491e7/0.725760e6 0 0 0 0;
-        0 0 0 0 0.59407e5/0.145152e6 0 0 0;
-        0 0 0 0 0 0.1548947e7/0.1209600e7 0 0;
-        0 0 0 0 0 0 0.3347963e7/0.3628800e7 0;
-        0 0 0 0 0 0 0 0.25641187e8/0.25401600e8;
-    ];
+    BP = 8;
+    if(m<2*BP)
+        error(['Operator requires at least ' num2str(2*BP) ' grid points']);
+    end
 
-    H(m-7:m,m-7:m) = fliplr(flipud(H(1:8,1:8)));
+    % Norm
+    Hv = ones(m,1);
+    Hv(1:8) = [0.7488203e7/0.25401600e8, 0.5539027e7/0.3628800e7, 0.308923e6/0.1209600e7, 0.1307491e7/0.725760e6, 0.59407e5/0.145152e6, 0.1548947e7/0.1209600e7, 0.3347963e7/0.3628800e7, 0.25641187e8/0.25401600e8];
+    Hv(m-7:m) = rot90(Hv(1:8),2);
+    Hv = h*Hv;
+    H = spdiag(Hv, 0);
+    HI = spdiag(1./Hv, 0);
 
-    e_1 = zeros(m,1);
-    e_1(1) = 1;
-    e_m = zeros(m,1);
-    e_m(m) = 1;
 
-    S_U = [-0.49e2/0.20e2 6 -0.15e2/0.2e1 0.20e2/0.3e1 -0.15e2/0.4e1 0.6e1/0.5e1 -0.1e1/0.6e1]/h;
-    S_1 = zeros(1,m);
-    S_1(1:7) = S_U;
-    S_m = zeros(1,m);
-    S_m(m-6:m) = fliplr(-S_U);
+    % Boundary operators
+    e_l = sparse(m,1);
+    e_l(1) = 1;
+    e_r = rot90(e_l, 2);
+
+    d1_l = sparse(m,1);
+    d1_l(1:7) = [-0.49e2/0.20e2 6 -0.15e2/0.2e1 0.20e2/0.3e1 -0.15e2/0.4e1 0.6e1/0.5e1 -0.1e1/0.6e1]/h;
+    d1_r = -rot90(d1_l);
 
-    S2_U = [0.203e3/0.45e2 -0.87e2/0.5e1 0.117e3/0.4e1 -0.254e3/0.9e1 0.33e2/0.2e1 -0.27e2/0.5e1 0.137e3/0.180e3]/h^2;
-    S2_1 = zeros(1,m);
-    S2_1(1:7) = S2_U;
-    S2_m = zeros(1,m);
-    S2_m(m-6:m) = fliplr(S2_U);
+    d2_l = sparse(m,1);
+    d2_l(1:7) = [0.203e3/0.45e2 -0.87e2/0.5e1 0.117e3/0.4e1 -0.254e3/0.9e1 0.33e2/0.2e1 -0.27e2/0.5e1 0.137e3/0.180e3]/h^2;
+    d2_r = rot90(d2_l, 2);
 
-    S3_U = [-0.49e2/0.8e1 29 -0.461e3/0.8e1 62 -0.307e3/0.8e1 13 -0.15e2/0.8e1]/h^3;
-    S3_1 = zeros(1,m);
-    S3_1(1:7) = S3_U;
-    S3_m = zeros(1,m);
-    S3_m(m-6:m) = fliplr(-S3_U);
+    d3_l = sparse(m,1);
+    d3_l(1:7) = [-0.49e2/0.8e1 29 -0.461e3/0.8e1 62 -0.307e3/0.8e1 13 -0.15e2/0.8e1]/h^3;
+    d3_r = -rot90(d3_l, 2);
 
-    H = h*H;
-    HI = inv(H);
+
 
     % Fourth derivative, 1th order accurate at first 8 boundary points (still
     % yield 5th order convergence if stable: for example u_tt = -u_xxxx
 
-    m5 = -0.41e2/0.7560e4;
-    m4 =  0.1261e4/0.15120e5;
-    m3 =  -0.541e3/0.840e3;
-    m2 =  0.4369e4/0.1260e4;
-    m1 =  -0.1669e4/0.180e3;
-    m0 =  0.1529e4/0.120e3;
-
-    M4 = m5*(diag(ones(m-5,1),5)+diag(ones(m-5,1),-5))+m4*(diag(ones(m-4,1),4)+diag(ones(m-4,1),-4))+m3*(diag(ones(m-3,1),3)+diag(ones(m-3,1),-3))+m2*(diag(ones(m-2,1),2)+diag(ones(m-2,1),-2))+m1*(diag(ones(m-1,1),1)+diag(ones(m-1,1),-1))+m0*diag(ones(m,1),0);
-
-    %M4 = (-1/6*(diag(ones(m-3,1),3)+diag(ones(m-3,1),-3) ) + 2*(diag(ones(m-2,1),2)+diag(ones(m-2,1),-2)) -13/2*(diag(ones(m-1,1),1)+diag(ones(m-1,1),-1)) + 28/3*diag(ones(m,1),0));
+    stencil = [-0.41e2/0.7560e4, 0.1261e4/0.15120e5, -0.541e3/0.840e3, 0.4369e4/0.1260e4, -0.1669e4/0.180e3, 0.1529e4/0.120e3, -0.1669e4/0.180e3, 0.4369e4/0.1260e4, -0.541e3/0.840e3, 0.1261e4/0.15120e5,-0.41e2/0.7560e4];
+    diags = -5:5;
 
     M4_U = [
         0.1031569831e10/0.155675520e9 -0.32874237931e11/0.1452971520e10 0.3069551773e10/0.90810720e8 -0.658395212131e12/0.21794572800e11 0.31068454007e11/0.1816214400e10 -0.39244130657e11/0.7264857600e10 0.1857767503e10/0.2724321600e10 0.1009939e7/0.49420800e8;
@@ -83,9 +68,8 @@
     ];
 
     M4(1:8,1:8) = M4_U;
+    M4(m-7:m,m-7:m) = rot90(M4_U, 2);
+    M4 = 1/h^3*M4;
 
-    M4(m-7:m,m-7:m) = flipud( fliplr( M4_U ) );
-    M4 = M4/h^3;
-
-    D4 = HI*(M4-e_1*S3_1+e_m*S3_m  + S_1'*S2_1-S_m'*S2_m);
+    D4=HI*(M4 - e_l*d3_l'+e_r*d3_r' + d1_l*d2_l'-d1_r*d2_r');
 end