Mercurial > repos > public > sbplib
diff +sbp/D2VariableCompatible.m @ 1198:2924b3a9b921 feature/d2_compatible
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 |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/D2VariableCompatible.m Fri Aug 16 14:30:28 2019 -0700 @@ -0,0 +1,81 @@ +classdef D2VariableCompatible < sbp.OpSet + properties + D1 % SBP operator approximating first derivative + H % Norm matrix + HI % H^-1 + Q % Skew-symmetric matrix + e_l % Left boundary operator + e_r % Right boundary operator + D2 % SBP operator for second derivative + M % Norm matrix, second derivative + d1_l % Left boundary first derivative + d1_r % Right boundary first derivative + m % Number of grid points. + h % Step size + x % grid + borrowing % Struct with borrowing limits for different norm matrices + end + + methods + function obj = D2VariableCompatible(m,lim,order) + + x_l = lim{1}; + x_r = lim{2}; + L = x_r-x_l; + obj.h = L/(m-1); + obj.x = linspace(x_l,x_r,m)'; + + switch order + + case 6 + + [obj.H, obj.HI, obj.D1, D2, ... + ~, obj.e_l, obj.e_r, ~, ~, ~, ~, ~,... + d1_l, d1_r] = ... + sbp.implementations.d4_variable_6(m, obj.h); + + case 4 + [obj.H, obj.HI, obj.D1, D2, obj.e_l,... + obj.e_r, d1_l, d1_r] = ... + sbp.implementations.d2_variable_4(m,obj.h); + case 2 + [obj.H, obj.HI, obj.D1, D2, obj.e_l,... + obj.e_r, d1_l, d1_r] = ... + sbp.implementations.d2_variable_2(m,obj.h); + + otherwise + error('Invalid operator order %d.',order); + end + obj.borrowing.H11 = obj.H(1,1)/obj.h; % First element in H/h, + obj.borrowing.M.d1 = obj.H(1,1)/obj.h; % First element in H/h is borrowing also for M + obj.borrowing.R.delta_D = inf; % Because delta_D is zero, one can borrow infinitely much. + % This sets penalties of the form 1/borrowing to 0, which is + % the desired behaviour. + obj.m = m; + obj.M = []; + + D1 = obj.D1; + e_r = obj.e_r; + e_l = obj.e_l; + + % D2 = Hinv * (-M + br*er*d1r^T - bl*el*d1l^T); + % Replace d1' by e'*D1 in D2. + D2_compatible = @(b) D2(b) - obj.HI*(b(m)*e_r*d1_r' - b(m)*e_r*e_r'*D1) ... + + obj.HI*(b(1)*e_l*d1_l' - b(1)*e_l*e_l'*D1); + + obj.D2 = D2_compatible; + obj.d1_l = (e_l'*D1)'; + obj.d1_r = (e_r'*D1)'; + + end + function str = string(obj) + str = [class(obj) '_' num2str(obj.order)]; + end + end + +end + + + + +