Mercurial > repos > public > sbplib

comparison +sbp/D2VariableCompatible.m @ 1198:2924b3a9b921 feature/d2_compatible

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add OpSet for fully compatible D2Variable, created from regular D2Variable by replacing d1 by first row of D1. Formal reduction by one order of accuracy at the boundary point.
author Martin Almquist <malmquist@stanford.edu>
date Fri, 16 Aug 2019 14:30:28 -0700
parents
children a3d9567d9004
comparison
equal deleted inserted replaced
1116:33c378e508d2 1198:2924b3a9b921
1 classdef D2VariableCompatible < sbp.OpSet
2 properties
3 D1 % SBP operator approximating first derivative
4 H % Norm matrix
5 HI % H^-1
6 Q % Skew-symmetric matrix
7 e_l % Left boundary operator
8 e_r % Right boundary operator
9 D2 % SBP operator for second derivative
10 M % Norm matrix, second derivative
11 d1_l % Left boundary first derivative
12 d1_r % Right boundary first derivative
13 m % Number of grid points.
14 h % Step size
15 x % grid
16 borrowing % Struct with borrowing limits for different norm matrices
17 end
18
19 methods
20 function obj = D2VariableCompatible(m,lim,order)
21
22 x_l = lim{1};
23 x_r = lim{2};
24 L = x_r-x_l;
25 obj.h = L/(m-1);
26 obj.x = linspace(x_l,x_r,m)';
27
28 switch order
29
30 case 6
31
32 [obj.H, obj.HI, obj.D1, D2, ...
33 ~, obj.e_l, obj.e_r, ~, ~, ~, ~, ~,...
34 d1_l, d1_r] = ...
35 sbp.implementations.d4_variable_6(m, obj.h);
36
37 case 4
38 [obj.H, obj.HI, obj.D1, D2, obj.e_l,...
39 obj.e_r, d1_l, d1_r] = ...
40 sbp.implementations.d2_variable_4(m,obj.h);
41 case 2
42 [obj.H, obj.HI, obj.D1, D2, obj.e_l,...
43 obj.e_r, d1_l, d1_r] = ...
44 sbp.implementations.d2_variable_2(m,obj.h);
45
46 otherwise
47 error('Invalid operator order %d.',order);
48 end
49 obj.borrowing.H11 = obj.H(1,1)/obj.h; % First element in H/h,
50 obj.borrowing.M.d1 = obj.H(1,1)/obj.h; % First element in H/h is borrowing also for M
51 obj.borrowing.R.delta_D = inf; % Because delta_D is zero, one can borrow infinitely much.
52 % This sets penalties of the form 1/borrowing to 0, which is
53 % the desired behaviour.
54 obj.m = m;
55 obj.M = [];
56
57 D1 = obj.D1;
58 e_r = obj.e_r;
59 e_l = obj.e_l;
60
61 % D2 = Hinv * (-M + br*er*d1r^T - bl*el*d1l^T);
62 % Replace d1' by e'*D1 in D2.
63 D2_compatible = @(b) D2(b) - obj.HI*(b(m)*e_r*d1_r' - b(m)*e_r*e_r'*D1) ...
64 + obj.HI*(b(1)*e_l*d1_l' - b(1)*e_l*e_l'*D1);
65
66 obj.D2 = D2_compatible;
67 obj.d1_l = (e_l'*D1)';
68 obj.d1_r = (e_r'*D1)';
69
70 end
71 function str = string(obj)
72 str = [class(obj) '_' num2str(obj.order)];
73 end
74 end
75
76 end
77
78
79
80
81