Mercurial > repos > public > sbplib

comparison +scheme/LaplaceCurvilinear.m @ 452:56eb7c088bd4 feature/grids

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Allow any coefficient in LaplaceCurvilinear
author Jonatan Werpers <jonatan@werpers.com>
date Wed, 10 May 2017 13:16:12 +0200
parents 8d455e49364f
children 08b6281ba2a9
comparison
equal deleted inserted replaced
451:4e266dfe9edc 452:56eb7c088bd4
7 7
8 order % Order accuracy for the approximation 8 order % Order accuracy for the approximation
9 9
10 D % non-stabalized scheme operator 10 D % non-stabalized scheme operator
11 M % Derivative norm 11 M % Derivative norm
12 c 12 a,b
13 J, Ji 13 J, Ji
14 a11, a12, a22 14 a11, a12, a22
15 15
16 H % Discrete norm 16 H % Discrete norm
17 Hi 17 Hi
34 y_u 34 y_u
35 y_v 35 y_v
36 end 36 end
37 37
38 methods 38 methods
39 function obj = LaplaceCurvilinear(g ,order, c, opSet) 39 % Implements a*div(b*grad(u)) as a SBP scheme
40 function obj = LaplaceCurvilinear(g ,order, a, b, opSet)
40 default_arg('opSet',@sbp.D2Variable); 41 default_arg('opSet',@sbp.D2Variable);
41 default_arg('c', 1); 42 default_arg('a', 1);
43 default_arg('b', 1);
44
45
46 if b ~=1
47 error('Not implemented yet')
48 end
42 49
43 assert(isa(g, 'grid.Curvilinear')) 50 assert(isa(g, 'grid.Curvilinear'))
44 51
45 m = g.size(); 52 m = g.size();
46 m_u = m(1); 53 m_u = m(1);
144 obj.m = m; 151 obj.m = m;
145 obj.h = [h_u h_v]; 152 obj.h = [h_u h_v];
146 obj.order = order; 153 obj.order = order;
147 obj.grid = g; 154 obj.grid = g;
148 155
149 obj.c = c; 156 obj.a = a;
157 obj.b = b;
150 obj.J = spdiags(J, 0, m_tot, m_tot); 158 obj.J = spdiags(J, 0, m_tot, m_tot);
151 obj.Ji = spdiags(1./J, 0, m_tot, m_tot); 159 obj.Ji = spdiags(1./J, 0, m_tot, m_tot);
152 obj.a11 = a11; 160 obj.a11 = a11;
153 obj.a12 = a12; 161 obj.a12 = a12;
154 obj.a22 = a22; 162 obj.a22 = a22;
155 obj.D = obj.Ji*c^2*(Duu + Duv + Dvu + Dvv); 163 obj.D = obj.Ji*a*(Duu + Duv + Dvu + Dvv);
156 obj.lambda = lambda; 164 obj.lambda = lambda;
157 165
158 obj.Dx = spdiag( y_v./J)*Du + spdiag(-y_u./J)*Dv; 166 obj.Dx = spdiag( y_v./J)*Du + spdiag(-y_u./J)*Dv;
159 obj.Dy = spdiag(-x_v./J)*Du + spdiag( x_u./J)*Dv; 167 obj.Dy = spdiag(-x_v./J)*Du + spdiag( x_u./J)*Dv;
160 168
197 tau = tuning * spdiag(tau); 205 tau = tuning * spdiag(tau);
198 sig1 = 1; 206 sig1 = 1;
199 207
200 penalty_parameter_1 = halfnorm_inv_n*(tau + sig1*halfnorm_inv_t*F*e'*halfnorm_t)*e; 208 penalty_parameter_1 = halfnorm_inv_n*(tau + sig1*halfnorm_inv_t*F*e'*halfnorm_t)*e;
201 209
202 closure = obj.Ji*obj.c^2 * penalty_parameter_1*e'; 210 closure = obj.Ji*obj.a * penalty_parameter_1*e';
203 penalty = -obj.Ji*obj.c^2 * penalty_parameter_1; 211 penalty = -obj.Ji*obj.a * penalty_parameter_1;
204 212
205 213
206 % Neumann boundary condition 214 % Neumann boundary condition
207 case {'N','n','neumann'} 215 case {'N','n','neumann'}
208 c = obj.c;
209
210 a_n = spdiags(coeff_n,0,length(coeff_n),length(coeff_n)); 216 a_n = spdiags(coeff_n,0,length(coeff_n),length(coeff_n));
211 a_t = spdiags(coeff_t,0,length(coeff_t),length(coeff_t)); 217 a_t = spdiags(coeff_t,0,length(coeff_t),length(coeff_t));
212 d = (a_n * d_n' + a_t*d_t')'; 218 d = (a_n * d_n' + a_t*d_t')';
213 219
214 tau1 = -s; 220 tau1 = -s;
215 tau2 = 0; 221 tau2 = 0;
216 tau = c.^2 * obj.Ji*(tau1*e + tau2*d); 222 tau = obj.a * obj.Ji*(tau1*e + tau2*d);
217 223
218 closure = halfnorm_inv*tau*d'; 224 closure = halfnorm_inv*tau*d';
219 penalty = -halfnorm_inv*tau; 225 penalty = -halfnorm_inv*tau;
220 226
221 % Characteristic boundary condition 227 % Characteristic boundary condition
222 case {'characteristic', 'char', 'c'} 228 case {'characteristic', 'char', 'c'}
223 default_arg('parameter', 1); 229 default_arg('parameter', 1);
224 beta = parameter; 230 beta = parameter;
225 c = obj.c;
226 231
227 a_n = spdiags(coeff_n,0,length(coeff_n),length(coeff_n)); 232 a_n = spdiags(coeff_n,0,length(coeff_n),length(coeff_n));
228 a_t = spdiags(coeff_t,0,length(coeff_t),length(coeff_t)); 233 a_t = spdiags(coeff_t,0,length(coeff_t),length(coeff_t));
229 d = s*(a_n * d_n' + a_t*d_t')'; % outward facing normal derivative 234 d = s*(a_n * d_n' + a_t*d_t')'; % outward facing normal derivative
230 235
231 tau = -c.^2 * 1/beta*obj.Ji*e; 236 tau = -obj.a * 1/beta*obj.Ji*e;
232 237
233 closure{1} = halfnorm_inv*tau/c*spdiag(scale_factor)*e'; 238 closure{1} = halfnorm_inv*tau*spdiag(scale_factor)*e';
234 closure{2} = halfnorm_inv*tau*beta*d'; 239 closure{2} = halfnorm_inv*tau*beta*d';
235 penalty = -halfnorm_inv*tau; 240 penalty = -halfnorm_inv*tau;
236 241
237 % Unknown, boundary condition 242 % Unknown, boundary condition
238 otherwise 243 otherwise
271 276
272 penalty_parameter_1 = halfnorm_inv_u_n*(e_u*tau + sig1*halfnorm_inv_u_t*F_u*e_u'*halfnorm_u_t*e_u); 277 penalty_parameter_1 = halfnorm_inv_u_n*(e_u*tau + sig1*halfnorm_inv_u_t*F_u*e_u'*halfnorm_u_t*e_u);
273 penalty_parameter_2 = halfnorm_inv_u_n * sig2 * e_u; 278 penalty_parameter_2 = halfnorm_inv_u_n * sig2 * e_u;
274 279
275 280
276 closure = obj.Ji*obj.c^2 * ( penalty_parameter_1*e_u' + penalty_parameter_2*F_u'); 281 closure = obj.Ji*obj.a * ( penalty_parameter_1*e_u' + penalty_parameter_2*F_u');
277 penalty = obj.Ji*obj.c^2 * (-penalty_parameter_1*e_v' + penalty_parameter_2*F_v'); 282 penalty = obj.Ji*obj.a * (-penalty_parameter_1*e_v' + penalty_parameter_2*F_v');
278 end 283 end
279 284
280 % Ruturns the boundary ops and sign for the boundary specified by the string boundary. 285 % Ruturns the boundary ops and sign for the boundary specified by the string boundary.
281 % The right boundary is considered the positive boundary 286 % The right boundary is considered the positive boundary
282 % 287 %