Mercurial > repos > public > sbplib
view +util/calc_borrowing.m @ 774:66eb4a2bbb72 feature/grids
Remove default scaling of the system.
The scaling doens't seem to help actual solutions. One example that fails in the flexural code.
With large timesteps the solutions seems to blow up. One particular example is profilePresentation
on the tdb_presentation_figures branch with k = 0.0005
author | Jonatan Werpers <jonatan@werpers.com> |
---|---|
date | Wed, 18 Jul 2018 15:42:52 -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