classdef ResidualViscosity < handle
    properties
        grid % grid
        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 and suppressed elsewhere.
        Mres % Coefficients for the residual viscosity
        Mfirst % Coefficients for the first order viscosity
        fRes % Function handle for computing the 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(grid, Df, waveSpeed, Cmax, Cres, h, normalization, postProcess)
            default_arg('normalization',@(v)abs(obj.minmaxDiffNeighborhood1d(v)-norm(v-mean(v),inf)));
            default_arg('postProcess','maximum neighbors');
            obj.Df = Df;
            obj.waveSpeed = waveSpeed;
            obj.h = h;
            obj.Cmax = Cmax;
            
            obj.Cres = Cres;
            obj.normalization = normalization;
            obj.grid = grid;
            obj.Mres = obj.Cres*obj.h^2;
            obj.Mfirst = obj.Cmax*obj.h;
            switch postProcess
                case 'none'
                    obj.fRes = @(v,dvdt) obj.Mres*abs(dvdt + obj.Df(v))./obj.normalization(v);
                case 'maximum neighbors'
                    switch grid.D()
                        case 1
                            obj.fRes = @(v,dvdt) movmax(obj.Mres*abs(dvdt + obj.Df(v))./obj.normalization(v),3)
                        case 2
                            obj.fRes = @obj.maxResidualNeighbors2d
                    end
            end
        end

        function viscosity = evaluateViscosity(obj, v, dvdt)
            firstOrderViscosity = obj.Mfirst*abs(obj.waveSpeed(v));
            viscosity = min(firstOrderViscosity, obj.fRes(v,dvdt));
        end

        function [viscosity, Df, firstOrderViscosity, residualViscosity] = evaluate(obj, v, dvdt)
            Df = obj.Df(v);
            firstOrderViscosity = obj.Mfirst*abs(obj.waveSpeed(v));
            residualViscosity = obj.fRes(v,dvdt);
            viscosity = min(firstOrderViscosity, residualViscosity);
        end

        function res = maxResidualNeighbors2d(obj,v,dvdt)
            res = obj.Mres*abs(dvdt + obj.Df(v))./obj.normalization(v);
            resMat = grid.funcToMatrix(obj.grid,res);
            res = reshape(max(movmax(resMat,3,1),movmax(resMat,3,2))',obj.grid.N(),1);
        end
    end
    methods (Static)
        function minmaxDiff = minmaxDiffNeighborhood1d(u)
            umax = movmax(u,3);
            umin = movmin(u,3);
            minmaxDiff = umax - umin;
        end
        
        function minmaxDiff = minmaxDiffNeighborhood2d(grid, u)
            uMatrix = grid.funcToMatrix(grid,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',g.N(),1);
        end
    end
end