classdef ResidualViscosity < handle
    properties
        g % 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
        % TODO: pass opt struct with waveSpeed, normalization etc.
        %  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, postProcess)
            %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.g = g;
            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 {'filt', 'filter'}
                    % TBD: Keep?
                    order = 
                    F = obj.shapiroFilter(obj.g, order);
                    obj.Mres = F*obj.Mres;
                    obj.fRes = @(v,dvdt) obj.Mres*abs(dvdt + obj.Df(v))./obj.normalization(v);
                case {'max', 'maximum neighbors'}
                    switch g.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)
            viscosity = min(obj.Mfirst*abs(obj.waveSpeed(v)), 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.g,res);
            res = reshape(max(movmax(resMat,3,1),movmax(resMat,3,2))',obj.g.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(g, u)
            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',g.N(),1);
        end
        function F = shapiroFilter(g, order)
            switch order
                case 2
                    F = spdiags(repmat([1 2 1],g.m(1),1), -1:1, g.m(1), g.m(1));
                case 4
                    F = spdiags(repmat([-1 4 10 4 -1],g.m(1),1), -2:2, g.m(1), g.m(1));
                case 6
                    F = spdiags(repmat([1 -6 15 44 15 -6 1],g.m(1),1), -3:3, g.m(1), g.m(1));
                case 8
                    F = spdiags(repmat([-1 8 -28 56 186 56 -28 8 -1],g.m(1),1), -4:4, g.m(1), g.m(1));
            end
            F = 1/2^order * F;
            % Fx = spdiags(repmat(1/(mid+2)*[1 mid 1],g.m(1),1), -1:1, g.m(1), g.m(1));
            % Fx(1,:) = 0;
            % Fx(1,1:2) = 1/(mid+1)*[mid 1];
            % Fx(end,:) = 0;
            % Fx(end,end-1:end) = 1/(mid+1)*[1 mid];
            % Fy = spdiags(repmat(1/(mid+2)*[1 mid 1],g.m(2),1), -1:1, g.m(2), g.m(2));
            % Fy(1,:) = 0;
            % Fy(1,1:2) = 1/(mid+1)*[mid 1];
            % Fy(end,:) = 0;
            % Fy(end,end-1:end) = 1/(mid+1)*[1 mid];
            % Ix = speye(g.m(1));
            % Iy = speye(g.m(1));
            % F = kron(Fx, Iy) + kron(Ix, Fy);
        end
    end

end