Mercurial > repos > public > sbplib

comparison +scheme/LaplaceCurvilinear.m @ 1033:037f203b9bf5 feature/burgers1d

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with branch feature/advectioRV to utilize the +rv package
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Thu, 17 Jan 2019 10:44:12 +0100
parents 706d1c2b4199
children a0b3161e44f3
comparison
equal deleted inserted replaced
854:18162a0a5bb5 1033:037f203b9bf5
36 du_e, dv_e 36 du_e, dv_e
37 du_s, dv_s 37 du_s, dv_s
38 du_n, dv_n 38 du_n, dv_n
39 gamm_u, gamm_v 39 gamm_u, gamm_v
40 lambda 40 lambda
41
41 end 42 end
42 43
43 methods 44 methods
44 % Implements a*div(b*grad(u)) as a SBP scheme 45 % Implements a*div(b*grad(u)) as a SBP scheme
45 % TODO: Implement proper H, it should be the real physical quadrature, the logic quadrature may be but in a separate variable (H_logic?) 46 % TODO: Implement proper H, it should be the real physical quadrature, the logic quadrature may be but in a separate variable (H_logic?)
51 52
52 if b ~=1 53 if b ~=1
53 error('Not implemented yet') 54 error('Not implemented yet')
54 end 55 end
55 56
56 assert(isa(g, 'grid.Curvilinear')) 57 % assert(isa(g, 'grid.Curvilinear'))
58 if isa(a, 'function_handle')
59 a = grid.evalOn(g, a);
60 a = spdiag(a);
61 end
57 62
58 m = g.size(); 63 m = g.size();
59 m_u = m(1); 64 m_u = m(1);
60 m_v = m(2); 65 m_v = m(2);
61 m_tot = g.N(); 66 m_tot = g.N();
266 otherwise 271 otherwise
267 error('No such boundary condition: type = %s',type); 272 error('No such boundary condition: type = %s',type);
268 end 273 end
269 end 274 end
270 275
271 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) 276 % type Struct that specifies the interface coupling.
277 % Fields:
278 % -- tuning: penalty strength, defaults to 1.2
279 % -- interpolation: type of interpolation, default 'none'
280 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type)
281
282 defaultType.tuning = 1.2;
283 defaultType.interpolation = 'none';
284 default_struct('type', defaultType);
285
286 switch type.interpolation
287 case {'none', ''}
288 [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type);
289 case {'op','OP'}
290 [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type);
291 otherwise
292 error('Unknown type of interpolation: %s ', type.interpolation);
293 end
294 end
295
296 function [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type)
297 tuning = type.tuning;
298
272 % u denotes the solution in the own domain 299 % u denotes the solution in the own domain
273 % v denotes the solution in the neighbour domain 300 % v denotes the solution in the neighbour domain
274 tuning = 1.2;
275 % tuning = 20.2;
276 [e_u, d_u, gamm_u, H_b_u, I_u] = obj.get_boundary_ops(boundary); 301 [e_u, d_u, gamm_u, H_b_u, I_u] = obj.get_boundary_ops(boundary);
277 [e_v, d_v, gamm_v, H_b_v, I_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); 302 [e_v, d_v, gamm_v, H_b_v, I_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary);
278 303
279 u = obj; 304 u = obj;
280 v = neighbour_scheme; 305 v = neighbour_scheme;
296 321
297 closure = obj.a*obj.Hi*( tau*e_u' + sig*d_u'); 322 closure = obj.a*obj.Hi*( tau*e_u' + sig*d_u');
298 penalty = obj.a*obj.Hi*(-tau*e_v' + sig*d_v'); 323 penalty = obj.a*obj.Hi*(-tau*e_v' + sig*d_v');
299 end 324 end
300 325
301 % Ruturns the boundary ops and sign for the boundary specified by the string boundary. 326 function [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type)
327
328 % TODO: Make this work for curvilinear grids
329 warning('LaplaceCurvilinear: Non-conforming grid interpolation only works for Cartesian grids.');
330
331 % User can request special interpolation operators by specifying type.interpOpSet
332 default_field(type, 'interpOpSet', @sbp.InterpOpsOP);
333 interpOpSet = type.interpOpSet;
334 tuning = type.tuning;
335
336
337 % u denotes the solution in the own domain
338 % v denotes the solution in the neighbour domain
339 [e_u, d_u, gamm_u, H_b_u, I_u] = obj.get_boundary_ops(boundary);
340 [e_v, d_v, gamm_v, H_b_v, I_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary);
341
342 % Find the number of grid points along the interface
343 m_u = size(e_u, 2);
344 m_v = size(e_v, 2);
345
346 Hi = obj.Hi;
347 a = obj.a;
348
349 u = obj;
350 v = neighbour_scheme;
351
352 b1_u = gamm_u*u.lambda(I_u)./u.a11(I_u).^2;
353 b2_u = gamm_u*u.lambda(I_u)./u.a22(I_u).^2;
354 b1_v = gamm_v*v.lambda(I_v)./v.a11(I_v).^2;
355 b2_v = gamm_v*v.lambda(I_v)./v.a22(I_v).^2;
356
357 tau_u = -1./(4*b1_u) -1./(4*b2_u);
358 tau_v = -1./(4*b1_v) -1./(4*b2_v);
359
360 tau_u = tuning * spdiag(tau_u);
361 tau_v = tuning * spdiag(tau_v);
362 beta_u = tau_v;
363
364 % Build interpolation operators
365 intOps = interpOpSet(m_u, m_v, obj.order, neighbour_scheme.order);
366 Iu2v = intOps.Iu2v;
367 Iv2u = intOps.Iv2u;
368
369 closure = a*Hi*e_u*tau_u*H_b_u*e_u' + ...
370 a*Hi*e_u*H_b_u*Iv2u.bad*beta_u*Iu2v.good*e_u' + ...
371 a*1/2*Hi*d_u*H_b_u*e_u' + ...
372 -a*1/2*Hi*e_u*H_b_u*d_u';
373
374 penalty = -a*Hi*e_u*tau_u*H_b_u*Iv2u.good*e_v' + ...
375 -a*Hi*e_u*H_b_u*Iv2u.bad*beta_u*e_v' + ...
376 -a*1/2*Hi*d_u*H_b_u*Iv2u.good*e_v' + ...
377 -a*1/2*Hi*e_u*H_b_u*Iv2u.bad*d_v';
378
379 end
380
381 % Returns the boundary ops and sign for the boundary specified by the string boundary.
302 % The right boundary is considered the positive boundary 382 % The right boundary is considered the positive boundary
303 % 383 %
304 % I -- the indecies of the boundary points in the grid matrix 384 % I -- the indices of the boundary points in the grid matrix
305 function [e, d, gamm, H_b, I] = get_boundary_ops(obj, boundary) 385 function [e, d, gamm, H_b, I] = get_boundary_ops(obj, boundary)
306
307 % gridMatrix = zeros(obj.m(2),obj.m(1));
308 % gridMatrix(:) = 1:numel(gridMatrix);
309
310 ind = grid.funcToMatrix(obj.grid, 1:prod(obj.m)); 386 ind = grid.funcToMatrix(obj.grid, 1:prod(obj.m));
311 387
312 switch boundary 388 switch boundary
313 case 'w' 389 case 'w'
314 e = obj.e_w; 390 e = obj.e_w;