classdef ResidualViscosity < handle
    properties
        D % Diff op approximating the gradient of the flux f(u)
        waveSpeed % Wave speed at each grid point, e.g f'(u). %TBD: Better naming?
        Cmax % Constant controlling relative amount of upwind dissipation
        Cres % Constant controling relative amount of upwind dissipation
        h % Length scale used for scaling the viscosity. Typically grid spacing.
        viscosity % Stores the computed viscosity.
        normalization % Function used to normalize the residual such that it is amplified in the
                      % shocks.

        % Convenice (for verification and plotting) TBD: Decide on if it should be kept.
        u_t % Stores the latest approximated time derivative of the solution.
        grad_f % Stores the latest approximated gradient of the flux
        residual % Stores the computed residual
    end

    methods
        %  TBD: Decide on how to treat waveSpeed. It would be nice to just pass a constant value without
        %       wrapping it in a function.
        function obj = ResidualViscosity(D, waveSpeed, Cmax, Cres, h, N, normalization)
            default_arg('normalization',@(v)norm(v-mean(v),inf));
            obj.D = D;
            obj.waveSpeed = waveSpeed;
            obj.h = h;
            obj.Cmax = Cmax;
            obj.Cres = Cres;
            obj.viscosity = zeros(N,1);
            obj.u_t = zeros(N,1);
            obj.grad_f = zeros(N,1);
            obj.residual = zeros(N,1);
            obj.normalization = normalization;
        end

        function obj = update(obj, v, dvdt)
            obj.u_t = dvdt;
            obj.grad_f = obj.D(v);
            obj.residual = obj.u_t + obj.grad_f;
            obj.viscosity = min(obj.Cmax*obj.h*abs(obj.waveSpeed(v)), obj.Cres*obj.h^2*abs(obj.residual)/obj.normalization(v));
        end

        function [residual, u_t, grad_f] = getResidual(obj)
            residual = obj.residual;
            u_t = obj.u_t;
            grad_f = obj.grad_f;
        end

        function viscosity = getViscosity(obj)
            viscosity = obj.viscosity;
        end
    end
end