Mercurial > repos > public > sbplib

comparison +sbp/+implementations/d4_compatible_4.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
comparison
equal deleted inserted replaced
266:bfa130b7abf6 267:f7ac3cd6eeaa
29 BP = 6; 29 BP = 6;
30 if(m<2*BP) 30 if(m<2*BP)
31 error(['Operator requires at least ' num2str(2*BP) ' grid points']); 31 error(['Operator requires at least ' num2str(2*BP) ' grid points']);
32 end 32 end
33 33
34 H=diag(ones(m,1),0); 34 H=speye(m,m);
35 H_U=[0.3e1 / 0.11e2 0 0 0 0 0; 0 0.2125516311e10 / 0.1311004640e10 0 0 0 0; 0 0 0.278735189e9 / 0.1966506960e10 0 0 0; 0 0 0 0.285925927e9 / 0.163875580e9 0 0; 0 0 0 0 0.1284335339e10 / 0.1966506960e10 0; 0 0 0 0 0 0.4194024163e10 / 0.3933013920e10;]; 35 H_U=[0.3e1 / 0.11e2 0 0 0 0 0; 0 0.2125516311e10 / 0.1311004640e10 0 0 0 0; 0 0 0.278735189e9 / 0.1966506960e10 0 0 0; 0 0 0 0.285925927e9 / 0.163875580e9 0 0; 0 0 0 0 0.1284335339e10 / 0.1966506960e10 0; 0 0 0 0 0 0.4194024163e10 / 0.3933013920e10;];
36 H(1:6,1:6)=H_U; 36 H(1:6,1:6)=H_U;
37 H(m-5:m,m-5:m)=fliplr(flipud(H_U)); 37 H(m-5:m,m-5:m)=rot90(H_U,2);
38 H=H*h; 38 H=H*h;
39 HI=inv(H); 39 HI=inv(H);
40 40
41 41
42 % First derivative SBP operator, 1st order accurate at first 6 boundary points 42 % First derivative SBP operator, 1st order accurate at first 6 boundary points
43 43
44 q2=-1/12;q1=8/12; 44 q2=-1/12;q1=8/12;
45 Q=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)); 45 % Q=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));
46 stencil = [-q2,-q1,0,q1,q2];
47 d = (length(stencil)-1)/2;
48 diags = -d:d;
49 Q = stripeMatrix(stencil, diags, m);
46 50
47 %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)); 51 %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));
48 52
49 Q_U = [0 0.9e1 / 0.11e2 -0.9e1 / 0.22e2 0.1e1 / 0.11e2 0 0; -0.9e1 / 0.11e2 0 0.2595224893e10 / 0.2622009280e10 -0.151435707e9 / 0.327751160e9 0.1112665611e10 / 0.2622009280e10 -0.1290899e7 / 0.9639740e7; 0.9e1 / 0.22e2 -0.2595224893e10 / 0.2622009280e10 0 0.1468436423e10 / 0.983253480e9 -0.1194603401e10 / 0.983253480e9 0.72033031e8 / 0.238364480e9; -0.1e1 / 0.11e2 0.151435707e9 / 0.327751160e9 -0.1468436423e10 / 0.983253480e9 0 0.439819541e9 / 0.327751160e9 -0.215942641e9 / 0.983253480e9; 0 -0.1112665611e10 / 0.2622009280e10 0.1194603401e10 / 0.983253480e9 -0.439819541e9 / 0.327751160e9 0 0.1664113643e10 / 0.2622009280e10; 0 0.1290899e7 / 0.9639740e7 -0.72033031e8 / 0.238364480e9 0.215942641e9 / 0.983253480e9 -0.1664113643e10 / 0.2622009280e10 0;]; 53 Q_U = [0 0.9e1 / 0.11e2 -0.9e1 / 0.22e2 0.1e1 / 0.11e2 0 0; -0.9e1 / 0.11e2 0 0.2595224893e10 / 0.2622009280e10 -0.151435707e9 / 0.327751160e9 0.1112665611e10 / 0.2622009280e10 -0.1290899e7 / 0.9639740e7; 0.9e1 / 0.22e2 -0.2595224893e10 / 0.2622009280e10 0 0.1468436423e10 / 0.983253480e9 -0.1194603401e10 / 0.983253480e9 0.72033031e8 / 0.238364480e9; -0.1e1 / 0.11e2 0.151435707e9 / 0.327751160e9 -0.1468436423e10 / 0.983253480e9 0 0.439819541e9 / 0.327751160e9 -0.215942641e9 / 0.983253480e9; 0 -0.1112665611e10 / 0.2622009280e10 0.1194603401e10 / 0.983253480e9 -0.439819541e9 / 0.327751160e9 0 0.1664113643e10 / 0.2622009280e10; 0 0.1290899e7 / 0.9639740e7 -0.72033031e8 / 0.238364480e9 0.215942641e9 / 0.983253480e9 -0.1664113643e10 / 0.2622009280e10 0;];
50 Q(1:6,1:6)=Q_U; 54 Q(1:6,1:6)=Q_U;
51 Q(m-5:m,m-5:m)=flipud( fliplr( -Q_U ) ); 55 Q(m-5:m,m-5:m)=rot90( -Q_U ,2 );
52 56
53 e_1=zeros(m,1);e_1(1)=1; 57 e_1=sparse(m,1);e_1(1)=1;
54 e_m=zeros(m,1);e_m(m)=1; 58 e_m=sparse(m,1);e_m(m)=1;
55 59
56 60
57 D1=HI*(Q-1/2*e_1*e_1'+1/2*e_m*e_m') ; 61 D1=H\(Q-1/2*(e_1*e_1')+1/2*(e_m*e_m')) ;
58 62
59 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 63 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
60 64
61 65
62 66
70 % 74 %
71 % M(m-5:m,m-5:m)=flipud( fliplr( M_U ) ); 75 % M(m-5:m,m-5:m)=flipud( fliplr( M_U ) );
72 % M=M/h; 76 % M=M/h;
73 % 77 %
74 S_U=[-0.11e2 / 0.6e1 3 -0.3e1 / 0.2e1 0.1e1 / 0.3e1;]/h; 78 S_U=[-0.11e2 / 0.6e1 3 -0.3e1 / 0.2e1 0.1e1 / 0.3e1;]/h;
75 S_1=zeros(1,m); 79 S_1=sparse(1,m);
76 S_1(1:4)=S_U; 80 S_1(1:4)=S_U;
77 S_m=zeros(1,m); 81 S_m=sparse(1,m);
78 82
79 S_m(m-3:m)=fliplr(-S_U); 83 S_m(m-3:m)=fliplr(-S_U);
80 84
81 % D2=HI*(-M-e_1*S_1+e_m*S_m); 85 % D2=H\(-M-e_1*S_1+e_m*S_m);
82 86
83 87
84 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 88 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
85 89
86 90
100 % Q3=Q3/h^2; 104 % Q3=Q3/h^2;
101 105
102 106
103 107
104 S2_U=[2 -5 4 -1;]/h^2; 108 S2_U=[2 -5 4 -1;]/h^2;
105 S2_1=zeros(1,m); 109 S2_1=sparse(1,m);
106 S2_1(1:4)=S2_U; 110 S2_1(1:4)=S2_U;
107 S2_m=zeros(1,m); 111 S2_m=sparse(1,m);
108 S2_m(m-3:m)=fliplr(S2_U); 112 S2_m(m-3:m)=fliplr(S2_U);
109 113
110 114
111 115
112 %D3=HI*(Q3 - e_1*S2_1 + e_m*S2_m +1/2*S_1'*S_1 -1/2*S_m'*S_m ) ; 116 %D3=H\(Q3 - e_1*S2_1 + e_m*S2_m +1/2*(S_1'*S_1) -1/2*(S_m'*S_m) ) ;
113 117
114 % Fourth derivative, 0th order accurate at first 6 boundary points (still 118 % Fourth derivative, 0th order accurate at first 6 boundary points (still
115 % yield 4th order convergence if stable: for example u_tt=-u_xxxx 119 % yield 4th order convergence if stable: for example u_tt=-u_xxxx
116 120
117 m3=-1/6;m2=2;m1=-13/2;m0=28/3; 121 m3=-1/6;m2=2;m1=-13/2;m0=28/3;
118 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); 122 % 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);
123
124 stencil = [m3,m2,m1,m0,m1,m2,m3];
125 d = (length(stencil)-1)/2;
126 diags = -d:d;
127 M4 = stripeMatrix(stencil, diags, m);
119 128
120 %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)); 129 %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));
121 130
122 M4_U=[0.227176919517319e15 / 0.94899692875680e14 -0.15262605263734e14 / 0.2965615402365e13 0.20205404771243e14 / 0.6778549491120e13 -0.3998303664097e13 / 0.23724923218920e14 0.1088305091927e13 / 0.94899692875680e14 -0.1686077077693e13 / 0.23724923218920e14; -0.15262605263734e14 / 0.2965615402365e13 0.280494781164181e15 / 0.23724923218920e14 -0.46417445546261e14 / 0.5931230804730e13 0.1705307929429e13 / 0.1694637372780e13 -0.553547394061e12 / 0.5931230804730e13 0.5615721694973e13 / 0.23724923218920e14; 0.20205404771243e14 / 0.6778549491120e13 -0.46417445546261e14 / 0.5931230804730e13 0.4135802350237e13 / 0.551742400440e12 -0.4140981465247e13 / 0.1078405600860e13 0.75538453067437e14 / 0.47449846437840e14 -0.4778134936391e13 / 0.11862461609460e14; -0.3998303664097e13 / 0.23724923218920e14 0.1705307929429e13 / 0.1694637372780e13 -0.4140981465247e13 / 0.1078405600860e13 0.20760974175677e14 / 0.2965615402365e13 -0.138330689701889e15 / 0.23724923218920e14 0.23711317526909e14 / 0.11862461609460e14; 0.1088305091927e13 / 0.94899692875680e14 -0.553547394061e12 / 0.5931230804730e13 0.75538453067437e14 / 0.47449846437840e14 -0.138330689701889e15 / 0.23724923218920e14 0.120223780251937e15 / 0.13557098982240e14 -0.151383731537477e15 / 0.23724923218920e14; -0.1686077077693e13 / 0.23724923218920e14 0.5615721694973e13 / 0.23724923218920e14 -0.4778134936391e13 / 0.11862461609460e14 0.23711317526909e14 / 0.11862461609460e14 -0.151383731537477e15 / 0.23724923218920e14 0.220304030094121e15 / 0.23724923218920e14;]; 131 M4_U=[0.227176919517319e15 / 0.94899692875680e14 -0.15262605263734e14 / 0.2965615402365e13 0.20205404771243e14 / 0.6778549491120e13 -0.3998303664097e13 / 0.23724923218920e14 0.1088305091927e13 / 0.94899692875680e14 -0.1686077077693e13 / 0.23724923218920e14; -0.15262605263734e14 / 0.2965615402365e13 0.280494781164181e15 / 0.23724923218920e14 -0.46417445546261e14 / 0.5931230804730e13 0.1705307929429e13 / 0.1694637372780e13 -0.553547394061e12 / 0.5931230804730e13 0.5615721694973e13 / 0.23724923218920e14; 0.20205404771243e14 / 0.6778549491120e13 -0.46417445546261e14 / 0.5931230804730e13 0.4135802350237e13 / 0.551742400440e12 -0.4140981465247e13 / 0.1078405600860e13 0.75538453067437e14 / 0.47449846437840e14 -0.4778134936391e13 / 0.11862461609460e14; -0.3998303664097e13 / 0.23724923218920e14 0.1705307929429e13 / 0.1694637372780e13 -0.4140981465247e13 / 0.1078405600860e13 0.20760974175677e14 / 0.2965615402365e13 -0.138330689701889e15 / 0.23724923218920e14 0.23711317526909e14 / 0.11862461609460e14; 0.1088305091927e13 / 0.94899692875680e14 -0.553547394061e12 / 0.5931230804730e13 0.75538453067437e14 / 0.47449846437840e14 -0.138330689701889e15 / 0.23724923218920e14 0.120223780251937e15 / 0.13557098982240e14 -0.151383731537477e15 / 0.23724923218920e14; -0.1686077077693e13 / 0.23724923218920e14 0.5615721694973e13 / 0.23724923218920e14 -0.4778134936391e13 / 0.11862461609460e14 0.23711317526909e14 / 0.11862461609460e14 -0.151383731537477e15 / 0.23724923218920e14 0.220304030094121e15 / 0.23724923218920e14;];
123 132
124 M4(1:6,1:6)=M4_U; 133 M4(1:6,1:6)=M4_U;
125 134
126 M4(m-5:m,m-5:m)=flipud( fliplr( M4_U ) ); 135 M4(m-5:m,m-5:m)=rot90( M4_U ,2 );
127 M4=M4/h^3; 136 M4=M4/h^3;
128 137
129 S3_U=[-1 3 -3 1;]/h^3; 138 S3_U=[-1 3 -3 1;]/h^3;
130 S3_1=zeros(1,m); 139 S3_1=sparse(1,m);
131 S3_1(1:4)=S3_U; 140 S3_1(1:4)=S3_U;
132 S3_m=zeros(1,m); 141 S3_m=sparse(1,m);
133 S3_m(m-3:m)=fliplr(-S3_U); 142 S3_m(m-3:m)=fliplr(-S3_U);
134 143
135 D4=HI*(M4-e_1*S3_1+e_m*S3_m + S_1'*S2_1-S_m'*S2_m); 144 D4=H\(M4-e_1*S3_1+e_m*S3_m + S_1'*S2_1-S_m'*S2_m);
136 145
137 S_1 = S_1'; 146 S_1 = S_1';
138 S_m = S_m'; 147 S_m = S_m';
139 S2_1 = S2_1'; 148 S2_1 = S2_1';
140 S2_m = S2_m'; 149 S2_m = S2_m';