feat: Pass cycleNumber to execute/output/collect calls (#89)

* Pass cycleNumber to execute/output/collect calls

* Update OBL examples for better comparison

* Update example doc

---------

Co-authored-by: alexbenedicto <[email protected]>

Author: av-novikov

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 )
 
 

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>
 

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 )
 
 

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,23 +103,23 @@
 			objectPath="ElementRegions/Region1"
 			fieldName="globalCompFraction"
 			component="0"
-			scale="0.7"/>
+			scale="0.276839"/>
 		<FieldSpecification
 			name="Region1InitCa"
 			initialCondition="1"
 			setNames="{ all }"
 			objectPath="ElementRegions/Region1"
 			fieldName="globalCompFraction"
 			component="1"
-			scale="0.0"/>
+			scale="1e-11"/>
 		<FieldSpecification
 			name="Region1InitC"
 			initialCondition="1"
 			setNames="{ all }"
 			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 }"

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 )

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" )

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" )

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 )
 
 

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 )
 
 

pygeos-tools/src/geos/pygeos_tools/solvers/Solver.py

@@ -701,34 +701,38 @@ 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
 
         Parameters
         ----------
             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:
         """Reinitialize Solver"""
         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
 
         Parameters
         ----------
             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

pygeos-tools/src/geos/pygeos_tools/solvers/WaveSolver.py

@@ -393,14 +393,16 @@ 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
 
         Parameters
         ----------
             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 )