view +time/SBPInTime.m @ 1037:2d7ba44340d0 feature/burgers1d

Pass scheme specific parameters as cell array. This will enabale constructDiffOps to be more general. In addition, allow for schemes returning function handles as diffOps, which is currently how non-linear schemes such as Burgers1d are implemented.
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Fri, 18 Jan 2019 09:02:02 +0100
parents 38173ea263ed
children 8894e9c49e40
line wrap: on
line source

classdef SBPInTime < time.Timestepper
    % The SBP in time method.
    % Implemented for v_t = A*v + f(t)
    %
    % Each "step" takes one block step and thus advances
    % k = k_local*(blockSize-1) in time.
    properties
        M     % System matrix
        L,U,P,Q % LU factorization of M
        A
        Et_r
        penalty
        f
        k_local % step size within a block
        k % Time size of a block  k/(blockSize-1) = k_local
        t
        v
        m
        n
        blockSize % number of points in each block
        order
        nodes
    end

    methods
        function obj = SBPInTime(A, f, k, t0, v0, TYPE, order, blockSize)

            default_arg('TYPE','gauss');

            if(strcmp(TYPE,'gauss'))
                default_arg('order',4)
                default_arg('blockSize',4)
            else
                default_arg('order', 8);
                default_arg('blockSize',time.SBPInTime.smallestBlockSize(order,TYPE));
            end

            obj.A = A;
            obj.f = f;
            obj.k_local = k/(blockSize-1);
            obj.k = k;
            obj.blockSize = blockSize;
            obj.t = t0;
            obj.m = length(v0);
            obj.n = 0;

            %==== Build the time discretization matrix =====%
            switch TYPE
                case 'equidistant'
                    ops = sbp.D2Standard(blockSize,{0,obj.k},order);
                case 'optimal'
                    ops = sbp.D1Nonequidistant(blockSize,{0,obj.k},order);
                case 'minimal'
                    ops = sbp.D1Nonequidistant(blockSize,{0,obj.k},order,'minimal');
                case 'gauss'
                    ops = sbp.D1Gauss(blockSize,{0,obj.k});
            end

            D1 = ops.D1;
            HI = ops.HI;
            e_l = ops.e_l;
            e_r = ops.e_r;
            obj.nodes = ops.x;

            Ix = speye(size(A));
            It = speye(blockSize,blockSize);

            obj.Et_r = kron(e_r,Ix);

            % Time derivative + penalty
            tau = 1;
            Mt = D1 + tau*HI*(e_l*e_l');

            % penalty to impose "data"
            penalty = tau*HI*e_l;
            obj.penalty = kron(penalty,Ix);

            Mx = kron(It,A);
            Mt = kron(Mt,Ix);
            obj.M = Mt - Mx;
            %==============================================%

            % LU factorization
            [obj.L,obj.U,obj.P,obj.Q] = lu(obj.M);

            % Pretend that the initial condition is the last level
            % of a previous step.
            obj.v =  1/(e_r'*e_r) * obj.Et_r * v0;

        end

        function [v,t] = getV(obj)
            v = obj.Et_r' * obj.v;
            t = obj.t;
        end

        function obj = step(obj)
            obj.v = time.sbp.sbpintime(obj.v, obj.t, obj.nodes,...
                              obj.penalty, obj.f, obj.blockSize,...
                              obj.Et_r,...
                              obj.L, obj.U, obj.P, obj.Q);
            obj.t = obj.t + obj.k;
            obj.n = obj.n + 1;
        end
    end


    methods(Static)
        function N = smallestBlockSize(order,TYPE)
            default_arg('TYPE','gauss')

            switch TYPE

                case 'equidistant'
                    switch order
                        case 2
                            N = 2;
                        case 4
                            N = 8;
                        case 6
                            N = 12;
                        case 8
                            N = 16;
                        case 10
                            N = 20;
                        case 12
                            N = 24;
                        otherwise
                            error('Operator does not exist');
                    end

                case 'optimal'

                    switch order
                        case 4
                            N = 8;
                        case 6
                            N = 12;
                        case 8
                            N = 16;
                        case 10
                            N = 20;
                        case 12
                            N = 24;
                        otherwise
                            error('Operator does not exist');
                    end

                case 'minimal'

                    switch order
                        case 4
                            N = 6;
                        case 6
                            N = 10;
                        case 8
                            N = 12;
                        case 10
                            N = 16;
                        case 12
                            N = 20;
                        otherwise
                            error('Operator does not exist');
                    end
                case 'gauss'
                    N = 4;
            end
        end
    end
end