classdef ResidualViscosity < handle
    properties
        Df % 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.
        normalization % Function used to normalize the residual such that it is amplified in the
                      % shocks.
    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(Df, waveSpeed, Cmax, Cres, h, normalization)
            default_arg('normalization',@(v)norm(v-mean(v),inf));
            obj.Df = Df;
            obj.waveSpeed = waveSpeed;
            obj.h = h;
            obj.Cmax = Cmax;
            obj.Cres = Cres;
            obj.normalization = normalization;
        end

        function [viscosity, Df, firstOrderViscosity, residualViscosity] = evaluate(obj, v, dvdt)
            Df = obj.Df(v);
            firstOrderViscosity = obj.Cmax*obj.h*abs(obj.waveSpeed(v));
            residualViscosity = obj.Cres*obj.h^2*abs(dvdt + Df)/obj.normalization(v);
            viscosity = min(firstOrderViscosity, residualViscosity);
            %viscosity = obj.Cmax*obj.h*abs(obj.waveSpeed(v));
        end
    end
end