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.
        boundaryIndices % indices of the boundary of the domain
    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(g, Df, waveSpeed, Cmax, Cres, h, normalization)
            default_arg('normalization',@(v)abs(obj.minmaxDiffNeighborhood1d(v)-norm(v-mean(v),inf)));
            obj.Df = Df;
            obj.waveSpeed = waveSpeed;
            obj.h = h;
            obj.Cmax = Cmax;
            obj.Cres = Cres;
            obj.normalization = normalization;
            obj.boundaryIndices = obj.getBoundaryIndices(g);
        end

        function viscosity = evaluateViscosity(obj, v, dvdt)
            viscosity = min(obj.Cmax*obj.h*abs(obj.waveSpeed(v)), obj.Cres*obj.h^2*abs(dvdt + obj.Df(v))./obj.normalization(v));
            % TODO: If the boundary conditions are implemented correctly, this might not be needed.
            viscosity(obj.boundaryIndices) = 0;
        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);
            % TODO: If the boundary conditions are implemented correctly, this might not be needed.
            viscosity(obj.boundaryIndices) = 0;
        end
    end
    methods (Static)

        function boundaryind = getBoundaryIndices(g)
            ind = grid.funcToMatrix(g, 1:g.N()); 
            switch g.D()
            case 1
                boundaryind = [ind(1) ind(end)];
            case 2
                boundaryind = [ind(1,:) ind(end,:), ind(:,1)' ind(:,end)'];
            case 3
                error;
            end
        end

        function minmaxDiff = minmaxDiffNeighborhood1d(u)
            umax = movmax(u,3);
            umin = movmin(u,3);
            minmaxDiff = umax - umin;
        end
        function minmaxDiff = minmaxDiffNeighborhood2d(g, u)
            m = g.m();
            uMatrix = grid.funcToMatrix(g,u);
            umax = max(movmax(uMatrix,3,1),movmax(uMatrix,3,2));
            umin = min(movmin(uMatrix,3,1),movmin(uMatrix,3,2));
            minmaxDiff = umax - umin;
            minmaxDiff = reshape(minmaxDiff',m(1)*m(2),1);
        end
    end
end