--- a/+scheme/Burgers1D.m	Mon Sep 03 14:50:27 2018 +0200
+++ b/+scheme/Burgers1D.m	Thu Sep 06 12:43:51 2018 +0200
@@ -4,6 +4,7 @@
         h % Grid spacing
         x % Grid
         order % Order accuracy for the approximation
+        params
 
         D % Non-stabalized scheme operator
         M % Derivative norm
@@ -13,21 +14,38 @@
         e_r
         d_l
         d_r
-        params % Parameters for the coefficient matrices
     end
 
     methods
-        function obj = Burgers1D(m, order, xlim, params)
-            [x, h] = util.get_grid(xlim{:},m);
-            ops = sbp.D2Variable(m, xlim, order);
+        function obj = Burgers1D(pde_form, operator_type, order, m, lim, params)
+            [x, h] = util.get_grid(lim{:},m);
+            default_arg('pde_form','skew-symmetric');
+            default_arg('operator_type','narrow');
+
+            switch operator_type
+                case 'narrow'
+                    ops = sbp.D2Variable(m, lim, order);
+                    D1 = ops.D1;
+                    D2 = ops.D2;  
+                otherwise
+                    error('Other operator types not yet supported', operator_type);
+            end
+
+            switch pde_form
+                case 'skew-symmetric'
+                    D = @(v, viscosity) -1/3*v.*D1*v - 1/3*D1*v.^2 + D2(obj.params.eps + viscosity)*v;
+                case 'conservative'
+                    D = @(v, viscosity) -1/2*D1*v.^2 + D2(obj.params.eps + viscosity)*v;
+            end
 
             obj.m = m;
             obj.h = h;
             obj.order = order;
             obj.x = x;
+            obj.params = params;
 
-            D1 = ops.D1;
-            obj.D = @(v)(-1/3*v.*D1*v - 1/3*D1*v.^2 + ops.D2(params.eps)*v);
+            %% TODO: Figure out how to evaluate viscosity as viscosity(v,t) here instead of parametrizing D on the viscosity.
+            obj.D = D;
             obj.M = ops.M;
             obj.H =  ops.H;
             obj.Hi = ops.HI;
@@ -35,7 +53,6 @@
             obj.e_r = ops.e_r;
             obj.d_l =  ops.d1_l;
             obj.d_r =  ops.d1_r;
-            obj.params = params;
         end
 
         % Closure functions return the opertors applied to the own doamin to close the boundary
@@ -53,7 +70,7 @@
                 % Stable robin-like boundary conditions ((u+-abs(u))*u/3 - eps*u_x)) with +- at left/right boundary
                 case {'R','robin'}
                     p = s*obj.Hi*e;
-                    closure = @(v) p*(e'*((v-s*abs(v))/3)*(e'*v) - e'*obj.params.eps*d'*v);
+                    closure = @(v, viscosity) p*(e'*((v-s*abs(v))/3)*(e'*v) - e'*(obj.params.eps + viscosity)*d'*v);
                     switch class(data)
                         case 'double'
                             penalty = s*p*data;

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+time/Rungekutta4RV.m	Thu Sep 06 12:43:51 2018 +0200
@@ -0,0 +1,54 @@
+classdef Rungekutta4RV < time.Timestepper
+    properties
+        F
+        k
+        t
+        v
+        m
+        n
+
+        % Additional members used for the RV update
+        RV
+    end
+
+
+    methods
+        function obj = Rungekutta4RV(F, k, t0, v0, RV)
+            obj.F = F;
+            obj.k = k;
+            obj.t = t0;
+            obj.v = v0;
+            obj.m = length(v0);
+            obj.n = 0;
+            obj.RV = RV;
+        end
+
+        function [v, t] = getV(obj)
+            v = obj.v;
+            t = obj.t;
+        end
+
+        function [residual, viscosity, t] = getRV(obj)
+            residual = obj.RV.getResidual();
+            viscosity = obj.RV.getViscosity();
+            t = obj.t;
+        end
+
+        function obj = step(obj)
+            v_prev = obj.v;
+            F = @(v,t) obj.F(v, t, obj.RV.getViscosity());
+            obj.v = time.rk4.rungekutta_4(obj.v, obj.t, obj.k, F);
+            obj.t = obj.t + obj.k;
+            obj.n = obj.n + 1;
+            obj.RV.update(obj.v, v_prev, obj.k);
+        end
+    end
+
+
+    methods (Static)
+        function k = getTimeStep(lambda)
+            k = rk4.get_rk4_time_step(lambda);
+        end
+    end
+
+end
\ No newline at end of file