comparison +scheme/LaplaceCurvilinear.m @ 566:9c98a0526afc feature/grids/laplace_refactor

Switch implementation of boundary and interface to SBP notation
author Jonatan Werpers <jonatan@werpers.com>
date Fri, 01 Sep 2017 10:43:43 +0200
parents f4b0d0e84305
children 33b962620e24
comparison
equal deleted inserted replaced
565:f4b0d0e84305 566:9c98a0526afc
188 s_w = sqrt((e_w'*x_v).^2 + (e_w'*y_v).^2); 188 s_w = sqrt((e_w'*x_v).^2 + (e_w'*y_v).^2);
189 s_e = sqrt((e_e'*x_v).^2 + (e_e'*y_v).^2); 189 s_e = sqrt((e_e'*x_v).^2 + (e_e'*y_v).^2);
190 s_s = sqrt((e_s'*x_u).^2 + (e_s'*y_u).^2); 190 s_s = sqrt((e_s'*x_u).^2 + (e_s'*y_u).^2);
191 s_n = sqrt((e_n'*x_u).^2 + (e_n'*y_u).^2); 191 s_n = sqrt((e_n'*x_u).^2 + (e_n'*y_u).^2);
192 192
193 obj.d_w = -1*(a11_w*obj.du_w' + a12_w*obj.dv_w')'; 193 obj.d_w = -1*(spdiag(1./s_w)*(a11_w*obj.du_w' + a12_w*obj.dv_w'))';
194 obj.d_e = (a11_e*obj.du_e' + a12_e*obj.dv_e')'; 194 obj.d_e = (spdiag(1./s_e)*(a11_e*obj.du_e' + a12_e*obj.dv_e'))';
195 obj.d_s = -1*(a22_s*obj.dv_s' + a12_s*obj.du_s')'; 195 obj.d_s = -1*(spdiag(1./s_s)*(a22_s*obj.dv_s' + a12_s*obj.du_s'))';
196 obj.d_n = (a22_n*obj.dv_n' + a12_n*obj.du_n')'; 196 obj.d_n = (spdiag(1./s_n)*(a22_n*obj.dv_n' + a12_n*obj.du_n'))';
197 197
198 obj.Dx = spdiag( y_v./J)*Du + spdiag(-y_u./J)*Dv; 198 obj.Dx = spdiag( y_v./J)*Du + spdiag(-y_u./J)*Dv;
199 obj.Dy = spdiag(-x_v./J)*Du + spdiag( x_u./J)*Dv; 199 obj.Dy = spdiag(-x_v./J)*Du + spdiag( x_u./J)*Dv;
200 200
201 %% Boundary inner products 201 %% Boundary inner products
231 % neighbour_boundary is a string specifying which boundary to interface to. 231 % neighbour_boundary is a string specifying which boundary to interface to.
232 function [closure, penalty] = boundary_condition(obj, boundary, type, parameter) 232 function [closure, penalty] = boundary_condition(obj, boundary, type, parameter)
233 default_arg('type','neumann'); 233 default_arg('type','neumann');
234 default_arg('parameter', []); 234 default_arg('parameter', []);
235 235
236 [e, d, s, gamm, halfnorm_inv , ~, ~, ~, scale_factor] = obj.get_boundary_ops(boundary); 236 [e, d, s, gamm, H_b, ~] = obj.get_boundary_ops(boundary);
237 switch type 237 switch type
238 % Dirichlet boundary condition 238 % Dirichlet boundary condition
239 case {'D','d','dirichlet'} 239 case {'D','d','dirichlet'}
240 tuning = 1.2; 240 tuning = 1.2;
241 % tuning = 20.2; 241 % tuning = 20.2;
242 [e, F, s, gamm, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t] = obj.get_boundary_ops(boundary); 242
243 243 b1 = gamm*obj.lambda./obj.a11.^2;
244 u = obj; 244 b2 = gamm*obj.lambda./obj.a22.^2;
245 245
246 b1 = gamm*u.lambda./u.a11.^2; 246 tau1 = tuning * spdiag(-1./b1 - 1./b2);
247 b2 = gamm*u.lambda./u.a22.^2; 247 tau2 = 1;
248 248
249 tau = -1./b1 - 1./b2; 249 tau = (tau1*e + tau2*d)*H_b;
250 tau = tuning * spdiag(tau); 250
251 sig1 = 1; 251 closure = obj.a*obj.Hi*tau*e';
252 252 penalty = -obj.a*obj.Hi*tau;
253 penalty_parameter_1 = halfnorm_inv_n*(tau + sig1*halfnorm_inv_t*F*e'*halfnorm_t)*e;
254
255 closure = obj.Ji*obj.a * penalty_parameter_1*e';
256 penalty = -obj.Ji*obj.a * penalty_parameter_1;
257 253
258 254
259 % Neumann boundary condition 255 % Neumann boundary condition
260 case {'N','n','neumann'} 256 case {'N','n','neumann'}
261 tau1 = -1; 257 tau1 = -1;
262 tau2 = 0; 258 tau2 = 0;
263 tau = obj.a*obj.Ji*(tau1*e + tau2*d); 259 tau = (tau1*e + tau2*d)*H_b;
264 260
265 closure = halfnorm_inv*tau*d'; 261 closure = obj.a*obj.Hi*tau*d';
266 penalty = -halfnorm_inv*tau; 262 penalty = -obj.a*obj.Hi*tau;
267 263
268 % Characteristic boundary condition
269 case {'characteristic', 'char', 'c'}
270 default_arg('parameter', 1);
271 beta = parameter;
272
273 tau = -obj.a * 1/beta*obj.Ji*e;
274
275 closure{1} = halfnorm_inv*tau*spdiag(scale_factor)*e';
276 closure{2} = halfnorm_inv*tau*beta*d';
277 penalty = -halfnorm_inv*tau;
278 264
279 % Unknown, boundary condition 265 % Unknown, boundary condition
280 otherwise 266 otherwise
281 error('No such boundary condition: type = %s',type); 267 error('No such boundary condition: type = %s',type);
282 end 268 end
285 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) 271 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
286 % u denotes the solution in the own domain 272 % u denotes the solution in the own domain
287 % v denotes the solution in the neighbour domain 273 % v denotes the solution in the neighbour domain
288 tuning = 1.2; 274 tuning = 1.2;
289 % tuning = 20.2; 275 % tuning = 20.2;
290 [e_u, F_u, s_u, gamm_u, halfnorm_inv_u_n, halfnorm_inv_u_t, halfnorm_u_t, I_u] = obj.get_boundary_ops(boundary); 276 [e_u, d_u, s_u, gamm_u, H_b_u, I_u] = obj.get_boundary_ops(boundary);
291 [e_v, F_v, s_v, gamm_v, halfnorm_inv_v_n, halfnorm_inv_v_t, halfnorm_v_t, I_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); 277 [e_v, d_v, s_v, gamm_v, H_b_v, I_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary);
292 278
293 u = obj; 279 u = obj;
294 v = neighbour_scheme; 280 v = neighbour_scheme;
295 281
296 b1_u = gamm_u*u.lambda(I_u)./u.a11(I_u).^2; 282 b1_u = gamm_u*u.lambda(I_u)./u.a11(I_u).^2;
297 b2_u = gamm_u*u.lambda(I_u)./u.a22(I_u).^2; 283 b2_u = gamm_u*u.lambda(I_u)./u.a22(I_u).^2;
298 b1_v = gamm_v*v.lambda(I_v)./v.a11(I_v).^2; 284 b1_v = gamm_v*v.lambda(I_v)./v.a11(I_v).^2;
299 b2_v = gamm_v*v.lambda(I_v)./v.a22(I_v).^2; 285 b2_v = gamm_v*v.lambda(I_v)./v.a22(I_v).^2;
300 286
301 tau = -1./(4*b1_u) -1./(4*b1_v) -1./(4*b2_u) -1./(4*b2_v); 287 tau1 = -1./(4*b1_u) -1./(4*b1_v) -1./(4*b2_u) -1./(4*b2_v);
302 tau = tuning * spdiag(tau); 288 tau1 = tuning * spdiag(tau1);
303 sig1 = 1/2; 289 tau2 = 1/2;
304 sig2 = -1/2; 290
305 291 sig1 = -1/2;
306 penalty_parameter_1 = halfnorm_inv_u_n*(e_u*tau + sig1*halfnorm_inv_u_t*F_u*e_u'*halfnorm_u_t*e_u); 292 sig2 = 0;
307 penalty_parameter_2 = halfnorm_inv_u_n * sig2 * e_u; 293
308 294 tau = (e_u*tau1 + tau2*d_u)*H_b_u;
309 295 sig = (sig1*e_u + sig2*d_u)*H_b_u;
310 closure = obj.Ji*obj.a * ( penalty_parameter_1*e_u' + penalty_parameter_2*F_u'); 296
311 penalty = obj.Ji*obj.a * (-penalty_parameter_1*e_v' + penalty_parameter_2*F_v'); 297 closure = obj.a*obj.Hi*( tau*e_u' + sig*d_u');
298 penalty = obj.a*obj.Hi*(-tau*e_v' + sig*d_v');
312 end 299 end
313 300
314 % Ruturns the boundary ops and sign for the boundary specified by the string boundary. 301 % Ruturns the boundary ops and sign for the boundary specified by the string boundary.
315 % The right boundary is considered the positive boundary 302 % The right boundary is considered the positive boundary
316 % 303 %
317 % I -- the indecies of the boundary points in the grid matrix 304 % I -- the indecies of the boundary points in the grid matrix
318 function [e, d, s, gamm, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t, I, scale_factor] = get_boundary_ops(obj, boundary) 305 function [e, d, s, gamm, H_b, I] = get_boundary_ops(obj, boundary)
319 306
320 % gridMatrix = zeros(obj.m(2),obj.m(1)); 307 % gridMatrix = zeros(obj.m(2),obj.m(1));
321 % gridMatrix(:) = 1:numel(gridMatrix); 308 % gridMatrix(:) = 1:numel(gridMatrix);
322 309
323 ind = grid.funcToMatrix(obj.grid, 1:prod(obj.m)); 310 ind = grid.funcToMatrix(obj.grid, 1:prod(obj.m));
324 311
325 switch boundary 312 switch boundary
326 case 'w' 313 case 'w'
327 e = obj.e_w; 314 e = obj.e_w;
328 d = obj.d_w; 315 d = obj.d_w;
316 H_b = obj.H_w;
329 s = -1; 317 s = -1;
330
331 I = ind(1,:); 318 I = ind(1,:);
332 scale_factor = sqrt(obj.x_v(I).^2 + obj.y_v(I).^2);
333 case 'e' 319 case 'e'
334 e = obj.e_e; 320 e = obj.e_e;
335 d = obj.d_e; 321 d = obj.d_e;
322 H_b = obj.H_e;
336 s = 1; 323 s = 1;
337
338 I = ind(end,:); 324 I = ind(end,:);
339 scale_factor = sqrt(obj.x_v(I).^2 + obj.y_v(I).^2);
340 case 's' 325 case 's'
341 e = obj.e_s; 326 e = obj.e_s;
342 d = obj.d_s; 327 d = obj.d_s;
328 H_b = obj.H_s;
343 s = -1; 329 s = -1;
344
345 I = ind(:,1)'; 330 I = ind(:,1)';
346 scale_factor = sqrt(obj.x_u(I).^2 + obj.y_u(I).^2);
347 case 'n' 331 case 'n'
348 e = obj.e_n; 332 e = obj.e_n;
349 d = obj.d_n; 333 d = obj.d_n;
334 H_b = obj.H_n;
350 s = 1; 335 s = 1;
351
352 I = ind(:,end)'; 336 I = ind(:,end)';
353 scale_factor = sqrt(obj.x_u(I).^2 + obj.y_u(I).^2);
354 otherwise 337 otherwise
355 error('No such boundary: boundary = %s',boundary); 338 error('No such boundary: boundary = %s',boundary);
356 end 339 end
357 340
358 switch boundary 341 switch boundary
359 case {'w','e'} 342 case {'w','e'}
360 halfnorm_inv_n = obj.Hiu;
361 halfnorm_inv_t = obj.Hiv;
362 halfnorm_t = obj.Hv;
363 gamm = obj.gamm_u; 343 gamm = obj.gamm_u;
364 case {'s','n'} 344 case {'s','n'}
365 halfnorm_inv_n = obj.Hiv;
366 halfnorm_inv_t = obj.Hiu;
367 halfnorm_t = obj.Hu;
368 gamm = obj.gamm_v; 345 gamm = obj.gamm_v;
369 end 346 end
370 end 347 end
371 348
372 function N = size(obj) 349 function N = size(obj)
373 N = prod(obj.m); 350 N = prod(obj.m);
374 end 351 end
375
376
377 end 352 end
378 end 353 end