Mercurial > repos > public > sbplib

diff +sbp/+implementations/d4_compatible_2.m @ 267:f7ac3cd6eeaa operator_remake

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Sparsified all implementation files, removed all matlab warnings, fixed small bugs on minimum grid points.
author Martin Almquist <martin.almquist@it.uu.se>
date Fri, 09 Sep 2016 14:53:41 +0200
parents bfa130b7abf6
children
line wrap: on
line diff
--- a/+sbp/+implementations/d4_compatible_2.m	Fri Sep 09 11:03:13 2016 +0200
+++ b/+sbp/+implementations/d4_compatible_2.m	Fri Sep 09 14:53:41 2016 +0200
@@ -26,12 +26,12 @@
     % Vi b?rjar med normen. Notera att alla SBP operatorer delar samma norm,
     % vilket ?r n?dv?ndigt f?r stabilitet
     
-    BP = 1;
+    BP = 4;
     if(m<2*BP)
         error(['Operator requires at least ' num2str(2*BP) ' grid points']);
     end
 
-    H=diag(ones(m,1),0);H(1,1)=1/2;H(m,m)=1/2;
+    H=speye(m,m);H(1,1)=1/2;H(m,m)=1/2;
 
 
     H=H*h;
@@ -41,34 +41,42 @@
     % First derivative SBP operator, 1st order accurate at first 6 boundary points
 
     q1=1/2;
-    Q=q1*(diag(ones(m-1,1),1)-diag(ones(m-1,1),-1));
+%     Q=q1*(diag(ones(m-1,1),1)-diag(ones(m-1,1),-1));
+    stencil = [-q1,0,q1];
+    d = (length(stencil)-1)/2;
+    diags = -d:d;
+    Q = stripeMatrix(stencil, diags, m);
 
     %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_1=zeros(m,1);e_1(1)=1;
-    e_m=zeros(m,1);e_m(m)=1;
+    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') ;
+    D1=H\(Q-1/2*(e_1*e_1')+1/2*(e_m*e_m')) ;
 
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 
 
     % Second derivative, 1st order accurate at first 6 boundary points
-    m1=-1;m0=2;
-    M=m1*(diag(ones(m-1,1),1)+diag(ones(m-1,1),-1))+m0*diag(ones(m,1),0);M(1,1)=1;M(m,m)=1;
-    M=M/h;
+%     m1=-1;m0=2;
+% %     M=m1*(diag(ones(m-1,1),1)+diag(ones(m-1,1),-1))+m0*diag(ones(m,1),0);M(1,1)=1;M(m,m)=1;
+%     stencil = [m2,m1,m0,m1,m2];
+%     d = (length(stencil)-1)/2;
+%     diags = -d:d;
+%     M = stripeMatrix(stencil, diags, m);
+%     M=M/h;
 
     S_U=[-1 1]/h;
-    S_1=zeros(1,m);
+    S_1=sparse(1,m);
     S_1(1:2)=S_U;
-    S_m=zeros(1,m);
+    S_m=sparse(1,m);
 
     S_m(m-1:m)=fliplr(-S_U);
 
-    D2=HI*(-M-e_1*S_1+e_m*S_m);
+%     D2=H\(-M-e_1*S_1+e_m*S_m);
 
 
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -77,28 +85,32 @@
 
     % Third derivative, 1st order accurate at first 6 boundary points
 
-    q2=1/2;q1=-1;
-    Q3=q2*(diag(ones(m-2,1),2)-diag(ones(m-2,1),-2))+q1*(diag(ones(m-1,1),1)-diag(ones(m-1,1),-1));
+%     q2=1/2;q1=-1;
+% %     Q3=q2*(diag(ones(m-2,1),2)-diag(ones(m-2,1),-2))+q1*(diag(ones(m-1,1),1)-diag(ones(m-1,1),-1));
+%     stencil = [-q2,-q1,0,q1,q2];
+%     d = (length(stencil)-1)/2;
+%     diags = -d:d;
+%     Q3 = stripeMatrix(stencil, diags, m);
 
     %QQ3=(-1/8*diag(ones(m-3,1),3) + 1*diag(ones(m-2,1),2) - 13/8*diag(ones(m-1,1),1) +13/8*diag(ones(m-1,1),-1) -1*diag(ones(m-2,1),-2) + 1/8*diag(ones(m-3,1),-3));
 
 
-    Q3_U = [0 -0.2e1 / 0.5e1 0.3e1 / 0.10e2 0.1e1 / 0.10e2; 0.2e1 / 0.5e1 0 -0.7e1 / 0.10e2 0.3e1 / 0.10e2; -0.3e1 / 0.10e2 0.7e1 / 0.10e2 0 -0.9e1 / 0.10e2; -0.1e1 / 0.10e2 -0.3e1 / 0.10e2 0.9e1 / 0.10e2 0;];
-    Q3(1:4,1:4)=Q3_U;
-    Q3(m-3:m,m-3:m)=flipud( fliplr( -Q3_U ) );
-    Q3=Q3/h^2;
+%     Q3_U = [0 -0.2e1 / 0.5e1 0.3e1 / 0.10e2 0.1e1 / 0.10e2; 0.2e1 / 0.5e1 0 -0.7e1 / 0.10e2 0.3e1 / 0.10e2; -0.3e1 / 0.10e2 0.7e1 / 0.10e2 0 -0.9e1 / 0.10e2; -0.1e1 / 0.10e2 -0.3e1 / 0.10e2 0.9e1 / 0.10e2 0;];
+%     Q3(1:4,1:4)=Q3_U;
+%     Q3(m-3:m,m-3:m)=rot90(  -Q3_U ,2 );
+%     Q3=Q3/h^2;
 
 
 
     S2_U=[1 -2 1;]/h^2;
-    S2_1=zeros(1,m);
+    S2_1=sparse(1,m);
     S2_1(1:3)=S2_U;
-    S2_m=zeros(1,m);
+    S2_m=sparse(1,m);
     S2_m(m-2:m)=fliplr(S2_U);
 
 
 
-    D3=HI*(Q3 - e_1*S2_1 + e_m*S2_m +1/2*S_1'*S_1 -1/2*S_m'*S_m ) ;
+%     D3=H\(Q3 - e_1*S2_1 + e_m*S2_m +1/2*(S_1'*S_1) -1/2*(S_m'*S_m) ) ;
 
     % Fourth derivative, 0th order accurate at first 6 boundary points (still
     % yield 4th order convergence if stable: for example u_tt=-u_xxxx
@@ -112,16 +124,16 @@
 
     M4(1:4,1:4)=M4_U;
 
-    M4(m-3:m,m-3:m)=flipud( fliplr( M4_U ) );
+    M4(m-3:m,m-3:m)=rot90(  M4_U ,2 );
     M4=M4/h^3;
 
     S3_U=[-1 3 -3 1;]/h^3;
-    S3_1=zeros(1,m);
+    S3_1=sparse(1,m);
     S3_1(1:4)=S3_U;
-    S3_m=zeros(1,m);
+    S3_m=sparse(1,m);
     S3_m(m-3:m)=fliplr(-S3_U);
 
-    D4=HI*(M4-e_1*S3_1+e_m*S3_m  + S_1'*S2_1-S_m'*S2_m);
+    D4=H\(M4-e_1*S3_1+e_m*S3_m  + S_1'*S2_1-S_m'*S2_m);