Mercurial > repos > public > sbplib

comparison +scheme/LaplaceCurvilinear.m @ 904:14b093a344eb feature/utux2D

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Outline new interface method in LaplaceCurvilinear
author Martin Almquist <malmquist@stanford.edu>
date Thu, 22 Nov 2018 22:03:06 -0800
parents f4595f14d696
children 0499239496cf
comparison
equal deleted inserted replaced
903:703183ed8c8b 904:14b093a344eb
274 otherwise 274 otherwise
275 error('No such boundary condition: type = %s',type); 275 error('No such boundary condition: type = %s',type);
276 end 276 end
277 end 277 end
278 278
279 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) 279 % type Struct that specifies the interface coupling.
280 % Fields:
281 % -- tuning: penalty strength, defaults to 1.2
282 % -- interpolation: struct of interpolation operators (empty for conforming grids)
283 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type)
284 if isempty(type)
285 [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary);
286 else
287 assertType(type, 'struct');
288 if isfield(type, 'I_local2neighbor')
289 [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type);
290 else
291 [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type);
292 end
293 end
294 end
295
296 function [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type)
297 % u denotes the solution in the own domain
298 % v denotes the solution in the neighbour domain
299
300 default_arg('type', struct);
301 default_field(type, 'tuning', 1.2);
302 tuning = type.tuning;
303
304 [e_u, d_u, gamm_u, H_b_u, I_u] = obj.get_boundary_ops(boundary);
305 [e_v, d_v, gamm_v, H_b_v, I_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary);
306
307 u = obj;
308 v = neighbour_scheme;
309
310 b1_u = gamm_u*u.lambda(I_u)./u.a11(I_u).^2;
311 b2_u = gamm_u*u.lambda(I_u)./u.a22(I_u).^2;
312 b1_v = gamm_v*v.lambda(I_v)./v.a11(I_v).^2;
313 b2_v = gamm_v*v.lambda(I_v)./v.a22(I_v).^2;
314
315 tau1 = -1./(4*b1_u) -1./(4*b1_v) -1./(4*b2_u) -1./(4*b2_v);
316 tau1 = tuning * spdiag(tau1);
317 tau2 = 1/2;
318
319 sig1 = -1/2;
320 sig2 = 0;
321
322 tau = (e_u*tau1 + tau2*d_u)*H_b_u;
323 sig = (sig1*e_u + sig2*d_u)*H_b_u;
324
325 closure = obj.a*obj.Hi*( tau*e_u' + sig*d_u');
326 penalty = obj.a*obj.Hi*(-tau*e_v' + sig*d_v');
327 end
328
329 function [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type)
330
331 default_field(type, 'tuning', 1.2);
332 tuning = type.tuning;
333
334 I_u2v_good = type.I_local2neighbor.good;
335 I_u2v_bad = type.I_local2neighbor.bad;
336 I_v2u_good = type.I_neighbor2local.good;
337 I_v2u_bad = type.I_neighbor2local.bad;
280 338
281 [e_u, d_u, gamm_u, H_b_u, I_u] = obj.get_boundary_ops(boundary); 339 [e_u, d_u, gamm_u, H_b_u, I_u] = obj.get_boundary_ops(boundary);
282 [e_v, d_v, gamm_v, H_b_v, I_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); 340 [e_v, d_v, gamm_v, H_b_v, I_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary);
283 Hi = obj.Hi; 341 Hi = obj.Hi;
284 a = obj.a; 342 a = obj.a;
285
286 m_u = length(H_b_u);
287 m_v = length(H_b_v);
288
289 grid_ratio = m_u/m_v;
290 if grid_ratio ~= 1
291
292 [ms, index] = sort([m_u, m_v]);
293 orders = [obj.order, neighbour_scheme.order];
294 orders = orders(index);
295
296 switch obj.interpolation_type
297 case 'MC'
298 interpOpSet = sbp.InterpMC(ms(1),ms(2),orders(1),orders(2));
299 if grid_ratio < 1
300 I_v2u_good = interpOpSet.IF2C;
301 I_v2u_bad = interpOpSet.IF2C;
302 I_u2v_good = interpOpSet.IC2F;
303 I_u2v_bad = interpOpSet.IC2F;
304 elseif grid_ratio > 1
305 I_v2u_good = interpOpSet.IC2F;
306 I_v2u_bad = interpOpSet.IC2F;
307 I_u2v_good = interpOpSet.IF2C;
308 I_u2v_bad = interpOpSet.IF2C;
309 end
310 case 'AWW'
311 %String 'C2F' indicates that ICF2 is more accurate.
312 interpOpSetF2C = sbp.InterpAWW(ms(1),ms(2),orders(1),orders(2),'F2C');
313 interpOpSetC2F = sbp.InterpAWW(ms(1),ms(2),orders(1),orders(2),'C2F');
314 if grid_ratio < 1
315 % Local is coarser than neighbour
316 I_v2u_good = interpOpSetF2C.IF2C;
317 I_v2u_bad = interpOpSetC2F.IF2C;
318 I_u2v_good = interpOpSetC2F.IC2F;
319 I_u2v_bad = interpOpSetF2C.IC2F;
320 elseif grid_ratio > 1
321 % Local is finer than neighbour
322 I_v2u_good = interpOpSetC2F.IC2F;
323 I_v2u_bad = interpOpSetF2C.IC2F;
324 I_u2v_good = interpOpSetF2C.IF2C;
325 I_u2v_bad = interpOpSetC2F.IF2C;
326 end
327 otherwise
328 error(['Interpolation type ' obj.interpolation_type ...
329 ' is not available.' ]);
330 end
331
332 else
333 % No interpolation required
334 I_v2u_good = speye(m_u,m_u);
335 I_v2u_bad = speye(m_u,m_u);
336 I_u2v_good = speye(m_u,m_u);
337 I_u2v_bad = speye(m_u,m_u);
338 end
339 343
340 % u denotes the solution in the own domain 344 % u denotes the solution in the own domain
341 % v denotes the solution in the neighbour domain 345 % v denotes the solution in the neighbour domain
342 tuning = 1.2; 346 tuning = 1.2;
343 % tuning = 20.2; 347 % tuning = 20.2;
363 367
364 penalty = -a*Hi*e_u*tau_u*H_b_u*I_v2u_good*e_v' + ... 368 penalty = -a*Hi*e_u*tau_u*H_b_u*I_v2u_good*e_v' + ...
365 -a*Hi*e_u*H_b_u*I_v2u_bad*beta_u*e_v' + ... 369 -a*Hi*e_u*H_b_u*I_v2u_bad*beta_u*e_v' + ...
366 -a*1/2*Hi*d_u*H_b_u*I_v2u_good*e_v' + ... 370 -a*1/2*Hi*d_u*H_b_u*I_v2u_good*e_v' + ...
367 -a*1/2*Hi*e_u*H_b_u*I_v2u_bad*d_v'; 371 -a*1/2*Hi*e_u*H_b_u*I_v2u_bad*d_v';
368 372
369 373
370 end 374 end
371 375
372 % Ruturns the boundary ops and sign for the boundary specified by the string boundary. 376 % Ruturns the boundary ops and sign for the boundary specified by the string boundary.
373 % The right boundary is considered the positive boundary 377 % The right boundary is considered the positive boundary