--- a/+scheme/Burgers1D.m	Tue Sep 11 16:29:00 2018 +0200
+++ b/+scheme/Burgers1D.m	Thu Sep 13 14:21:47 2018 +0200
@@ -6,7 +6,6 @@
         params
 
         D % Non-stabalized scheme operator
-        M % Derivative norm
         H % Discrete norm
         Hi % Norm inverse
         e_l
@@ -16,44 +15,63 @@
     end
 
     methods
-        function obj = Burgers1D(grid, pde_form, operator_type, order, params)
+        function obj = Burgers1D(grid, pde_form, operator_type, order, dissipation, params)
             assert(grid.D == 1);
             assert(grid.size() == length(params.eps));
             m = grid.size();
-            h = grid.scaling();
             lim = grid.lim{1}; % Ugly, and only applicable for cartesian grids.
             default_arg('pde_form','skew-symmetric');
             default_arg('operator_type','narrow');
+            default_arg('dissipation','on');
 
             switch operator_type
                 case 'narrow'
                     ops = sbp.D4Variable(m, lim, order);
                     D1 = ops.D1;
-                    D2 = ops.D2;  
+                    D2 = ops.D2;
+                    DissipationOp = -1*sbp.dissipationOperator(m, order, ops.HI);
+                    d_l = ops.d1_l';
+                    d_r = ops.d1_r';
+                case 'upwind'
+                    ops = sbp.D1Upwind(m, lim, order);
+                    D1 = (ops.Dp + ops.Dm)/2;
+                    %D2eps = @(eps) ops.Dp*diag(eps)*ops.Dm;
+                    %D2eps = @(eps) ops.Dm*diag(eps)*ops.Dp;
+                    D2 = @(eps) (ops.Dp*diag(eps)*ops.Dm + ops.Dm*diag(eps)*ops.Dp)/2;
+                    DissipationOp = (ops.Dp-ops.Dm)/2;
+                    d_l = D1;
+                    d_r = D1;
                 otherwise
                     error('Other operator types not yet supported', operator_type);
             end
 
+            %% TODO: Figure out how to evaluate viscosity as viscosity(v,t) here instead of parametrizing D on the viscosity.
             switch pde_form
                 case 'skew-symmetric'
-                    D = @(v, viscosity) -1/3*v.*D1*v - 1/3*D1*v.^2 + D2(params.eps + viscosity)*v;
+                    if (strcmp(dissipation,'on'))
+                        D = @(v, viscosity) - 1/3*D1*v.^2 + (-1/3*v.*D1 + D2(params.eps + viscosity) + max(abs(v))*DissipationOp)*v;
+                    else
+                        D = @(v, viscosity) - 1/3*D1*v.^2 + (-1/3*v.*D1 + D2(params.eps + viscosity))*v;
+                    end
                 case 'conservative'
-                    D = @(v, viscosity) -1/2*D1*v.^2 + D2(params.eps + viscosity)*v;
+                    if (strcmp(dissipation,'on'))
+                        D = @(v, viscosity) -1/2*D1*v.^2 + (D2(params.eps + viscosity) + max(abs(v))*DissipationOp)*v;
+                    else
+                        D = @(v, viscosity) -1/2*D1*v.^2 + D2(params.eps + viscosity)*v;
+                    end
             end
 
             obj.grid = grid;
             obj.order = order;
             obj.params = params;
 
-            %% 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;
             obj.e_l = ops.e_l;
             obj.e_r = ops.e_r;
-            obj.d_l =  ops.d1_l;
-            obj.d_r =  ops.d1_r;
+            obj.d_l =  d_l;
+            obj.d_r =  d_r;
         end
 
         % Closure functions return the opertors applied to the own doamin to close the boundary
@@ -71,7 +89,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, viscosity) p*(e'*((v-s*abs(v))/3)*(e'*v) - e'*(obj.params.eps + viscosity)*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;