Mercurial > repos > public > sbplib

comparison +sbp/+implementations/d4_variable_2.m @ 313:52b4cdf27633 feature/beams

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Cleaning order 2.
author Jonatan Werpers <jonatan@werpers.com>
date Fri, 23 Sep 2016 21:58:52 +0200
parents 9230c056a574
children 88584b0cfba1
comparison
equal deleted inserted replaced
312:9230c056a574 313:52b4cdf27633
8 %%% %%% 8 %%% %%%
9 %%% Datum: 2013-11-11 %%% 9 %%% Datum: 2013-11-11 %%%
10 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 10 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
11 11
12 BP = 4; 12 BP = 4;
13 if(m<2*BP) 13 if(m < 2*BP)
14 error(['Operator requires at least ' num2str(2*BP) ' grid points']); 14 error('Operator requires at least %d grid points', 2*BP);
15 end 15 end
16 16
17 % Norm
18 Hv = ones(m,1);
19 Hv(1) = 1/2;
20 Hv(m) = 1/2;
21 Hv = h*Hv;
22 H = spdiag(Hv, 0);
23 HI = spdiag(1./Hv, 0);
17 24
18 H=speye(m,m); 25 % Boundary operators
19 H(1,1)=1/2; 26 e_l = sparse(m,1);
20 H(m,m)=1/2; 27 e_l(1) = 1;
28 e_r = rot90(e_l, 2);
21 29
22 H=H*h; 30 d1_l = sparse(m,1);
23 HI=inv(H); 31 d1_l(1:3) = 1/h*[-3/2 2 -1/2];
32 d1_r = -rot90(d1_l);
33
34 d2_l = sparse(m,1);
35 d2_l(1:3) = 1/h^2*[1 -2 1];
36 d2_r = rot90(d2_l, 2);
37
38 d3_l = sparse(m,1);
39 d3_l(1:4) = 1/h^3*[-1 3 -3 1];
40 d3_r = -rot90(d3_l, 2)
24 41
25 42
26 % First derivative SBP operator, 1st order accurate at first 6 boundary points 43 % First derivative SBP operator, 1st order accurate at first 6 boundary points
27 44 stencil = [-1/2, 0, 1/2];
28 q1=1/2; 45 diags = [-1 0 1];
29 % Q=q1*(diag(ones(m-1,1),1)-diag(ones(m-1,1),-1));
30 stencil = [-q1,0,q1];
31 d = (length(stencil)-1)/2;
32 diags = -d:d;
33 Q = stripeMatrix(stencil, diags, m); 46 Q = stripeMatrix(stencil, diags, m);
34 47
35 %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 D1 = HI*(Q-1/2*(e_1*e_1') + 1/2*(e_m*e_m'));
36
37 e_1=sparse(m,1);
38 e_1(1)=1;
39 e_m=sparse(m,1);
40 e_m(m)=1;
41
42 D1=HI*(Q-1/2*(e_1*e_1')+1/2*(e_m*e_m')) ;
43
44 49
45 % Second derivative, 1st order accurate at first boundary points 50 % Second derivative, 1st order accurate at first boundary points
46 51 M = sparse(m,m);
47 % below for constant coefficients
48 % m1=-1;m0=2;
49 % 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;
50 % M=M/h;
51 %D2=HI*(-M-e_1*S_1+e_m*S_m);
52
53 % Below for variable coefficients
54 % Require a vector c with the koeffients
55
56 S_U=[-3/2 2 -1/2]/h;
57 S_1=sparse(1,m);
58 S_1(1:3)=S_U;
59 S_m=sparse(1,m);
60 S_m(m-2:m)=fliplr(-S_U);
61
62 S_1 = S_1';
63 S_m = S_m';
64
65 M=sparse(m,m);
66 e_1 = sparse(e_1);
67 e_m = sparse(e_m);
68 S_1 = sparse(S_1);
69 S_m = sparse(S_m);
70 52
71 scheme_width = 3; 53 scheme_width = 3;
72 scheme_radius = (scheme_width-1)/2; 54 scheme_radius = (scheme_width-1)/2;
73 r = (1+scheme_radius):(m-scheme_radius); 55 r = (1+scheme_radius):(m-scheme_radius);
74 56
75 function D2 = D2_fun(c) 57 function D2 = D2_fun(c)
76
77 Mm1 = -c(r-1)/2 - c(r)/2; 58 Mm1 = -c(r-1)/2 - c(r)/2;
78 M0 = c(r-1)/2 + c(r) + c(r+1)/2; 59 M0 = c(r-1)/2 + c(r) + c(r+1)/2;
79 Mp1 = -c(r)/2 - c(r+1)/2; 60 Mp1 = -c(r)/2 - c(r+1)/2;
80 61
81 M(r,:) = spdiags([Mm1 M0 Mp1],0:2*scheme_radius,length(r),m); 62 M(r,:) = spdiags([Mm1 M0 Mp1],0:2*scheme_radius,length(r),m);
82 63
64 M(1:2,1:2) = [c(1)/2 + c(2)/2 -c(1)/2 - c(2)/2; -c(1)/2 - c(2)/2 c(1)/2 + c(2) + c(3)/2;];
65 M(m-1:m,m-1:m) = [c(m-2)/2 + c(m-1) + c(m)/2 -c(m-1)/2 - c(m)/2; -c(m-1)/2 - c(m)/2 c(m-1)/2 + c(m)/2;];
66 M = 1/h*M;
83 67
84 M(1:2,1:2)=[c(1)/2 + c(2)/2 -c(1)/2 - c(2)/2; -c(1)/2 - c(2)/2 c(1)/2 + c(2) + c(3)/2;]; 68 D2 = HI*(-M - c(1)*e_1*d1_l' + c(m)*e_r*d1_r');
85 M(m-1:m,m-1:m)=[c(m-2)/2 + c(m-1) + c(m)/2 -c(m-1)/2 - c(m)/2; -c(m-1)/2 - c(m)/2 c(m-1)/2 + c(m)/2;];
86 M=M/h;
87
88 D2=HI*(-M-c(1)*e_1*S_1'+c(m)*e_m*S_m');
89 end 69 end
90 D2 = @D2_fun; 70 D2 = @D2_fun;
91 71
92 72 % Fourth derivative, 0th order accurate at first 6 boundary points
93 73 stencil = [1, -4, 6, -4, 1];
94 74 diags = -2:2;
95
96 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
97
98
99
100 % Third derivative, 1st order accurate at first 6 boundary points
101
102 q2=1/2;q1=-1;
103 % 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));
104 stencil = [-q2,-q1,0,q1,q2];
105 d = (length(stencil)-1)/2;
106 diags = -d:d;
107 Q3 = stripeMatrix(stencil, diags, m);
108
109 %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));
110
111
112 Q3_U = [
113 0 -0.13e2/0.16e2 0.7e1/0.8e1 -0.1e1/0.16e2;
114 0.13e2/0.16e2 0 -0.23e2/0.16e2 0.5e1/0.8e1;
115 -0.7e1/0.8e1 0.23e2/0.16e2 0 -0.17e2/0.16e2;
116 0.1e1/0.16e2 -0.5e1/0.8e1 0.17e2/0.16e2 0;
117 ];
118 Q3(1:4,1:4)=Q3_U;
119 Q3(m-3:m,m-3:m)=rot90( -Q3_U ,2 );
120 Q3=Q3/h^2;
121
122
123
124 S2_U=[1 -2 1;]/h^2;
125 S2_1=sparse(1,m);
126 S2_1(1:3)=S2_U;
127 S2_m=sparse(1,m);
128 S2_m(m-2:m)=fliplr(S2_U);
129 S2_1 = S2_1';
130 S2_m = S2_m';
131
132
133
134 D3=HI*(Q3 - e_1*S2_1' + e_m*S2_m' +1/2*(S_1*S_1') -1/2*(S_m*S_m') ) ;
135
136 % Fourth derivative, 0th order accurate at first 6 boundary points (still
137 % yield 4th order convergence if stable: for example u_tt=-u_xxxx
138
139 m2=1;m1=-4;m0=6;
140 % M4=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);
141 stencil = [m2,m1,m0,m1,m2];
142 d = (length(stencil)-1)/2;
143 diags = -d:d;
144 M4 = stripeMatrix(stencil, diags, m); 75 M4 = stripeMatrix(stencil, diags, m);
145 76
146 %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)); 77 M4_U = [
147 78 0.13e2/0.10e2 -0.12e2/0.5e1 0.9e1/0.10e2 0.1e1/0.5e1;
148 M4_U=[ 79 -0.12e2/0.5e1 0.26e2/0.5e1 -0.16e2/0.5e1 0.2e1/0.5e1;
149 0.13e2/0.10e2 -0.12e2/0.5e1 0.9e1/0.10e2 0.1e1/0.5e1; 80 0.9e1/0.10e2 -0.16e2/0.5e1 0.47e2/0.10e2 -0.17e2/0.5e1;
150 -0.12e2/0.5e1 0.26e2/0.5e1 -0.16e2/0.5e1 0.2e1/0.5e1; 81 0.1e1/0.5e1 0.2e1/0.5e1 -0.17e2/0.5e1 0.29e2/0.5e1;
151 0.9e1/0.10e2 -0.16e2/0.5e1 0.47e2/0.10e2 -0.17e2/0.5e1;
152 0.1e1/0.5e1 0.2e1/0.5e1 -0.17e2/0.5e1 0.29e2/0.5e1;
153 ]; 82 ];
154 83
155 M4(1:4,1:4)=M4_U; 84 M4(1:4,1:4) = M4_U;
156 85 M4(m-3:m,m-3:m) = rot90(M4_U, 2);
157 M4(m-3:m,m-3:m)=rot90( M4_U ,2 ); 86 M4 = 1/h^3*M4;
158 M4=M4/h^3;
159
160 S3_U=[-1 3 -3 1;]/h^3;
161 S3_1=sparse(1,m);
162 S3_1(1:4)=S3_U;
163 S3_m=sparse(1,m);
164 S3_m(m-3:m)=fliplr(-S3_U);
165 S3_1 = S3_1';
166 S3_m = S3_m';
167 87
168 D4=HI*(M4-e_1*S3_1'+e_m*S3_m' + S_1*S2_1'-S_m*S2_m'); 88 D4=HI*(M4-e_1*S3_1'+e_m*S3_m' + S_1*S2_1'-S_m*S2_m');
169 end 89 end