Mercurial > repos > public > sbplib

comparison +sbp/D2Nonequidistant.m @ 1325:1b0f2415237f feature/D2_boundary_opt

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add variable coefficient boundary-optimized second derivatives.
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Sun, 13 Feb 2022 19:32:34 +0100
parents
children 7ab7d42a5b24
comparison
equal deleted inserted replaced
1301:8978521b0f06 1325:1b0f2415237f
1 classdef D2Nonequidistant < sbp.OpSet
2 % Implements the boundary optimized variable coefficient
3 % second derivative.
4 %
5 % The boundary closure uses the first and last rows of the
6 % boundary-optimized D1 operator, i.e. the operators are
7 % fully compatible.
8 properties
9 H % Norm matrix
10 HI % H^-1
11 D1 % SBP operator approximating first derivative
12 D2 % SBP operator approximating first derivative
13 DI % Dissipation operator
14 e_l % Left boundary operator
15 e_r % Right boundary operator
16 d1_l % Left boundary first derivative
17 d1_r % Right boundary first derivative
18 m % Number of grid points.
19 h % Step size
20 x % grid
21 borrowing % Struct with borrowing limits for different norm matrices
22 options % Struct holding options used to create the operator
23 end
24
25 methods
26
27 function obj = D2Nonequidistant(m, lim, order, options)
28 % m - number of gridpoints
29 % lim - cell array holding the limits of the domain
30 % order - order of the operator
31 % options - struct holding options used to construct the operator
32 % struct.stencil_width: {'minimal', 'nonminimal', 'wide'}
33 % minimal: minimal compatible stencil width (default)
34 % nonminimal: a few additional stencil points compared to minimal
35 % wide: wide stencil obtained by applying D1 twice
36 % struct.AD: {'op', 'upwind'}
37 % 'op': order-preserving AD (preserving interior stencil order) (default)
38 % 'upwind': upwind AD (order-1 upwind interior stencil)
39 % struct.variable_coeffs: {true, false}
40 % true: obj.D2 is a function handle D2(c) returning a matrix
41 % for coefficient vector c (default)
42 % false: obj.D2 is a matrix.
43 default_arg('options', struct);
44 default_field(options,'stencil_width','minimal');
45 default_field(options,'AD','op');
46 default_field(options,'variable_coeffs',true);
47 [x, h] = sbp.grid.accurateBoundaryOptimizedGrid(lim, m, order);
48 switch order
49 case 4
50 [obj.H, obj.HI, obj.D1, obj.D2, obj.DI] = sbp.implementations.d2_noneq_variable_4(m, h, options);
51 case 6
52 % [obj.H, obj.HI, obj.D1, D2, M, Q, obj.e_l, obj.e_r, obj.d1_l, obj.d1_r] = sbp.implementations.d2_noneq_6(m,h,2);
53 % obj.D2 = @(c)D2;
54 [obj.H, obj.HI, obj.D1, obj.D2, obj.DI] = sbp.implementations.d2_noneq_variable_6(m, h, options);
55 case 8
56 % [obj.H, obj.HI, obj.D1, D2, M, Q, obj.e_l, obj.e_r, obj.d1_l, obj.d1_r] = sbp.implementations.d2_noneq_8(m,h,1);
57 % obj.D2 = @(c)D2;
58 [obj.H, obj.HI, obj.D1, obj.D2, obj.DI] = sbp.implementations.d2_noneq_variable_8(m, h, options);
59 case 10
60 [obj.H, obj.HI, obj.D1, obj.D2, obj.DI] = sbp.implementations.d2_noneq_variable_10(m, h, options);
61 case 12
62 [obj.H, obj.HI, obj.D1, obj.D2, obj.DI] = sbp.implementations.d2_noneq_variable_12(m, h, options);
63 otherwise
64 error('Invalid operator order %d.', order);
65 end
66
67 if ~options.variable_coeffs
68 obj.D2 = obj.D2(ones(m,1));
69 end
70
71 % Boundary operators
72 obj.e_l = sparse(m, 1); obj.e_l(1) = 1;
73 obj.e_r = sparse(m, 1); obj.e_r(m) = 1;
74 obj.d1_l = (obj.e_l' * obj.D1)';
75 obj.d1_r = (obj.e_r' * obj.D1)';
76
77 % Borrowing coefficients
78 obj.borrowing.H11 = obj.H(1, 1) / h; % First element in H/h,
79 obj.borrowing.M.d1 = obj.H(1, 1) / h; % First element in H/h is borrowing also for M
80 obj.borrowing.R.delta_D = inf;
81
82 % grid data
83 obj.x = x;
84 obj.h = h;
85 obj.m = m;
86
87 % misc
88 obj.options = options;
89 end
90
91 function str = string(obj)
92 str = [class(obj) '_' num2str(obj.order)];
93 end
94
95 end
96
97 end