Mercurial > repos > public > sbplib
view +util/calc_borrowing.m @ 1031:2ef20d00b386 feature/advectionRV
For easier comparison, return both the first order and residual viscosity when evaluating the residual. Add the first order and residual viscosity to the state of the RungekuttaRV time steppers
author | Vidar Stiernström <vidar.stiernstrom@it.uu.se> |
---|---|
date | Thu, 17 Jan 2019 10:25:06 +0100 |
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