view +time/Rungekutta.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 4e5e53d6336c
children
line wrap: on
line source

classdef Rungekutta < time.Timestepper
    properties
        F       % RHS of the ODE
        k       % Time step
        t       % Time point
        v       % Solution vector
        n       % Time level
        scheme  % The scheme used for the time stepping, e.g rk4, rk6 etc.
    end


    methods
        % Timesteps v_t = F(v,t), using RK with specfied order from t = t0 with
        % timestep k and initial conditions v = v0
        function obj = Rungekutta(F, k, t0, v0, order)
            default_arg('order',4);
            obj.F = F;
            obj.k = k;
            obj.t = t0;
            obj.v = v0;
            obj.n = 0;
            % TBD: Order 4 is also implemented in the butcher tableau, but the rungekutta_4.m implementation
            % might be slightly more efficient. Need to do some profiling before deciding whether or not to keep it.
            if (order == 4)
                obj.scheme = @time.rk.rungekutta_4;
            else
                % Extract the coefficients for the specified order
                % used for the RK updates from the Butcher tableua.
                [s,a,b,c] = time.rk.butcherTableau(order);
                coeffs = struct('s',s,'a',a,'b',b,'c',c);
                obj.scheme = @(v,t,dt,F) time.rk.rungekutta(v, t , dt, F, coeffs);
            end
        end

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

        function obj = step(obj)
            obj.v = obj.scheme(obj.v, obj.t, obj.k, obj.F);
            obj.t = obj.t + obj.k;
            obj.n = obj.n + 1;
        end
    end
end