--- a/+rv/constructDiffOps.m	Mon Jun 24 17:44:21 2019 +0200
+++ b/+rv/constructDiffOps.m	Tue Jun 25 13:24:01 2019 +0200
@@ -1,12 +1,20 @@
-function [D_rv, D_flux, DvDt, solutionPenalties, residualPenalties] = constructDiffOps(scheme, g, order, schemeParams, opSet, BCs)
+function [D_rv, D_flux, DvDt, solutionPenalties, residualPenalties] = constructDiffOps(scheme, g, schemeOrder, residualOrder, schemeParams, opSet, BCs)
     %% DiffOps for solution vector
-    [D, solutionPenalties] = constructFluxDiffOp(scheme, g, order, schemeParams, opSet, BCs);
-    D2 = constructSymmetricD2Operator(g, order, opSet);
+    [D, solutionPenalties] = constructFluxDiffOp(scheme, g, schemeOrder, schemeParams, opSet, BCs);
+    D2 = constructSymmetricD2Operator(g, schemeOrder, opSet);
     D_rv = @(v,Viscosity)(D(v) + D2(Viscosity)*v);
 
     %% DiffOps for residual viscosity
-    [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
+    [D_flux, residualPenalties] = constructFluxDiffOp(scheme, g, residualOrder, schemeParams, opSet, BCs);
+    % TODO: Construct D_flux without closures when using bdfs.
+    % diffOp = scheme(g, residualOrder, schemeParams{:}, opSet);
+    % if ~isa(diffOp.D, 'function_handle')
+    %     D_flux = @(v) diffOp.D*v;
+    % else
+    %     D_flux = diffOp.D;
+    % end
+
+    % DiffOp for flux in residual viscosity. Due to sign conventions of the implemented schemes, we need to
     % change the sign.
     D_flux = @(v) -D_flux(v);
     % DiffOp for time derivative in residual viscosity
@@ -37,8 +45,6 @@
 end
 
 function D2 = constructSymmetricD2Operator(g, order, opSet)
-
-
     m = g.size();
     ops = cell(g.D(),1);
     I = cell(g.D(),1);
@@ -61,7 +67,8 @@
                 D2 = @(Viscosity) M*Viscosity*Dp;
             else
                 % TODO: Fix closure for D2Variable
-                D2 = @(Viscosity)ops{1}.D2(Viscosity);
+                % TODO: Fix Viscosity not being vector
+                D2 = @(Viscosity)ops{1}.D2(diag(Viscosity));
             end
         case 2
             % TODO: