diff --git a/docs/pygeos_tools_docs/example.rst b/docs/pygeos_tools_docs/example.rst
index a217eb8b..dda8cb9e 100644
--- a/docs/pygeos_tools_docs/example.rst
+++ b/docs/pygeos_tools_docs/example.rst
@@ -117,12 +117,12 @@ This block contains a ``forceDt`` attribute that will be used later to choose as
 
 The "outputs" event triggers the output of the vtk files. The attribute "timeFrequency" has the same value as "forceDt"
 so we can use the same timestep for the solver and the outputs.
-To start, we will set the time to 0.0 and trigger one output of the vtk files.
+To start, we will set the time to 0.0 and the cycle number to 0.
 
 .. code-block:: python
 
    time = 0.0
-   solver.outputVtk( time )
+   cycle = 0
 
 
 ------------------------------------------------------------------
@@ -135,9 +135,12 @@ Once done, the simulation is ended by calling the ``cleanup`` method.
 .. code-block:: python
 
    while time < solver.maxTime:
-       solver.execute( time )
-       solver.outputVtk( time )
+       solver.outputVtk( time, cycle )
+       solver.execute( time, cycle )
        time += solver.dt
+       cycle += 1
+
+   solver.outputVtk( time, cycle )
    solver.cleanup( time )
 
 
diff --git a/pygeos-tools/examples/reactiveCompositionalMultiphaseOBL_modeling/2ph_comp/input_file_adaptive.xml b/pygeos-tools/examples/reactiveCompositionalMultiphaseOBL_modeling/2ph_comp/input_file_adaptive.xml
index ec4e1ea2..a664970c 100644
--- a/pygeos-tools/examples/reactiveCompositionalMultiphaseOBL_modeling/2ph_comp/input_file_adaptive.xml
+++ b/pygeos-tools/examples/reactiveCompositionalMultiphaseOBL_modeling/2ph_comp/input_file_adaptive.xml
@@ -16,7 +16,9 @@
 			newtonTol="0.0001"
 			newtonMaxIter="25"/>
 		<LinearSolverParameters
-			directParallel="0"/>
+            solverType="fgmres"
+            preconditionerType="mgr"
+            krylovTol="1.0e-5"/>
 		</ReactiveCompositionalMultiphaseOBL>
 	</Solvers>
 
diff --git a/pygeos-tools/examples/reactiveCompositionalMultiphaseOBL_modeling/2ph_comp/main.py b/pygeos-tools/examples/reactiveCompositionalMultiphaseOBL_modeling/2ph_comp/main.py
index 2cf7cc7e..c33a4780 100644
--- a/pygeos-tools/examples/reactiveCompositionalMultiphaseOBL_modeling/2ph_comp/main.py
+++ b/pygeos-tools/examples/reactiveCompositionalMultiphaseOBL_modeling/2ph_comp/main.py
@@ -98,7 +98,7 @@ def run_darts_model( xml_name: str, darts_model=None ):
     cycle: int = 0
 
     solver.setDt( 86400.0 )
-    solver.outputVtk( time )
+
     while time < solver.maxTime:
         if time < 604800:
             solver.setDt( 86400.0 )
@@ -108,10 +108,14 @@ def run_darts_model( xml_name: str, darts_model=None ):
         if rank == 0:
             if solver.dt is not None:
                 print( f"time = {time:.3f}s, dt = {solver.getDt():.4f}, iter = {cycle + 1}" )
-        solver.execute( time )
+
+        solver.outputVtk( time, cycle )
+        solver.execute( time, cycle )
+
         time += solver.getDt()
-        solver.outputVtk( time )
         cycle += 1
+
+    solver.outputVtk( time, cycle )
     solver.cleanup( time )
 
 
diff --git a/pygeos-tools/examples/reactiveCompositionalMultiphaseOBL_modeling/carbonated_water/1d_setup.xml b/pygeos-tools/examples/reactiveCompositionalMultiphaseOBL_modeling/carbonated_water/1d_setup.xml
index a2e012dc..9e714abc 100644
--- a/pygeos-tools/examples/reactiveCompositionalMultiphaseOBL_modeling/carbonated_water/1d_setup.xml
+++ b/pygeos-tools/examples/reactiveCompositionalMultiphaseOBL_modeling/carbonated_water/1d_setup.xml
@@ -18,7 +18,9 @@
 			newtonTol="0.0001"
 			newtonMaxIter="25"/>
 		<LinearSolverParameters
-			directParallel="0"/>
+            solverType="fgmres"
+            preconditionerType="mgr"
+            krylovTol="1.0e-5"/>
 		</ReactiveCompositionalMultiphaseOBL>
 	</Solvers>
 
@@ -101,7 +103,7 @@
 			objectPath="ElementRegions/Region1"
 			fieldName="globalCompFraction"
 			component="0"
-			scale="0.7"/>
+			scale="0.276839"/>
 		<FieldSpecification
 			name="Region1InitCa"
 			initialCondition="1"
@@ -109,7 +111,7 @@
 			objectPath="ElementRegions/Region1"
 			fieldName="globalCompFraction"
 			component="1"
-			scale="0.0"/>
+			scale="1e-11"/>
 		<FieldSpecification
 			name="Region1InitC"
 			initialCondition="1"
@@ -117,7 +119,7 @@
 			objectPath="ElementRegions/Region1"
 			fieldName="globalCompFraction"
 			component="2"
-			scale="0.0"/>
+			scale="1e-11"/>
 		<FieldSpecification
 			name="Region1InitO"
 			initialCondition="1"
@@ -139,7 +141,7 @@
 			name="inj1Pressure"
 			objectPath="ElementRegions/Region1"
 			fieldName="pressure"
-			scale="10010000.0"
+			scale="10005000.0"
 			setNames="{ inj1 }"/>
 		<FieldSpecification
 			name="inj1Calcite"
@@ -219,21 +221,21 @@
 			objectPath="ElementRegions/Region1"
 			fieldName="globalCompFraction"
 			component="0"
-			scale="0.7"/>
+			scale="0.276839"/>
 		<FieldSpecification
 			name="prd1Ca"
 			setNames="{ prd1 }"
 			objectPath="ElementRegions/Region1"
 			fieldName="globalCompFraction"
 			component="1"
-			scale="0.0"/>
+			scale="1e-11"/>
 		<FieldSpecification
 			name="prd1C"
 			setNames="{ prd1 }"
 			objectPath="ElementRegions/Region1"
 			fieldName="globalCompFraction"
 			component="2"
-			scale="0.0"/>
+			scale="1e-11"/>
 		<FieldSpecification
 			name="prd1O"
 			setNames="{ prd1 }"
diff --git a/pygeos-tools/examples/reactiveCompositionalMultiphaseOBL_modeling/carbonated_water/main.py b/pygeos-tools/examples/reactiveCompositionalMultiphaseOBL_modeling/carbonated_water/main.py
index 58a677ba..4ebb62be 100644
--- a/pygeos-tools/examples/reactiveCompositionalMultiphaseOBL_modeling/carbonated_water/main.py
+++ b/pygeos-tools/examples/reactiveCompositionalMultiphaseOBL_modeling/carbonated_water/main.py
@@ -51,7 +51,6 @@ def run_darts_model( domain: str, xml_name: str, darts_model=None ):
     cycle: int = 0
     solver.setDt( 8.64 )
 
-    solver.outputVtk( time )
     while time < solver.maxTime:
         # choose new timestep
         if domain == '1D':
@@ -84,16 +83,20 @@ def run_darts_model( domain: str, xml_name: str, darts_model=None ):
                 solver.setDt( 100.0 )
         if rank == 0:
             print( f"time = {time:.3f}s, dt = {solver.getDt():.4f}, step = {cycle + 1}" )
-        # run simulation
-        solver.execute( time )
-        time += solver.getDt()
+
         if cycle % 5 == 0:
-            solver.outputVtk( time )
+            solver.outputVtk( time, cycle )
+
+        solver.execute( time, cycle )
+
+        time += solver.getDt()
         cycle += 1
+
+    solver.outputVtk( time, cycle )
     solver.cleanup( time )
 
 
 if __name__ == "__main__":
     darts_model = Model()
-    # run_darts_model( domain='1D', xml_name="1d_setup.xml", darts_model=darts_model )
-    run_darts_model( domain='2D', xml_name="2d_setup.xml", darts_model=darts_model )
+    run_darts_model( domain='1D', xml_name="1d_setup.xml", darts_model=darts_model )
+    # run_darts_model( domain='2D', xml_name="2d_setup.xml", darts_model=darts_model )
diff --git a/pygeos-tools/examples/solvers/acoustic_modeling.py b/pygeos-tools/examples/solvers/acoustic_modeling.py
index 07827260..92732003 100644
--- a/pygeos-tools/examples/solvers/acoustic_modeling.py
+++ b/pygeos-tools/examples/solvers/acoustic_modeling.py
@@ -149,12 +149,17 @@ def main():
         while time < solver.maxTime:
             if rank == 0 and cycle % 100 == 0:
                 print( f"time = {time:.3f}s, dt= {solver.dt:.4f}, iter = {cycle + 1}" )
-            solver.execute( time )
+
             if cycle % 50 == 0:
-                solver.outputVtk( time )
+                solver.outputVtk( time, cycle )
+
+            solver.execute( time, cycle )
+
             time += solver.dt
             cycle += 1
 
+        solver.outputVtk( time, cycle )
+
         shot.flag = "Done"
         if rank == 0:
             print( f"Shot {shot.id} done" )
diff --git a/pygeos-tools/examples/solvers/elastic_modeling.py b/pygeos-tools/examples/solvers/elastic_modeling.py
index 10de4918..aa489623 100644
--- a/pygeos-tools/examples/solvers/elastic_modeling.py
+++ b/pygeos-tools/examples/solvers/elastic_modeling.py
@@ -82,12 +82,17 @@ def main():
         while t < solver.maxTime:
             if rank == 0 and cycle % 100 == 0:
                 print( f"time = {t:.3f}s, dt = {solver.dt:.4f}, iter = {cycle + 1}" )
-            solver.execute( t )
+
             if cycle % 100 == 0:
-                solver.outputVtk( t )
+                solver.outputVtk( t, cycle )
+
+            solver.execute( t, cycle )
+
             t += solver.dt
             cycle += 1
 
+        solver.outputVtk( t, cycle )
+
         shot.flag = "Done"
         if rank == 0:
             print( f"Shot {shot.id} done" )
diff --git a/pygeos-tools/examples/solvers/geomechanics_modeling.py b/pygeos-tools/examples/solvers/geomechanics_modeling.py
index a3a288d7..14516ae1 100644
--- a/pygeos-tools/examples/solvers/geomechanics_modeling.py
+++ b/pygeos-tools/examples/solvers/geomechanics_modeling.py
@@ -53,15 +53,16 @@ def main():
     time: float = 0.0
     cycle: int = 0
 
-    solver.outputVtk( time )
     while time < solver.maxTime:
         if rank == 0:
             if solver.dt is not None:
                 print( f"time = {time:.3f}s, dt = {solver.dt:.4f}, iter = {cycle + 1}" )
-        solver.execute( time )
-        solver.outputVtk( time )
+        solver.outputVtk( time, cycle )
+        solver.execute( time, cycle )
         time += solver.dt
         cycle += 1
+
+    solver.outputVtk( time, cycle )
     solver.cleanup( time )
 
 
diff --git a/pygeos-tools/examples/solvers/reservoir_modeling.py b/pygeos-tools/examples/solvers/reservoir_modeling.py
index 1e075cf8..885cde06 100644
--- a/pygeos-tools/examples/solvers/reservoir_modeling.py
+++ b/pygeos-tools/examples/solvers/reservoir_modeling.py
@@ -53,17 +53,18 @@ def main():
     time: float = 0.0
     cycle: int = 0
 
-    solver.outputVtk( time )
     while time < solver.maxTime:
         if rank == 0:
             if solver.dt is not None:
                 print( f"time = {time:.3f}s, dt = {solver.dt:.4f}, iter = {cycle + 1}" )
-        solver.execute( time )
-        solver.outputVtk( time )
+        solver.outputVtk( time, cycle )
+        solver.execute( time, cycle )
         pressure = solver.getPressures()
         print( pressure )
         time += solver.dt
         cycle += 1
+
+    solver.outputVtk( time, cycle )
     solver.cleanup( time )
 
 
diff --git a/pygeos-tools/src/geos/pygeos_tools/solvers/Solver.py b/pygeos-tools/src/geos/pygeos_tools/solvers/Solver.py
index d518496b..e562836d 100644
--- a/pygeos-tools/src/geos/pygeos_tools/solvers/Solver.py
+++ b/pygeos-tools/src/geos/pygeos_tools/solvers/Solver.py
@@ -701,7 +701,7 @@ def cleanup( self: Self, time: float ) -> None:
         self.solver.cleanup( time )
 
     @required_attributes( "solver" )
-    def execute( self: Self, time: float ) -> None:
+    def execute( self: Self, time: float, cycleNumber: int ) -> None:
         """
         Do one solver iteration
 
@@ -709,8 +709,10 @@ def execute( self: Self, time: float ) -> None:
         ----------
             time : float
                 Current time of simulation
+            cycleNumber : int
+                Current cycle number
         """
-        self.solver.execute( time, self.dt )
+        self.solver.execute( time, self.dt, cycleNumber )
 
     @required_attributes( "solver" )
     def reinitSolver( self: Self ) -> None:
@@ -718,7 +720,7 @@ def reinitSolver( self: Self ) -> None:
         self.solver.reinit()
 
     @required_attributes( "vtkOutputs" )
-    def outputVtk( self: Self, time: float ) -> None:
+    def outputVtk( self: Self, time: float, cycleNumber: int ) -> None:
         """
         Trigger the VTK output
 
@@ -726,9 +728,11 @@ def outputVtk( self: Self, time: float ) -> None:
         ----------
             time : float
                 Current time of simulation
+            cycleNumber : int
+                Current cycle number
         """
         for vtkOutput in self.vtkOutputs:
-            vtkOutput.output( time, self.dt )
+            vtkOutput.output( time, self.dt, cycleNumber )
 
     """
     Update methods when initializing or reinitializing the solver
diff --git a/pygeos-tools/src/geos/pygeos_tools/solvers/WaveSolver.py b/pygeos-tools/src/geos/pygeos_tools/solvers/WaveSolver.py
index 9daf59c0..2e05d06f 100644
--- a/pygeos-tools/src/geos/pygeos_tools/solvers/WaveSolver.py
+++ b/pygeos-tools/src/geos/pygeos_tools/solvers/WaveSolver.py
@@ -393,7 +393,7 @@ def filterSource( self: Self, fmax: Union[ str, float ] ) -> None:
         self.setGeosWrapperValueByTargetKey( "Events/minTime", self.minTimeSim )
         self.sourceValue = np.real( y[ max( i1 - d, 0 ):min( i4 + d, n ), : ] )
 
-    def outputWaveField( self: Self, time: float ) -> None:
+    def outputWaveField( self: Self, time: float, cycleNumber: int ) -> None:
         """
         Trigger the wavefield output
 
@@ -401,6 +401,8 @@ def outputWaveField( self: Self, time: float ) -> None:
         ----------
             time : float
                 Current time of simulation
+            cycleNumber : int
+                Current cycle number
         """
-        self.collections[ 0 ].collect( time, self.dt )
-        self.hdf5Outputs[ 0 ].output( time, self.dt )
+        self.collections[ 0 ].collect( time, self.dt, cycleNumber )
+        self.hdf5Outputs[ 0 ].output( time, self.dt, cycleNumber )