Mercurial > repos > public > sbplib

changeset 314:88584b0cfba1 feature/beams

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Corrections and clean up order 4
author Jonatan Werpers <jonatan@werpers.com>
date Fri, 23 Sep 2016 22:55:30 +0200
parents 52b4cdf27633
children 297d2cbfbe15
files +sbp/+implementations/d1_upwind_4.m +sbp/+implementations/d4_variable_2.m +sbp/+implementations/d4_variable_4.m
diffstat 3 files changed, 49 insertions(+), 83 deletions(-) [+]
line wrap: on
line diff
--- a/+sbp/+implementations/d1_upwind_4.m	Fri Sep 23 21:58:52 2016 +0200
+++ b/+sbp/+implementations/d1_upwind_4.m	Fri Sep 23 22:55:30 2016 +0200
@@ -23,7 +23,7 @@
     ];
 
     Qp(1:4,1:4)=Q_U;
-    Qp(m-3:m,m-3:m)=rot90(Q_U,2)'; %%% This is different from standard SBP
+    Qp(m-3:m,m-3:m)=rot90(Q_U, 2)'; %%% This is different from standard SBP
 
     Qm=-Qp';
 
--- a/+sbp/+implementations/d4_variable_2.m	Fri Sep 23 21:58:52 2016 +0200
+++ b/+sbp/+implementations/d4_variable_2.m	Fri Sep 23 22:55:30 2016 +0200
@@ -37,7 +37,7 @@
 
     d3_l = sparse(m,1);
     d3_l(1:4) = 1/h^3*[-1 3 -3 1];
-    d3_r = -rot90(d3_l, 2)
+    d3_r = -rot90(d3_l, 2);
 
 
     % First derivative SBP operator, 1st order accurate at first 6 boundary points
@@ -45,7 +45,7 @@
     diags = [-1 0 1];
     Q = stripeMatrix(stencil, diags, m);
 
-    D1 = HI*(Q-1/2*(e_1*e_1') + 1/2*(e_m*e_m'));
+    D1 = HI*(Q - 1/2*e_l*e_l' + 1/2*e_r*e_r');
 
     % Second derivative, 1st order accurate at first boundary points
     M = sparse(m,m);
@@ -85,5 +85,5 @@
     M4(m-3:m,m-3:m) = rot90(M4_U, 2);
     M4 = 1/h^3*M4;
 
-    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
--- a/+sbp/+implementations/d4_variable_4.m	Fri Sep 23 21:58:52 2016 +0200
+++ b/+sbp/+implementations/d4_variable_4.m	Fri Sep 23 22:55:30 2016 +0200
@@ -13,56 +13,53 @@
         error(['Operator requires at least ' num2str(2*BP) ' grid points']);
     end
 
-    H = speye(m,m);
-    H(1:4,1:4) = diag([17/48 59/48 43/48 49/48]);
-    H(m-3:m,m-3:m) = rot90(diag([17/48 59/48 43/48 49/48]),2);
-    H = H*h;
-    HI = inv(H);
-    HI = sparse(HI);
+    % Norm
+    Hv = ones(m,1);
+    Hv(1:4) = [17/48 59/48 43/48 49/48];
+    Hv(m-3:m) = rot90(Hv(1:4),2);
+    Hv = h*Hv;
+    H = spdiag(Hv, 0);
+    HI = spdiag(1./Hv, 0);
 
 
-%     Q=(-1/12*diag(ones(m-2,1),2)+8/12*diag(ones(m-1,1),1)-8/12*diag(ones(m-1,1),-1)+1/12*diag(ones(m-2,1),-2));
-      e = ones(m,1);
-%       Q=spdiags([e -8*e 0*e 8*e -e], -2:2, m, m)/12;
-%     Q_U = [0 0.59e2/0.96e2 -0.1e1/0.12e2 -0.1e1/0.32e2; -0.59e2/0.96e2 0 0.59e2/0.96e2 0; 0.1e1/0.12e2 -0.59e2/0.96e2 0 0.59e2/0.96e2; 0.1e1/0.32e2 0 -0.59e2/0.96e2 0;];
-%     Q(1:4,1:4)=Q_U;
-%     Q(m-3:m,m-3:m)=rot90( -Q_U(1:4,1:4) ,2 );
-
-    e_1 = sparse(m,1);
-    e_1(1) = 1;
-    e_m = sparse(m,1);
-    e_m(m) = 1;
-
-%     D1=HI*(Q-1/2*(e_1*e_1')+1/2*(e_m*e_m')) ;
+    % Boundary operators
+    e_l = sparse(m,1);
+    e_l(1) = 1;
+    e_r = rot90(e_l, 2);
 
-    M_U = [
-        0.9e1/0.8e1 -0.59e2/0.48e2 0.1e1/0.12e2 0.1e1/0.48e2;
-        -0.59e2/0.48e2 0.59e2/0.24e2 -0.59e2/0.48e2 0;
-        0.1e1/0.12e2 -0.59e2/0.48e2 0.55e2/0.24e2 -0.59e2/0.48e2;
-        0.1e1/0.48e2 0 -0.59e2/0.48e2 0.59e2/0.24e2;
-    ];
-%     M=-(-1/12*diag(ones(m-2,1),2)+16/12*diag(ones(m-1,1),1)+16/12*diag(ones(m-1,1),-1)-1/12*diag(ones(m-2,1),-2)-30/12*diag(ones(m,1),0));
-    M = -spdiags([-e 16*e -30*e 16*e -e], -2:2, m, m)/12;
-
-    M(1:4,1:4) = M_U;
+    d1_l = sparse(m,1);
+    d1_l(1:4) = 1/h*[-11/6 3 -3/2 1/3];
+    d1_r = -rot90(d1_l);
 
-    M(m-3:m,m-3:m) = rot90(  M_U ,2 );
-    M=M/h;
+    d2_l = sparse(m,1);
+    d2_l(1:4) = 1/h^2*[2 -5 4 -1];
+    d2_r = rot90(d2_l, 2);
 
-    S_U=[-0.11e2/0.6e1 3 -0.3e1/0.2e1 0.1e1/0.3e1;]/h;
-    S_1=sparse(1,m);
-    S_1(1:4)=S_U;
-    S_m=sparse(1,m);
-    S_m(m-3:m)=fliplr(-S_U);
-    S_1 = S_1';
-    S_m = S_m';
+    d3_l = sparse(m,1);
+    d3_l(1:4) = 1/h^3*[-1 3 -3 1];
+    d3_r = -rot90(d3_l, 2);
 
 
+    % First derivative SBP operator,
+    stencil = [1/12 -2/3 0 2/3 -1/12];
+    diags = [-1 0 1];
+
+    Q_U = [
+        0 0.59e2/0.96e2 -0.1e1/0.12e2 -0.1e1/0.32e2;
+        -0.59e2/0.96e2 0 0.59e2/0.96e2 0;
+        0.1e1/0.12e2 -0.59e2/0.96e2 0 0.59e2/0.96e2;
+        0.1e1/0.32e2 0 -0.59e2/0.96e2 0;
+    ];
+
+    Q = stripeMatrix(stencil, diags, m);
+    Q(1:4,1:4)=Q_U;
+    Q(m-3:m,m-3:m) = -rot90(Q_U, 2);
+
+    D1 = HI*(Q - 1/2*e_l*e_l' + 1/2*e_r*e_r');
+
+
+    % Second derivative
     M=sparse(m,m);
-    e_1 = sparse(e_1);
-    e_m = sparse(e_m);
-    S_1 = sparse(S_1);
-    S_m = sparse(S_m);
 
     scheme_width = 5;
     scheme_radius = (scheme_width-1)/2;
@@ -109,9 +106,6 @@
         % Kan man skriva det som en multiplikation av en 3-dim matris?
         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-
-
-
         M(1:6,1:6) = [
             0.12e2/0.17e2 * c(1) + 0.59e2/0.192e3 * c(2) + 0.27010400129e11/0.345067064608e12 * c(3) + 0.69462376031e11/0.2070402387648e13 * c(4) -0.59e2/0.68e2 * c(1) - 0.6025413881e10/0.21126554976e11 * c(3) - 0.537416663e9/0.7042184992e10 * c(4) 0.2e1/0.17e2 * c(1) - 0.59e2/0.192e3 * c(2) + 0.213318005e9/0.16049630912e11 * c(4) + 0.2083938599e10/0.8024815456e10 * c(3) 0.3e1/0.68e2 * c(1) - 0.1244724001e10/0.21126554976e11 * c(3) + 0.752806667e9/0.21126554976e11 * c(4) 0.49579087e8/0.10149031312e11 * c(3) - 0.49579087e8/0.10149031312e11 * c(4) -c(4)/0.784e3 + c(3)/0.784e3;
             -0.59e2/0.68e2 * c(1) - 0.6025413881e10/0.21126554976e11 * c(3) - 0.537416663e9/0.7042184992e10 * c(4) 0.3481e4/0.3264e4 * c(1) + 0.9258282831623875e16/0.7669235228057664e16 * c(3) + 0.236024329996203e15/0.1278205871342944e16 * c(4) -0.59e2/0.408e3 * c(1) - 0.29294615794607e14/0.29725717938208e14 * c(3) - 0.2944673881023e13/0.29725717938208e14 * c(4) -0.59e2/0.1088e4 * c(1) + 0.260297319232891e15/0.2556411742685888e16 * c(3) - 0.60834186813841e14/0.1278205871342944e16 * c(4) -0.1328188692663e13/0.37594290333616e14 * c(3) + 0.1328188692663e13/0.37594290333616e14 * c(4) -0.8673e4/0.2904112e7 * c(3) + 0.8673e4/0.2904112e7 * c(4);
@@ -136,26 +130,12 @@
     end
     D2 = @D2_fun;
 
-    S2_U=[2 -5 4 -1;]/h^2;
-    S2_1=sparse(1,m);
-    S2_1(1:4)=S2_U;
-    S2_m=sparse(1,m);
-    S2_m(m-3:m)=fliplr(S2_U);
-    S2_1 = S2_1';
-    S2_m = S2_m';
 
-    m3 = -1/6;
-    m2 = 2;
-    m1 = -13/2;
-    m0 = 28/3;
-%     M4=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);
-    stencil = [m3,m2,m1,m0,m1,m2,m3];
-    d = (length(stencil)-1)/2;
-    diags = -d:d;
+    % Fourth derivative
+    stencil = [-1/6,2,-13/2, 28/3,-13/2,2,-1/6];
+    diags = -3:3;
     M4 = stripeMatrix(stencil, diags, m);
 
-    %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));
-
     M4_U = [
         0.5762947e7/0.2316384e7 -0.6374287e7/0.1158192e7 0.573947e6/0.165456e6 -0.124637e6/0.289548e6 0.67979e5/0.2316384e7 -0.60257e5/0.1158192e7;
         -0.6374287e7/0.1158192e7 0.30392389e8/0.2316384e7 -0.2735053e7/0.289548e6 0.273109e6/0.165456e6 0.83767e5/0.1158192e7 0.245549e6/0.2316384e7;
@@ -166,22 +146,8 @@
     ];
 
     M4(1:6,1:6) = M4_U;
-
-    M4(m-5:m,m-5:m) = rot90(  M4_U ,2 );
-    M4 = M4/h^3;
+    M4(m-5:m,m-5:m) = rot90(M4_U, 2);
+    M4 = 1/h^3*M4;
 
-    S3_U = [-1 3 -3 1;]/h^3;
-    S3_1 = sparse(1,m);
-    S3_1(1:4)=S3_U;
-    S3_m = sparse(1,m);
-    S3_m(m-3:m) = fliplr(-S3_U);
-    S3_1 = S3_1';
-    S3_m = S3_m';
-
-    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