diff --git a/lib/solvers/LowestOrderLocalMg.m b/lib/solvers/LowestOrderLocalMg.m
index cde8fb1d573b76a5d2de1c0927f0222c1637b3a9..af461468a373e915b4e6311ad4c124673e485c52 100644
--- a/lib/solvers/LowestOrderLocalMg.m
+++ b/lib/solvers/LowestOrderLocalMg.m
@@ -25,6 +25,7 @@ classdef LowestOrderLocalMg < MGSolver
         freeVertices
         freeVerticesOld
         changedPatches
+        otherSolver
     end
     
     %% event data
@@ -55,9 +56,11 @@ classdef LowestOrderLocalMg < MGSolver
             obj.blf = blf;
             
             obj.listenerHandle = mesh.listener('IsAboutToRefine', @obj.getChangedPatches);
+            obj.otherSolver = OptimalVcycleMultigridSolver(P1JacobiSmoother(fes, blf, P));
         end
 
         function setupSystemMatrix(obj, A)
+            obj.otherSolver.setupSystemMatrix(A);
             obj.nLevels = obj.nLevels + 1;
             
             L = obj.nLevels;
@@ -76,6 +79,18 @@ classdef LowestOrderLocalMg < MGSolver
 
         function setupRhs(obj, b, varargin)
             setupRhs@MGSolver(obj, b, varargin{:});
+            obj.otherSolver.setupRhs(b, varargin{:});
+        end
+
+        function step(obj)
+            obj.otherSolver.applyStoppingCriterion();
+            obj.otherSolver.step();
+            step@MGSolver(obj);
+
+            assert(all(abs(obj.Cresidual - obj.otherSolver.Cresidual) < 1e-8, 'all'))
+            assert(all(abs(obj.residualCNorm - obj.otherSolver.residualCNorm.^2) < 1e-8, 'all'))
+            assert(all(abs(obj.iterationCount - obj.otherSolver.iterationCount) < 1e-8, 'all'))
+            assert(all(abs(obj.x - obj.otherSolver.x) < 1e-8, 'all'))
         end
 
         % Geometric MultiGrid
diff --git a/lib/solvers/chooseIterativeSolver.m b/lib/solvers/chooseIterativeSolver.m
index fa30db22d2f6701c3fc1c43c78c4b31d43c63937..02d735c3b8fb17f82f0f0c2bae9c2a5658f8383a 100644
--- a/lib/solvers/chooseIterativeSolver.m
+++ b/lib/solvers/chooseIterativeSolver.m
@@ -53,13 +53,13 @@ switch class
         switch variant
             case {"", "lowOrderVcycle"}
                 if order == 1
-                    solver = OptimalVcycleMultigridSolver(P1JacobiSmoother(fes, blf, P));
+                    solver = LowestOrderLocalMg(fes, blf, P);
                 else
                     solver = LocalMgLowOrderVcycle(fes, blf);
                 end
             case "highOrderVcycle"
                 if order == 1
-                    solver = OptimalVcycleMultigridSolver(P1JacobiSmoother(fes, blf, P));
+                    solver = LowestOrderLocalMg(fes, blf, P);
                 else
                     solver = LocalMgHighOrderVcycle(fes, blf, P);
                 end