Mercurial > repos > public > sbplib

changeset 316:203afa156f59 feature/beams

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Collected boundary operators.
author Jonatan Werpers <jonatan@werpers.com>
date Fri, 23 Sep 2016 23:10:44 +0200
parents 297d2cbfbe15
children c7ac7e12de8a
files +sbp/+implementations/d4_variable_6.m +sbp/+implementations/d4_variable_6_2.m +sbp/+implementations/d4_variable_6_3.m +sbp/+implementations/d4_variable_6_min_boundary_points.m +sbp/+implementations/d4_variable_8_higher_boundary_order.m +sbp/+implementations/d4_variable_8_min_boundary_points.m
diffstat 6 files changed, 93 insertions(+), 122 deletions(-) [+]
line wrap: on
line diff
--- a/+sbp/+implementations/d4_variable_6.m	Fri Sep 23 22:58:45 2016 +0200
+++ b/+sbp/+implementations/d4_variable_6.m	Fri Sep 23 23:10:44 2016 +0200
@@ -16,15 +16,6 @@
     %%% DI=D4*B*H*D4                            %%%
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-
-    %m=10;          %problemstorlek
-    %h=1/(m-1);
-
-    % Variable koefficicients are stored in vector: c, size m,
-    % with the unknown stored as c(1), c(2), ..., c_m
-    % x=1:h:m*h;x=x';
-    % c=x.^0;
-
     BP = 8;
     if(m<2*BP)
         error(['Operator requires at least ' num2str(2*BP) ' grid points']);
@@ -76,7 +67,24 @@
     S_m(m-4:m)=fliplr(-S_U);
     S_1 = S_1';
     S_m = S_m';
-
+    e_1 = sparse(e_1);
+    e_m = sparse(e_m);
+    S_1 = sparse(S_1);
+    S_m = sparse(S_m);
+    S2_U=[0.35e2/0.12e2 -0.26e2/0.3e1 0.19e2/0.2e1 -0.14e2/0.3e1 0.11e2/0.12e2;]/h^2;
+    S2_1=sparse(1,m);
+    S2_1(1:5)=S2_U;
+    S2_m=sparse(1,m);
+    S2_m(m-4:m)=fliplr(S2_U);
+    S2_1 = S2_1';
+    S2_m = S2_m';
+    S3_U = [-5/2 9 -12 7 -3/2]/h^3;
+    S3_1 = sparse(1,m);
+    S3_1(1:5)=S3_U;
+    S3_m = sparse(1,m);
+    S3_m(m-4:m) = fliplr(-S3_U);
+    S3_1 = S3_1';
+    S3_m = S3_m';
 
 
     %DS=sparse(m,m);
@@ -90,10 +98,6 @@
 
 
     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 = 7;
     scheme_radius = (scheme_width-1)/2;
@@ -145,13 +149,7 @@
     end
     D2 = @D2_fun;
 
-    S2_U=[0.35e2/0.12e2 -0.26e2/0.3e1 0.19e2/0.2e1 -0.14e2/0.3e1 0.11e2/0.12e2;]/h^2;
-    S2_1=sparse(1,m);
-    S2_1(1:5)=S2_U;
-    S2_m=sparse(1,m);
-    S2_m(m-4:m)=fliplr(S2_U);
-    S2_1 = S2_1';
-    S2_m = S2_m';
+
 
 
 
@@ -189,13 +187,7 @@
     M4(m-7:m,m-7:m) = rot90(  M4_U ,2 );
     M4 = M4/h^3;
 
-    S3_U = [-5/2 9 -12 7 -3/2]/h^3;
-    S3_1 = sparse(1,m);
-    S3_1(1:5)=S3_U;
-    S3_m = sparse(1,m);
-    S3_m(m-4: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');
 
--- a/+sbp/+implementations/d4_variable_6_2.m	Fri Sep 23 22:58:45 2016 +0200
+++ b/+sbp/+implementations/d4_variable_6_2.m	Fri Sep 23 23:10:44 2016 +0200
@@ -42,12 +42,17 @@
     S_m = zeros(1,m);
     S_m(m-5:m) = fliplr(-S_U);
 
-
+    S2_U = [0.15e2/0.4e1 -0.77e2/0.6e1 0.107e3/0.6e1 -13 0.61e2/0.12e2 -0.5e1/0.6e1;]/h^2;
+    S2_1 = zeros(1,m);
+    S2_1(1:6) = S2_U;
+    S2_m = zeros(1,m);
+    S2_m(m-5:m) = fliplr(S2_U);
 
-    %DS = zeros(m,m);
-    %DS(1,1:5) = -[-25/12, 4, -3, 4/3, -1/4];
-    %DS(m,m-4:m) = fliplr(-[-25/12, 4, -3, 4/3, -1/4]);
-    %DS = diag(c)*DS/h;
+    S3_U = [-0.17e2/0.4e1 0.71e2/0.4e1 -0.59e2/0.2e1 0.49e2/0.2e1 -0.41e2/0.4e1 0.7e1/0.4e1;]/h^3;
+    S3_1 = zeros(1,m);
+    S3_1(1:6) = S3_U;
+    S3_m = zeros(1,m);
+    S3_m(m-5:m) = fliplr(-S3_U);
 
 
     H = h*H;
@@ -92,15 +97,6 @@
     %
     % D2=HI*(-M-diag(c)*e_1*S_1+diag(c)*e_m*S_m);
 
-    S2_U = [0.15e2/0.4e1 -0.77e2/0.6e1 0.107e3/0.6e1 -13 0.61e2/0.12e2 -0.5e1/0.6e1;]/h^2;
-    S2_1 = zeros(1,m);
-    S2_1(1:6) = S2_U;
-    S2_m = zeros(1,m);
-    S2_m(m-5:m) = fliplr(S2_U);
-
-
-
-
 
     % Fourth derivative, 1th order accurate at first 8 boundary points (still
     % yield 5th order convergence if stable: for example u_tt = -u_xxxx
@@ -129,13 +125,5 @@
     M4(m-5:m,m-5:m) = flipud( fliplr( M4_U ) );
     M4 = M4/h^3;
 
-    S3_U = [-0.17e2/0.4e1 0.71e2/0.4e1 -0.59e2/0.2e1 0.49e2/0.2e1 -0.41e2/0.4e1 0.7e1/0.4e1;]/h^3;
-    S3_1 = zeros(1,m);
-    S3_1(1:6) = S3_U;
-    S3_m = zeros(1,m);
-    S3_m(m-5:m) = fliplr(-S3_U);
-
     D4 = HI*(M4-e_1*S3_1+e_m*S3_m  + S_1'*S2_1-S_m'*S2_m);
-
-
 end
--- a/+sbp/+implementations/d4_variable_6_3.m	Fri Sep 23 22:58:45 2016 +0200
+++ b/+sbp/+implementations/d4_variable_6_3.m	Fri Sep 23 23:10:44 2016 +0200
@@ -48,6 +48,18 @@
     S_m(m-5:m)=fliplr(-S_U);
 
 
+    S2_U = [0.15e2/0.4e1 -0.77e2/0.6e1 0.107e3/0.6e1 -13 0.61e2/0.12e2 -0.5e1/0.6e1;]/h^2;
+    S2_1 = zeros(1,m);
+    S2_1(1:6) = S2_U;
+    S2_m = zeros(1,m);
+    S2_m(m-5:m) = fliplr(S2_U);
+
+
+    S3_U = [-0.17e2/0.4e1 0.71e2/0.4e1 -0.59e2/0.2e1 0.49e2/0.2e1 -0.41e2/0.4e1 0.7e1/0.4e1;]/h^3;
+    S3_1 = zeros(1,m);
+    S3_1(1:6) = S3_U;
+    S3_m = zeros(1,m);
+    S3_m(m-5:m) = fliplr(-S3_U);
 
     %DS=zeros(m,m);
     %DS(1,1:5)=-[-25/12, 4, -3, 4/3, -1/4];
@@ -97,15 +109,6 @@
     %
     % D2=HI*(-M-diag(c)*e_1*S_1+diag(c)*e_m*S_m);
 
-    S2_U = [0.15e2/0.4e1 -0.77e2/0.6e1 0.107e3/0.6e1 -13 0.61e2/0.12e2 -0.5e1/0.6e1;]/h^2;
-    S2_1 = zeros(1,m);
-    S2_1(1:6) = S2_U;
-    S2_m = zeros(1,m);
-    S2_m(m-5:m) = fliplr(S2_U);
-
-
-
-
 
     % Fourth derivative, 1th order accurate at first 8 boundary points (still
     % yield 5th order convergence if stable: for example u_tt=-u_xxxx
@@ -135,11 +138,5 @@
     M4(m-6:m,m-6:m) = flipud( fliplr( M4_U ) );
     M4 = M4/h^3;
 
-    S3_U = [-0.17e2/0.4e1 0.71e2/0.4e1 -0.59e2/0.2e1 0.49e2/0.2e1 -0.41e2/0.4e1 0.7e1/0.4e1;]/h^3;
-    S3_1 = zeros(1,m);
-    S3_1(1:6) = S3_U;
-    S3_m = zeros(1,m);
-    S3_m(m-5:m) = fliplr(-S3_U);
-
     D4 = HI*(M4-e_1*S3_1+e_m*S3_m  + S_1'*S2_1-S_m'*S2_m);
 end
--- a/+sbp/+implementations/d4_variable_6_min_boundary_points.m	Fri Sep 23 22:58:45 2016 +0200
+++ b/+sbp/+implementations/d4_variable_6_min_boundary_points.m	Fri Sep 23 23:10:44 2016 +0200
@@ -22,6 +22,27 @@
     H(1:6,1:6) = diag([13649/43200,12013/8640,2711/4320,5359/4320,7877/8640, 43801/43200]);
     H(m-5:m,m-5:m) = fliplr(flipud(diag([13649/43200,12013/8640, 2711/4320,5359/4320,7877/8640,43801/43200])));
 
+    e_1 = zeros(m,1);e_1(1)=1;
+    e_m = zeros(m,1);e_m(m)=1;
+
+    S_U = [-25/12, 4, -3, 4/3, -1/4]/h;
+    S_1 = zeros(1,m);
+    S_1(1:5) = S_U;
+    S_m = zeros(1,m);
+    S_m(m-4:m) = fliplr(-S_U);
+
+    S2_U = [0.35e2/0.12e2 -0.26e2/0.3e1 0.19e2/0.2e1 -0.14e2/0.3e1 0.11e2/0.12e2;]/h^2;
+    S2_1 = zeros(1,m);
+    S2_1(1:5) = S2_U;
+    S2_m = zeros(1,m);
+    S2_m(m-4:m) = fliplr(S2_U);
+
+    S3_U = [-0.5e1/0.2e1 9 -12 7 -0.3e1/0.2e1;]/h^3;
+    S3_1 = zeros(1,m);
+    S3_1(1:5) = S3_U;
+    S3_m = zeros(1,m);
+    S3_m(m-4:m) = fliplr(-S3_U);
+
 
     x1=0.70127127127127;
 
@@ -39,21 +60,6 @@
     D1(m-5:m,m-8:m) = flipud( fliplr(-D1(1:6,1:9)));
     D1 = D1/h;
 
-    e_1 = zeros(m,1);e_1(1)=1;
-    e_m = zeros(m,1);e_m(m)=1;
-
-    S_U = [-25/12, 4, -3, 4/3, -1/4]/h;
-    S_1 = zeros(1,m);
-    S_1(1:5) = S_U;
-    S_m = zeros(1,m);
-    S_m(m-4:m) = fliplr(-S_U);
-
-
-    %DS=zeros(m,m);
-    %DS(1,1:5)=-[-25/12, 4, -3, 4/3, -1/4];
-    %DS(m,m-4:m)=fliplr(-[-25/12, 4, -3, 4/3, -1/4]);
-    %DS=diag(c)*DS/h;
-
 
     H = h*H;
     HI = inv(H);
@@ -97,15 +103,6 @@
 
     D2 = HI*(-M-diag(c)*e_1*S_1+diag(c)*e_m*S_m);
 
-    S2_U = [0.35e2/0.12e2 -0.26e2/0.3e1 0.19e2/0.2e1 -0.14e2/0.3e1 0.11e2/0.12e2;]/h^2;
-    S2_1 = zeros(1,m);
-    S2_1(1:5) = S2_U;
-    S2_m = zeros(1,m);
-    S2_m(m-4:m) = fliplr(S2_U);
-
-
-
-
 
     % Fourth derivative, 1th order accurate at first 8 boundary points (still
     % yield 5th order convergence if stable: for example u_tt=-u_xxxx
@@ -134,12 +131,6 @@
     M4(m-5:m,m-5:m) = flipud( fliplr( M4_U ) );
     M4 = M4/h^3;
 
-    S3_U = [-0.5e1/0.2e1 9 -12 7 -0.3e1/0.2e1;]/h^3;
-    S3_1 = zeros(1,m);
-    S3_1(1:5) = S3_U;
-    S3_m = zeros(1,m);
-    S3_m(m-4:m) = fliplr(-S3_U);
-
     D4 = HI*(M4-e_1*S3_1+e_m*S3_m  + S_1'*S2_1-S_m'*S2_m);
 
 end
--- a/+sbp/+implementations/d4_variable_8_higher_boundary_order.m	Fri Sep 23 22:58:45 2016 +0200
+++ b/+sbp/+implementations/d4_variable_8_higher_boundary_order.m	Fri Sep 23 23:10:44 2016 +0200
@@ -31,8 +31,10 @@
 
     H(m-7:m,m-7:m) = fliplr(flipud(H(1:8,1:8)));
 
-    e_1 = zeros(m,1);e_1(1) = 1;
-    e_m = zeros(m,1);e_m(m) = 1;
+    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);
@@ -40,6 +42,18 @@
     S_m = zeros(1,m);
     S_m(m-6:m) = fliplr(-S_U);
 
+    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);
+
+    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);
+
     H = h*H;
     HI = inv(H);
 
@@ -81,13 +95,6 @@
     %
     % D2 = HI*(-M-diag(c)*e_1*S_1+diag(c)*e_m*S_m);
 
-    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);
-
-
     % Fourth derivative, 1th order accurate at first 8 boundary points (still
     % yield 5th order convergence if stable: for example u_tt = -u_xxxx
 
@@ -118,11 +125,5 @@
     M4(m-7:m,m-7:m) = flipud( fliplr( M4_U ) );
     M4 = M4/h^3;
 
-    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);
-
     D4 = HI*(M4-e_1*S3_1+e_m*S3_m  + S_1'*S2_1-S_m'*S2_m);
 end
--- a/+sbp/+implementations/d4_variable_8_min_boundary_points.m	Fri Sep 23 22:58:45 2016 +0200
+++ b/+sbp/+implementations/d4_variable_8_min_boundary_points.m	Fri Sep 23 23:10:44 2016 +0200
@@ -32,6 +32,20 @@
     S_m = zeros(1,m);
     S_m(m-5:m) = fliplr(-S_U);
 
+    S2_U = [0.15e2/0.4e1 -0.77e2/0.6e1 0.107e3/0.6e1 -13 0.61e2/0.12e2 -0.5e1/0.6e1;]/h^2;
+    S2_1 = zeros(1,m);
+    S2_1(1:6) = S2_U;
+    S2_m = zeros(1,m);
+    S2_m(m-5:m) = fliplr(S2_U);
+
+    S3_U = [-0.17e2/0.4e1 0.71e2/0.4e1 -0.59e2/0.2e1 0.49e2/0.2e1 -0.41e2/0.4e1 0.7e1/0.4e1;]/h^3;
+    S3_1 = zeros(1,m);
+    S3_1(1:6) = S3_U;
+    S3_m = zeros(1,m);
+    S3_m(m-5:m) = fliplr(-S3_U);
+
+
+
     H = h*H;
     HI = inv(H);
 
@@ -74,12 +88,6 @@
     %
     % D2=HI*(-M-diag(c)*e_1*S_1+diag(c)*e_m*S_m);
 
-    S2_U = [0.15e2/0.4e1 -0.77e2/0.6e1 0.107e3/0.6e1 -13 0.61e2/0.12e2 -0.5e1/0.6e1;]/h^2;
-    S2_1 = zeros(1,m);
-    S2_1(1:6) = S2_U;
-    S2_m = zeros(1,m);
-    S2_m(m-5:m) = fliplr(S2_U);
-
 
 
     % Fourth derivative, 1th order accurate at first 8 boundary points (still
@@ -112,11 +120,5 @@
     M4(m-7:m,m-7:m) = flipud( fliplr( M4_U ) );
     M4 = M4/h^3;
 
-    S3_U = [-0.17e2/0.4e1 0.71e2/0.4e1 -0.59e2/0.2e1 0.49e2/0.2e1 -0.41e2/0.4e1 0.7e1/0.4e1;]/h^3;
-    S3_1 = zeros(1,m);
-    S3_1(1:6) = S3_U;
-    S3_m = zeros(1,m);
-    S3_m(m-5:m) = fliplr(-S3_U);
-
     D4=HI*(M4-e_1*S3_1+e_m*S3_m  + S_1'*S2_1-S_m'*S2_m);
 end