Mercurial > repos > public > sbplib
view +util/calc_borrowing.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 | d24869abc7cd |
children |
line wrap: on
line source
function calc_borrowing(m, h) default_arg('m',100); default_arg('h',1); operators = { { 'd4_lonely', getM4_lonely, { {4, 'min_boundary_points'}, {6, 'min_boundary_points'}, {6, '2'}, {6, '3'}, {8, 'min_boundary_points'}, {8, 'higher_boundary_order'}, } }, { 'd4_variable', { {2}, {4}, {6}, } } % BORKEN BAD IDEA } for i = 1:operators baseName = operators{i}{1}; postFixes = operators{i}{2}; for pf = postFixes [a2, a3] = borrowFromD4(m, h, l{:}); end end lonely = { {4, 'min_boundary_points'}, {6, 'min_boundary_points'}, {6, '2'}, {6, '3'}, {8, 'min_boundary_points'}, {8, 'higher_boundary_order'}, }; for i = 1:length(lonely) l = lonely{i}; [a2, a3] = d4_lonely(m, h, l{:}); fprintf('d4_lonely %d %s\n', l{:}) fprintf('\t alpha_II = %f\n', a2) fprintf('\t alpha_III = %f\n', a3) fprintf('\n') end variable = { {2}, {4}, {6}, }; for i = 1:length(variable) l = variable{i}; [a2, a3] = d4_variable(m, h, l{:}); fprintf('d4_variable %d\n', l{:}) fprintf('\t alpha_II = %f\n', a2) fprintf('\t alpha_III = %f\n', a3) fprintf('\n') end %% 4th order non-compatible [H, HI, D1, D2, D3, D4, e_1, e_m, M, M4,Q, Q3, S2_1, S2_m, S3_1, S3_m, S_1, S_m] = sbp.higher4(m,h); S1 = S_1*S_1' + S_m*S_m'; S2 = S2_1*S2_1' + S2_m*S2_m'; S3 = S3_1*S3_1' + S3_m*S3_m'; alpha_I = util.matrixborrow(M4, h^-1*S1 ); alpha_II = util.matrixborrow(M4, h*S2 ); alpha_III = util.matrixborrow(M4, h^3*S3); fprintf('4th order non-compatible\n') fprintf('alpha_I1: %.10f\n',alpha_I) fprintf('alpha_II: %.10f\n',alpha_II) fprintf('alpha_III: %.10f\n',alpha_III) fprintf('\n') %% 6th order non-compatible [H, HI, D1, D2, D3, D4, e_1, e_m, M, M4,Q, Q3, S2_1, S2_m, S3_1, S3_m, S_1, S_m] = sbp.higher6(m,h); S1 = S_1*S_1' + S_m*S_m'; S2 = S2_1*S2_1' + S2_m*S2_m'; S3 = S3_1*S3_1' + S3_m*S3_m'; alpha_II = util.matrixborrow(M4, h*S2 ); alpha_III = util.matrixborrow(M4, h^3*S3); fprintf('6th order non-compatible\n') fprintf('alpha_II: %.10f\n',alpha_II) fprintf('alpha_III: %.10f\n',alpha_III) fprintf('\n') %% 2nd order compatible [H, HI, D1, D4, e_1, e_m, M4, Q, S2_1, S2_m, S3_1, S3_m, S_1, S_m] = sbp.higher_compatible2(m,h); S1 = S_1*S_1' + S_m*S_m'; S2 = S2_1*S2_1' + S2_m*S2_m'; S3 = S3_1*S3_1' + S3_m*S3_m'; alpha_II = util.matrixborrow(M4, h*S2 ); alpha_III = util.matrixborrow(M4, h^3*S3); fprintf('2nd order compatible\n') fprintf('alpha_II: %.10f\n',alpha_II) fprintf('alpha_III: %.10f\n',alpha_III) fprintf('\n') %% 4th order compatible [H, HI, D1, D4, e_1, e_m, M4, Q, S2_1, S2_m, S3_1, S3_m, S_1, S_m] = sbp.higher_compatible4(m,h); S1 = S_1*S_1' + S_m*S_m'; S2 = S2_1*S2_1' + S2_m*S2_m'; S3 = S3_1*S3_1' + S3_m*S3_m'; alpha_II = util.matrixborrow(M4, h*S2 ); alpha_III = util.matrixborrow(M4, h^3*S3); fprintf('4th order compatible\n') fprintf('alpha_II: %.10f\n',alpha_II) fprintf('alpha_III: %.10f\n',alpha_III) fprintf('\n') %% 6th order compatible [H, HI, D1, D4, e_1, e_m, M4, Q, S2_1, S2_m, S3_1, S3_m, S_1, S_m] = sbp.higher_compatible6(m,h); S1 = S_1*S_1' + S_m*S_m'; S2 = S2_1*S2_1' + S2_m*S2_m'; S3 = S3_1*S3_1' + S3_m*S3_m'; alpha_II = util.matrixborrow(M4, h*S2 ); alpha_III = util.matrixborrow(M4, h^3*S3); fprintf('6th order compatible\n') fprintf('alpha_II: %.10f\n',alpha_II) fprintf('alpha_III: %.10f\n',alpha_III) fprintf('\n') % Ordinary for order = [2 4 6 8 10] op = sbp.Ordinary(m,h, order); S_1 = op.boundary.S_1; S_m = op.boundary.S_m; M = op.norms.M; S1 = S_1*S_1' + S_m*S_m'; alpha = util.matrixborrow(M, h*S1); fprintf('%dth order Ordinary\n', order) fprintf('alpha: %.10f\n', alpha) fprintf('\n') end end function [alpha_II, alpha_III] = d4_lonely(m, h, order, modifier) default_arg('modifier', []) func = sprintf('sbp.implementations.d4_lonely_%d', order); if ~isempty(modifier) func = sprintf('%s_%s', func, modifier); end funcCall = sprintf('%s(%s,%s)', func, toString(m), toString(h)); [H, HI, D4, e_l, e_r, M4, d2_l, d2_r, d3_l, d3_r, d1_l, d1_r] = eval(funcCall); d2d2 = d2_l*d2_l' + d2_r*d2_r'; alpha_II = util.matrixborrow(M4, h*d2d2); d3d3 = d3_l*d3_l' + d3_r*d3_r'; alpha_III = util.matrixborrow(M4, h^3*d3d3); end function [alpha_II, alpha_III] = d4_variable(m, h, order) default_arg('modifier', []) func = sprintf('sbp.implementations.d4_variable_%d', order); funcCall = sprintf('%s(%s,%s)', func, toString(m), toString(h)); [H, HI, D1, D2, D4, e_l, e_r, M4, d2_l, d2_r, d3_l, d3_r, d1_l, d1_r] = eval(funcCall); d2d2 = d2_l*d2_l' + d2_r*d2_r'; alpha_II = util.matrixborrow(M4, h*d2d2); d3d3 = d3_l*d3_l' + d3_r*d3_r'; alpha_III = util.matrixborrow(M4, h^3*d3d3); end function [d2_l, d2_r, d3_l, d3_r, M4] = getM4_lonely(m, h, order, modifier) fStr = getFunctionCallStr('d4_lonely', {order, modifier}, {m ,h}); [H, HI, D4, e_l, e_r, M4, d2_l, d2_r, d3_l, d3_r, d1_l, d1_r] = eval(funcCall); end % Calculates the borrowing constants for a D4 operator. % getM4 is a function handle on the form % [d2_l, d2_r, d3_l, d3_r, M4] = getM4(m,h) function [a2, a3] = borrowFromD4(m, h, getM4) [d2_l, d2_r, d3_l, d3_r, M4] = getM4(m, h); d2d2 = d2_l*d2_l' + d2_r*d2_r'; a2 = util.matrixborrow(M4, h*d2d2); d3d3 = d3_l*d3_l' + d3_r*d3_r'; a3 = util.matrixborrow(M4, h^3*d3d3); end function funcCallStr = getFunctionCallStr(baseName, postFix, parameters) default_arg('postFix', []) default_arg('parameters', []) funcCallStr = sprintf('sbp.implementations.%s', baseName); for i = 1:length(postFix) if ischar(postFix{i}) funcCallStr = [funcCallStr '_' postFix{i}]; else funcCallStr = [funcCallStr '_' toString(postFix{i})]; end end if isempty(parameters) return end funcCallStr = [funcCallStr '(' toString(parameters{1})]; for i = 2:length(parameters) funcCallStr = [funcCallStr ', ' toString(parameters{i})]; end funcCallStr = [funcCallStr ')'; end