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