Mercurial > repos > public > sbplib
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 |