Mercurial > repos > public > sbplib

comparison +scheme/LaplaceCurvilinear.m @ 756:f891758ad7a4 feature/d1_staggered

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with feature/utux2d.
author Martin Almquist <malmquist@stanford.edu>
date Sat, 16 Jun 2018 14:30:45 -0700
parents f4595f14d696
children 14b093a344eb
comparison
equal deleted inserted replaced
755:14f0058356f2 756:f891758ad7a4
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
42 interpolation_type
41 end 43 end
42 44
43 methods 45 methods
44 % Implements a*div(b*grad(u)) as a SBP scheme 46 % 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?) 47 % TODO: Implement proper H, it should be the real physical quadrature, the logic quadrature may be but in a separate variable (H_logic?)
46 48
47 function obj = LaplaceCurvilinear(g ,order, a, b, opSet) 49 function obj = LaplaceCurvilinear(g ,order, a, b, opSet, interpolation_type)
48 default_arg('opSet',@sbp.D2Variable); 50 default_arg('opSet',@sbp.D2Variable);
49 default_arg('a', 1); 51 default_arg('a', 1);
50 default_arg('b', 1); 52 default_arg('b', 1);
53 default_arg('interpolation_type','AWW');
51 54
52 if b ~=1 55 if b ~=1
53 error('Not implemented yet') 56 error('Not implemented yet')
54 end 57 end
55 58
56 assert(isa(g, 'grid.Curvilinear')) 59 % assert(isa(g, 'grid.Curvilinear'))
60 if isa(a, 'function_handle')
61 a = grid.evalOn(g, a);
62 a = spdiag(a);
63 end
57 64
58 m = g.size(); 65 m = g.size();
59 m_u = m(1); 66 m_u = m(1);
60 m_v = m(2); 67 m_v = m(2);
61 m_tot = g.N(); 68 m_tot = g.N();
207 % Misc. 214 % Misc.
208 obj.m = m; 215 obj.m = m;
209 obj.h = [h_u h_v]; 216 obj.h = [h_u h_v];
210 obj.order = order; 217 obj.order = order;
211 obj.grid = g; 218 obj.grid = g;
219 obj.interpolation_type = interpolation_type;
212 220
213 obj.a = a; 221 obj.a = a;
214 obj.b = b; 222 obj.b = b;
215 obj.a11 = a11; 223 obj.a11 = a11;
216 obj.a12 = a12; 224 obj.a12 = a12;
267 error('No such boundary condition: type = %s',type); 275 error('No such boundary condition: type = %s',type);
268 end 276 end
269 end 277 end
270 278
271 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) 279 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
280
281 [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);
283 Hi = obj.Hi;
284 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
272 % u denotes the solution in the own domain 340 % u denotes the solution in the own domain
273 % v denotes the solution in the neighbour domain 341 % v denotes the solution in the neighbour domain
274 tuning = 1.2; 342 tuning = 1.2;
275 % tuning = 20.2; 343 % tuning = 20.2;
276 [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);
278
279 u = obj; 344 u = obj;
280 v = neighbour_scheme; 345 v = neighbour_scheme;
281 346
282 b1_u = gamm_u*u.lambda(I_u)./u.a11(I_u).^2; 347 b1_u = gamm_u*u.lambda(I_u)./u.a11(I_u).^2;
283 b2_u = gamm_u*u.lambda(I_u)./u.a22(I_u).^2; 348 b2_u = gamm_u*u.lambda(I_u)./u.a22(I_u).^2;
284 b1_v = gamm_v*v.lambda(I_v)./v.a11(I_v).^2; 349 b1_v = gamm_v*v.lambda(I_v)./v.a11(I_v).^2;
285 b2_v = gamm_v*v.lambda(I_v)./v.a22(I_v).^2; 350 b2_v = gamm_v*v.lambda(I_v)./v.a22(I_v).^2;
286 351
287 tau1 = -1./(4*b1_u) -1./(4*b1_v) -1./(4*b2_u) -1./(4*b2_v); 352 tau_u = -1./(4*b1_u) -1./(4*b2_u);
288 tau1 = tuning * spdiag(tau1); 353 tau_v = -1./(4*b1_v) -1./(4*b2_v);
289 tau2 = 1/2; 354
290 355 tau_u = tuning * spdiag(tau_u);
291 sig1 = -1/2; 356 tau_v = tuning * spdiag(tau_v);
292 sig2 = 0; 357 beta_u = tau_v;
293 358
294 tau = (e_u*tau1 + tau2*d_u)*H_b_u; 359 closure = a*Hi*e_u*tau_u*H_b_u*e_u' + ...
295 sig = (sig1*e_u + sig2*d_u)*H_b_u; 360 a*Hi*e_u*H_b_u*I_v2u_bad*beta_u*I_u2v_good*e_u' + ...
296 361 a*1/2*Hi*d_u*H_b_u*e_u' + ...
297 closure = obj.a*obj.Hi*( tau*e_u' + sig*d_u'); 362 -a*1/2*Hi*e_u*H_b_u*d_u';
298 penalty = obj.a*obj.Hi*(-tau*e_v' + sig*d_v'); 363
364 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' + ...
366 -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';
368
369
299 end 370 end
300 371
301 % Ruturns the boundary ops and sign for the boundary specified by the string boundary. 372 % Ruturns the boundary ops and sign for the boundary specified by the string boundary.
302 % The right boundary is considered the positive boundary 373 % The right boundary is considered the positive boundary
303 % 374 %