Mercurial > repos > public > sbplib

comparison +scheme/Heat2dCurvilinear.m @ 1100:27aaf8646a80 feature/timesteppers

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with default
author Jonatan Werpers <jonatan@werpers.com>
date Tue, 09 Apr 2019 21:48:30 +0200
parents 8d73fcdb07a5
children
comparison
equal deleted inserted replaced
1099:d4fe089b2c4a 1100:27aaf8646a80
1 classdef Heat2dCurvilinear < scheme.Scheme 1 classdef Heat2dCurvilinear < scheme.Scheme
2 2
3 % Discretizes the Laplacian with variable coefficent, curvilinear, 3 % Discretizes the Laplacian with variable coefficent, curvilinear,
4 % in the Heat equation way (i.e., the discretization matrix is not necessarily 4 % in the Heat equation way (i.e., the discretization matrix is not necessarily
5 % symmetric) 5 % symmetric)
6 % u_t = div * (kappa * grad u ) 6 % u_t = div * (kappa * grad u )
7 % opSet should be cell array of opSets, one per dimension. This 7 % opSet should be cell array of opSets, one per dimension. This
8 % is useful if we have periodic BC in one direction. 8 % is useful if we have periodic BC in one direction.
9 9
10 properties 10 properties
11 m % Number of points in each direction, possibly a vector 11 m % Number of points in each direction, possibly a vector
27 27
28 H, Hi % Inner products 28 H, Hi % Inner products
29 e_l, e_r 29 e_l, e_r
30 d1_l, d1_r % Normal derivatives at the boundary 30 d1_l, d1_r % Normal derivatives at the boundary
31 alpha % Vector of borrowing constants 31 alpha % Vector of borrowing constants
32 32
33 % Boundary inner products 33 % Boundary inner products
34 H_boundary_l, H_boundary_r 34 H_boundary_l, H_boundary_r
35 35
36 % Metric coefficients 36 % Metric coefficients
37 b % Cell matrix of size dim x dim 37 b % Cell matrix of size dim x dim
38 J, Ji 38 J, Ji
39 beta % Cell array of scale factors 39 beta % Cell array of scale factors
107 xmax = max(ops{1}.x); 107 xmax = max(ops{1}.x);
108 ymax = max(ops{2}.x); 108 ymax = max(ops{2}.x);
109 opSetMetric{1} = sbp.D2Variable(m(1), {0, xmax}, order); 109 opSetMetric{1} = sbp.D2Variable(m(1), {0, xmax}, order);
110 opSetMetric{2} = sbp.D2Variable(m(2), {0, ymax}, order); 110 opSetMetric{2} = sbp.D2Variable(m(2), {0, ymax}, order);
111 D1Metric{1} = kron(opSetMetric{1}.D1, I{2}); 111 D1Metric{1} = kron(opSetMetric{1}.D1, I{2});
112 D1Metric{2} = kron(I{1}, opSetMetric{2}.D1); 112 D1Metric{2} = kron(I{1}, opSetMetric{2}.D1);
113 113
114 x_xi = D1Metric{1}*x; 114 x_xi = D1Metric{1}*x;
115 x_eta = D1Metric{2}*x; 115 x_eta = D1Metric{2}*x;
116 y_xi = D1Metric{1}*y; 116 y_xi = D1Metric{1}*y;
117 y_eta = D1Metric{2}*y; 117 y_eta = D1Metric{2}*y;
155 obj.d1_r{2} = kron(I{1},d1_r{2}); 155 obj.d1_r{2} = kron(I{1},d1_r{2});
156 156
157 % D2 coefficients 157 % D2 coefficients
158 kappa_coeff = cell(dim,dim); 158 kappa_coeff = cell(dim,dim);
159 for j = 1:dim 159 for j = 1:dim
160 obj.D2_kappa{j} = sparse(m_tot,m_tot); 160 obj.D2_kappa{j} = sparse(m_tot,m_tot);
161 kappa_coeff{j} = sparse(m_tot,1); 161 kappa_coeff{j} = sparse(m_tot,1);
162 for i = 1:dim 162 for i = 1:dim
163 kappa_coeff{j} = kappa_coeff{j} + b{i,j}*J*b{i,j}*kappa; 163 kappa_coeff{j} = kappa_coeff{j} + b{i,j}*J*b{i,j}*kappa;
164 end 164 end
165 end 165 end
268 function [closure, penalty] = boundary_condition(obj, boundary, type, symmetric, tuning) 268 function [closure, penalty] = boundary_condition(obj, boundary, type, symmetric, tuning)
269 default_arg('type','Neumann'); 269 default_arg('type','Neumann');
270 default_arg('symmetric', false); 270 default_arg('symmetric', false);
271 default_arg('tuning',1.2); 271 default_arg('tuning',1.2);
272 272
273 % j is the coordinate direction of the boundary 273 % nj: outward unit normal component.
274 % nj: outward unit normal component.
275 % nj = -1 for west, south, bottom boundaries 274 % nj = -1 for west, south, bottom boundaries
276 % nj = 1 for east, north, top boundaries 275 % nj = 1 for east, north, top boundaries
277 [j, nj] = obj.get_boundary_number(boundary); 276 nj = obj.getBoundarySign(boundary);
278 switch nj 277
279 case 1 278 Hi = obj.Hi;
280 e = obj.e_r{j}; 279 [e, flux] = obj.getBoundaryOperator({'e', 'flux'}, boundary);
281 flux = obj.flux_r{j}; 280 H_gamma = obj.getBoundaryQuadrature(boundary);
282 H_gamma = obj.H_boundary_r{j}; 281 alpha = obj.getBoundaryBorrowing(boundary);
283 case -1
284 e = obj.e_l{j};
285 flux = obj.flux_l{j};
286 H_gamma = obj.H_boundary_l{j};
287 end
288 282
289 Hi = obj.Hi; 283 Hi = obj.Hi;
290 Ji = obj.Ji; 284 Ji = obj.Ji;
291 KAPPA = obj.KAPPA; 285 KAPPA = obj.KAPPA;
292 kappa_gamma = e'*KAPPA*e; 286 kappa_gamma = e'*KAPPA*e;
293 h = obj.h(j);
294 alpha = h*obj.alpha(j);
295 287
296 switch type 288 switch type
297 289
298 % Dirichlet boundary condition 290 % Dirichlet boundary condition
299 case {'D','d','dirichlet','Dirichlet'} 291 case {'D','d','dirichlet','Dirichlet'}
300 292
301 if ~symmetric 293 if ~symmetric
302 closure = -Ji*Hi*flux'*e*H_gamma*(e' ); 294 closure = -Ji*Hi*flux'*e*H_gamma*(e' );
303 penalty = Ji*Hi*flux'*e*H_gamma; 295 penalty = Ji*Hi*flux'*e*H_gamma;
304 else 296 else
305 closure = Ji*Hi*flux'*e*H_gamma*(e' )... 297 closure = Ji*Hi*flux'*e*H_gamma*(e' )...
306 -tuning*2/alpha*Ji*Hi*e*kappa_gamma*H_gamma*(e' ) ; 298 -tuning*2/alpha*Ji*Hi*e*kappa_gamma*H_gamma*(e' ) ;
307 penalty = -Ji*Hi*flux'*e*H_gamma ... 299 penalty = -Ji*Hi*flux'*e*H_gamma ...
308 +tuning*2/alpha*Ji*Hi*e*kappa_gamma*H_gamma; 300 +tuning*2/alpha*Ji*Hi*e*kappa_gamma*H_gamma;
309 end 301 end
310 302
311 % Normal flux boundary condition 303 % Normal flux boundary condition
312 case {'N','n','neumann','Neumann'} 304 case {'N','n','neumann','Neumann'}
313 closure = -Ji*Hi*e*H_gamma*(e'*flux ); 305 closure = -Ji*Hi*e*H_gamma*(e'*flux );
314 penalty = Ji*Hi*e*H_gamma; 306 penalty = Ji*Hi*e*H_gamma;
315 307
316 % Unknown boundary condition 308 % Unknown boundary condition
317 otherwise 309 otherwise
318 error('No such boundary condition: type = %s',type); 310 error('No such boundary condition: type = %s',type);
319 end 311 end
323 % u denotes the solution in the own domain 315 % u denotes the solution in the own domain
324 % v denotes the solution in the neighbour domain 316 % v denotes the solution in the neighbour domain
325 error('Interface not implemented'); 317 error('Interface not implemented');
326 end 318 end
327 319
328 % Returns the coordinate number and outward normal component for the boundary specified by the string boundary. 320 % Returns the boundary operator op for the boundary specified by the string boundary.
329 function [j, nj] = get_boundary_number(obj, boundary) 321 % op -- string or a cell array of strings
330 322 % boundary -- string
331 switch boundary 323 function varargout = getBoundaryOperator(obj, op, boundary)
332 case {'w','W','west','West', 'e', 'E', 'east', 'East'} 324 assertIsMember(boundary, {'w', 'e', 's', 'n'})
333 j = 1; 325
334 case {'s','S','south','South', 'n', 'N', 'north', 'North'} 326 if ~iscell(op)
335 j = 2; 327 op = {op};
336 otherwise 328 end
337 error('No such boundary: boundary = %s',boundary); 329
338 end 330 for i = 1:numel(op)
339 331 switch op{i}
340 switch boundary
341 case {'w','W','west','West','s','S','south','South'}
342 nj = -1;
343 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
344 nj = 1;
345 end
346 end
347
348 % Returns the coordinate number and outward normal component for the boundary specified by the string boundary.
349 function [return_op] = get_boundary_operator(obj, op, boundary)
350
351 switch boundary
352 case {'w','W','west','West', 'e', 'E', 'east', 'East'}
353 j = 1;
354 case {'s','S','south','South', 'n', 'N', 'north', 'North'}
355 j = 2;
356 otherwise
357 error('No such boundary: boundary = %s',boundary);
358 end
359
360 switch op
361 case 'e' 332 case 'e'
362 switch boundary 333 switch boundary
363 case {'w','W','west','West','s','S','south','South'} 334 case 'w'
364 return_op = obj.e_l{j}; 335 e = obj.e_l{1};
365 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} 336 case 'e'
366 return_op = obj.e_r{j}; 337 e = obj.e_r{1};
338 case 's'
339 e = obj.e_l{2};
340 case 'n'
341 e = obj.e_r{2};
367 end 342 end
343 varargout{i} = e;
344
368 case 'd' 345 case 'd'
369 switch boundary 346 switch boundary
370 case {'w','W','west','West','s','S','south','South'} 347 case 'w'
371 return_op = obj.d1_l{j}; 348 d = obj.d1_l{1};
372 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} 349 case 'e'
373 return_op = obj.d1_r{j}; 350 d = obj.d1_r{1};
351 case 's'
352 d = obj.d1_l{2};
353 case 'n'
354 d = obj.d1_r{2};
374 end 355 end
375 otherwise 356 varargout{i} = d;
376 error(['No such operator: operatr = ' op]); 357
377 end 358 case 'flux'
378 359 switch boundary
360 case 'w'
361 flux = obj.flux_l{1};
362 case 'e'
363 flux = obj.flux_r{1};
364 case 's'
365 flux = obj.flux_l{2};
366 case 'n'
367 flux = obj.flux_r{2};
368 end
369 varargout{i} = flux;
370 end
371 end
372 end
373
374 % Returns square boundary quadrature matrix, of dimension
375 % corresponding to the number of boundary points
376 %
377 % boundary -- string
378 function H_b = getBoundaryQuadrature(obj, boundary)
379 assertIsMember(boundary, {'w', 'e', 's', 'n'})
380
381 switch boundary
382 case 'w'
383 H_b = obj.H_boundary_l{1};
384 case 'e'
385 H_b = obj.H_boundary_r{1};
386 case 's'
387 H_b = obj.H_boundary_l{2};
388 case 'n'
389 H_b = obj.H_boundary_r{2};
390 end
391 end
392
393 % Returns the boundary sign. The right boundary is considered the positive boundary
394 % boundary -- string
395 function s = getBoundarySign(obj, boundary)
396 assertIsMember(boundary, {'w', 'e', 's', 'n'})
397
398 switch boundary
399 case {'e','n'}
400 s = 1;
401 case {'w','s'}
402 s = -1;
403 end
404 end
405
406 % Returns borrowing constant gamma*h
407 % boundary -- string
408 function gamm = getBoundaryBorrowing(obj, boundary)
409 assertIsMember(boundary, {'w', 'e', 's', 'n'})
410
411 switch boundary
412 case {'w','e'}
413 gamm = obj.h(1)*obj.alpha(1);
414 case {'s','n'}
415 gamm = obj.h(2)*obj.alpha(2);
416 end
379 end 417 end
380 418
381 function N = size(obj) 419 function N = size(obj)
382 N = prod(obj.m); 420 N = prod(obj.m);
383 end 421 end