Mercurial > repos > public > sbplib

comparison +scheme/LaplaceCurvilinear.m @ 720:07f8311374c6 feature/utux2D

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add interpolation to LaplaceCurveilinear.
author Martin Almquist <malmquist@stanford.edu>
date Wed, 14 Mar 2018 21:01:34 -0700
parents 33b962620e24
children f4595f14d696
comparison
equal deleted inserted replaced
719:b3f8fb9cefd2 720:07f8311374c6
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'))
57 60
58 m = g.size(); 61 m = g.size();
59 m_u = m(1); 62 m_u = m(1);
60 m_v = m(2); 63 m_v = m(2);
61 m_tot = g.N(); 64 m_tot = g.N();
207 % Misc. 210 % Misc.
208 obj.m = m; 211 obj.m = m;
209 obj.h = [h_u h_v]; 212 obj.h = [h_u h_v];
210 obj.order = order; 213 obj.order = order;
211 obj.grid = g; 214 obj.grid = g;
215 obj.interpolation_type = interpolation_type;
212 216
213 obj.a = a; 217 obj.a = a;
214 obj.b = b; 218 obj.b = b;
215 obj.a11 = a11; 219 obj.a11 = a11;
216 obj.a12 = a12; 220 obj.a12 = a12;
267 error('No such boundary condition: type = %s',type); 271 error('No such boundary condition: type = %s',type);
268 end 272 end
269 end 273 end
270 274
271 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) 275 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
276
277 [e_u, d_u, gamm_u, H_b_u, I_u] = obj.get_boundary_ops(boundary);
278 [e_v, d_v, gamm_v, H_b_v, I_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary);
279 Hi = obj.Hi;
280 a = obj.a;
281
282 m_u = length(H_b_u);
283 m_v = length(H_b_v);
284
285 grid_ratio = m_u/m_v;
286 if grid_ratio ~= 1
287
288 [ms, index] = sort([m_u, m_v]);
289 orders = [obj.order, neighbour_scheme.order];
290 orders = orders(index);
291
292 switch obj.interpolation_type
293 case 'MC'
294 interpOpSet = sbp.InterpMC(ms(1),ms(2),orders(1),orders(2));
295 if grid_ratio < 1
296 I_v2u_good = interpOpSet.IF2C;
297 I_v2u_bad = interpOpSet.IF2C;
298 I_u2v_good = interpOpSet.IC2F;
299 I_u2v_bad = interpOpSet.IC2F;
300 elseif grid_ratio > 1
301 I_v2u_good = interpOpSet.IC2F;
302 I_v2u_bad = interpOpSet.IC2F;
303 I_u2v_good = interpOpSet.IF2C;
304 I_u2v_bad = interpOpSet.IF2C;
305 end
306 case 'AWW'
307 %String 'C2F' indicates that ICF2 is more accurate.
308 interpOpSetF2C = sbp.InterpAWW(ms(1),ms(2),orders(1),orders(2),'F2C');
309 interpOpSetC2F = sbp.InterpAWW(ms(1),ms(2),orders(1),orders(2),'C2F');
310 if grid_ratio < 1
311 % Local is coarser than neighbour
312 I_v2u_good = interpOpSetF2C.IF2C;
313 I_v2u_bad = interpOpSetC2F.IF2C;
314 I_u2v_good = interpOpSetC2F.IC2F;
315 I_u2v_bad = interpOpSetF2C.IC2F;
316 elseif grid_ratio > 1
317 % Local is finer than neighbour
318 I_v2u_good = interpOpSetC2F.IC2F;
319 I_v2u_bad = interpOpSetF2C.IC2F;
320 I_u2v_good = interpOpSetF2C.IF2C;
321 I_u2v_bad = interpOpSetC2F.IF2C;
322 end
323 otherwise
324 error(['Interpolation type ' obj.interpolation_type ...
325 ' is not available.' ]);
326 end
327
328 else
329 % No interpolation required
330 I_v2u_good = speye(m_u,m_u);
331 I_v2u_bad = speye(m_u,m_u);
332 I_u2v_good = speye(m_u,m_u);
333 I_u2v_bad = speye(m_u,m_u);
334 end
335
272 % u denotes the solution in the own domain 336 % u denotes the solution in the own domain
273 % v denotes the solution in the neighbour domain 337 % v denotes the solution in the neighbour domain
274 tuning = 1.2; 338 tuning = 1.2;
275 % tuning = 20.2; 339 % 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; 340 u = obj;
280 v = neighbour_scheme; 341 v = neighbour_scheme;
281 342
282 b1_u = gamm_u*u.lambda(I_u)./u.a11(I_u).^2; 343 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; 344 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; 345 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; 346 b2_v = gamm_v*v.lambda(I_v)./v.a22(I_v).^2;
286 347
287 tau1 = -1./(4*b1_u) -1./(4*b1_v) -1./(4*b2_u) -1./(4*b2_v); 348 tau_u = -1./(4*b1_u) -1./(4*b2_u);
288 tau1 = tuning * spdiag(tau1); 349 tau_v = -1./(4*b1_v) -1./(4*b2_v);
289 tau2 = 1/2; 350
290 351 tau_u = tuning * spdiag(tau_u);
291 sig1 = -1/2; 352 tau_v = tuning * spdiag(tau_v);
292 sig2 = 0; 353 beta_u = tau_v;
293 354
294 tau = (e_u*tau1 + tau2*d_u)*H_b_u; 355 closure = a*Hi*e_u*tau_u*H_b_u*e_u' + ...
295 sig = (sig1*e_u + sig2*d_u)*H_b_u; 356 a*Hi*e_u*H_b_u*I_v2u_bad*beta_u*I_u2v_good*e_u' + ...
296 357 a*1/2*Hi*d_u*H_b_u*e_u' + ...
297 closure = obj.a*obj.Hi*( tau*e_u' + sig*d_u'); 358 -a*1/2*Hi*e_u*H_b_u*d_u';
298 penalty = obj.a*obj.Hi*(-tau*e_v' + sig*d_v'); 359
360 penalty = -a*Hi*e_u*tau_u*H_b_u*I_v2u_good*e_v' + ...
361 -a*Hi*e_u*H_b_u*I_v2u_bad*beta_u*e_v' + ...
362 -a*1/2*Hi*d_u*H_b_u*I_v2u_good*e_v' + ...
363 -a*1/2*Hi*e_u*H_b_u*I_v2u_bad*d_v';
364
365
299 end 366 end
300 367
301 % Ruturns the boundary ops and sign for the boundary specified by the string boundary. 368 % Ruturns the boundary ops and sign for the boundary specified by the string boundary.
302 % The right boundary is considered the positive boundary 369 % The right boundary is considered the positive boundary
303 % 370 %