--- a/+rv/constructDiffOps.m	Tue Mar 05 10:57:26 2019 +0100
+++ b/+rv/constructDiffOps.m	Wed Mar 06 09:45:52 2019 +0100
@@ -1,11 +1,11 @@
 function [D_rv, D_flux, DvDt, solutionPenalties, residualPenalties] = constructDiffOps(scheme, g, order, schemeParams, opSet, BCs)
     %% DiffOps for solution vector
-    [D, solutionPenalties] = constructTotalFluxDiffOp(scheme, g, order, schemeParams, opSet, BCs);
+    [D, solutionPenalties] = constructFluxDiffOp(scheme, g, order, schemeParams, opSet, BCs);
     D2 = constructSymmetricD2Operator(g, order, opSet);
-    D_rv = @(v,viscosity)(D(v) + D2(v, viscosity));
+    D_rv = @(v,Viscosity)(D(v) + D2(Viscosity)*v);
 
     %% DiffOps for residual viscosity
-    [D_flux, residualPenalties] = constructTotalFluxDiffOp(scheme, g, max(order-2,2), schemeParams, opSet, BCs);
+    [D_flux, residualPenalties] = constructFluxDiffOp(scheme, g, max(order-2,2), schemeParams, opSet, BCs);
     % DiffOp for flux in residual viscosity. Due to sign conventions of the implemnted schemes, we need to
     % change the sign.
     D_flux = @(v) -D_flux(v);
@@ -13,7 +13,7 @@
     DvDt = D;
 end
 
-function [D, penalties] = constructTotalFluxDiffOp(scheme, g, order, schemeParams, opSet, BCs)
+function [D, penalties] = constructFluxDiffOp(scheme, g, order, schemeParams, opSet, BCs)
     diffOp = scheme(g, order, schemeParams{:}, opSet);
     [D, penalties]  = addClosuresToDiffOp(diffOp, BCs);
 end
@@ -48,10 +48,8 @@
        I{i} = speye(m(i));
     end
 
-    % TBD: How is this generalized to a loop over dimensions or similar?
     switch g.D()
         case 1
-            
             e_r = ops{1}.e_r;
             e_l = ops{1}.e_l;
             Hi = ops{1}.HI;
@@ -59,9 +57,11 @@
             if isequal(opSet,@sbp.D1Upwind)
                 Dm = ops{1}.Dm;
                 Dp = ops{1}.Dp;
-                D2 = @(viscosity) Dm*spdiag(viscosity)*Dp-Hi*(B*spdiag(viscosity)*Dp);
+                M =  Dm - Hi*B;
+                D2 = @(Viscosity) M*Viscosity*Dp;
             else
-                D2 = @(viscosity)ops{1}.D2(viscosity);
+                % TODO: Fix closure for D2Variable
+                D2 = @(Viscosity)ops{1}.D2(Viscosity);
             end
         case 2
             % TODO: 
@@ -75,7 +75,8 @@
             Dp_x  = kron(ops{1}.Dp,I{2});
             H_x  = kron(ops{1}.HI,I{2});
             B_x = e_e*e_e' - e_w*e_w';
-            D2_x = @(viscosity) Dm_x*spdiag(viscosity)*Dp_x-H_x*(B_x*spdiag(viscosity)*Dp_x);
+            M_x = Dm_x-H_x*B_x;
+            D2_x = @(Viscosity) M_x*Viscosity*Dp_x;
 
             e_n = kron(I{1},ops{2}.e_r);
             e_s = kron(I{1},ops{2}.e_l);
@@ -83,10 +84,10 @@
             Dp_y  = kron(I{1},ops{2}.Dp);
             H_y = kron(I{1},ops{2}.HI);
             B_y = e_n*e_n' - e_s*e_s';
-            D2_y = @(viscosity) Dm_y*spdiag(viscosity)*Dp_y-H_y*(B_y*spdiag(viscosity)*Dp_y);
-            D2 = @(viscosity)D2_x(viscosity) + D2_y(viscosity);
+            M_y = Dm_y-H_y*B_y;
+            D2_y = @(Viscosity) M_y*Viscosity*Dp_y;
+            D2 = @(Viscosity) D2_x(Viscosity) + D2_y(Viscosity);
         otherwise
             error('3D not yet implemented')
     end
-    D2 = @(v, viscosity) D2(viscosity)*v;
 end
\ No newline at end of file