diff --git a/docs/pygeos_tools_docs/example.rst b/docs/pygeos_tools_docs/example.rst
index a217eb8be..dda8cb9e8 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 ec4e1ea22..a664970c3 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 2cf7cc7ef..c33a47802 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 a2e012dc3..9e714abc7 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 58a677ba3..4ebb62bef 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 07827260c..927320036 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 10de49187..aa489623e 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 a3a288d7e..14516ae10 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 1e075cf87..885cde06c 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 d518496bf..e562836de 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 9daf59c08..2e05d06f7 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 )