--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+rv/ResidualViscosity.m	Thu Jan 17 10:44:12 2019 +0100
@@ -0,0 +1,33 @@
+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
\ No newline at end of file