feat: Python interfaces to geos from Makutu repo (refactoring)

docs/pygeos_tools_docs/Example/reservoir.rst

@@ -15,7 +15,7 @@ The example python script for this documentation is the following:
 
 .. code-block:: console
 
-   pygeos-tools/src/solvers_examples/reservoir_modeling.py
+   pygeos-tools/examples/solvers/reservoir_modeling.py
 
 
 ------------------------------------------------------------------
@@ -152,7 +152,7 @@ Using the same python used to build your GEOS installation with, run this comman
 
 .. code-block:: console
 
-   python pygeos-tools/src/solvers_examples/reservoir_modeling.py
+   python pygeos-tools/examples/solvers/reservoir_modeling.py
           --xml /path/to/your/GEOS/src/inputFiles/compositionalMultiphaseFlow/2ph_cap_1d_ihu.xml
 
 

pygeos-tools/examples/obl/2ph_comp/input_file_adaptative.xml

@@ -0,0 +1,213 @@
+<?xml version="1.0"?>
+<Problem>
+	<Solvers>
+		<ReactiveCompositionalMultiphaseOBL
+			name="compflow"
+			logLevel="1"
+			discretization="fluidTPFA"
+			targetRegions="{Region1}"
+			enableEnergyBalance="0"
+			maxCompFractionChange="1"
+			numComponents="3"
+			numPhases="2"
+			transMultExp="1">
+		<NonlinearSolverParameters
+			timeStepDecreaseFactor="0.5"
+			newtonTol="0.0001"
+			newtonMaxIter="25"/>
+		<LinearSolverParameters
+			directParallel="0"/>
+		</ReactiveCompositionalMultiphaseOBL>
+	</Solvers>
+
+	<Mesh>
+		<InternalMesh
+			name="mesh1"
+			elementTypes="{C3D8}"
+			xCoords="{0.0, 500.0}"
+			yCoords="{0.0, 500.0}"
+			zCoords="{0.0, 20.0}"
+			nx="{25}"
+			ny="{25}"
+			nz="{1}"
+			cellBlockNames="{block1}"/>
+	</Mesh>
+
+	<Geometry>
+		<Box
+			name="inj1"
+			xMin="{-0.2, -0.2, -0.2}"
+			xMax="{20.2, 20.2, 20.2}"/>
+		<Box
+			name="prd1"
+			xMin="{479.8, 479.8, -0.1999999999999993}"
+			xMax="{500.2, 500.2, 20.2}"/>
+	</Geometry>
+
+	<Constitutive>
+		<CompressibleSolidConstantPermeability
+			name="rock"
+			solidModelName="nullSolid"
+			porosityModelName="rockPorosity"
+			permeabilityModelName="rockPerm"/>
+		<NullModel
+			name="nullSolid"/>
+		<PressurePorosity
+			name="rockPorosity"
+			defaultReferencePorosity="0.3"
+			referencePressure="5000000.0"
+			compressibility="1e-09"/>
+ 		<ConstantPermeability
+			name="rockPerm"
+			permeabilityComponents="{9.869e-14, 9.869e-14, 9.869e-14}"/>
+		<OBLFluid
+			name="fluid"
+			interpolatorMode="adaptive"
+			interpolatorType="multilinear"/>
+	</Constitutive>
+
+	<Events
+		maxTime="31536000.0">
+		<!--PeriodicEvent
+			name="outputs"
+			timeFrequency="2592000.0"
+			target="/Outputs/vtkOutput"/-->
+		<PeriodicEvent
+			name="solver_1"
+			forceDt="86400.0"
+			endTime="604800.0"
+			target="/Solvers/compflow"/>
+		<PeriodicEvent
+			name="solver_2"
+			forceDt="1209600.0"
+			beginTime="604800.0"
+			target="/Solvers/compflow"/>
+	</Events>
+
+	<Outputs>
+		<VTK name="vtkOutput"/>
+	</Outputs>
+
+	<NumericalMethods>
+		<FiniteVolume>
+			<TwoPointFluxApproximation name="fluidTPFA"/>
+		</FiniteVolume>
+	</NumericalMethods>
+
+	<ElementRegions>
+		<CellElementRegion
+			name="Region1"
+			cellBlocks="{block1}"
+			materialList="{rock, fluid}"/>
+	</ElementRegions>
+
+	<FieldSpecifications>
+		<FieldSpecification
+			name="Region1InitialPressure"
+			initialCondition="1"
+			setNames="{ all }"
+			objectPath="ElementRegions/Region1"
+			fieldName="pressure"
+			scale="5000000.0"/>
+		<FieldSpecification
+			name="Region1InitialTemp"
+			initialCondition="1"
+			setNames="{ all }"
+			objectPath="ElementRegions/Region1"
+			fieldName="temperature"
+			scale="348.15"/>
+		<FieldSpecification
+			name="Region1InitCO2"
+			initialCondition="1"
+			setNames="{ all }"
+			objectPath="ElementRegions/Region1"
+			fieldName="globalCompFraction"
+			component="0"
+			scale="0.1"/>
+		<FieldSpecification
+			name="Region1InitC1"
+			initialCondition="1"
+			setNames="{ all }"
+			objectPath="ElementRegions/Region1"
+			fieldName="globalCompFraction"
+			component="1"
+			scale="0.2"/>
+		<FieldSpecification
+			name="Region1InitH2O"
+			initialCondition="1"
+			setNames="{ all }"
+			objectPath="ElementRegions/Region1"
+			fieldName="globalCompFraction"
+			component="2"
+			scale="0.7"/>
+
+		<FieldSpecification
+			name="inj1Pressure"
+			objectPath="ElementRegions/Region1"
+			fieldName="pressure"
+			scale="14000000.0"
+			setNames="{ inj1 }"/>
+		<FieldSpecification
+			name="inj1TermTemp"
+			setNames="{ inj1 }"
+			objectPath="ElementRegions/Region1"
+			fieldName="temperature"
+			scale="348.15"/>
+		<FieldSpecification
+			name="inj1CO2"
+			setNames="{ inj1 }"
+			objectPath="ElementRegions/Region1"
+			fieldName="globalCompFraction"
+			component="0"
+			scale="0.99999998"/>
+		<FieldSpecification
+			name="inj1C1"
+			setNames="{ inj1 }"
+			objectPath="ElementRegions/Region1"
+			fieldName="globalCompFraction"
+			component="1"
+			scale="1e-08"/>
+		<FieldSpecification
+			name="inj1H2O"
+			setNames="{ inj1 }"
+			objectPath="ElementRegions/Region1"
+			fieldName="globalCompFraction"
+			component="2"
+			scale="1e-08"/>
+
+		<FieldSpecification
+			name="prd1Pressure"
+			objectPath="ElementRegions/Region1"
+			fieldName="pressure"
+			scale="5000000.0"
+			setNames="{ prd1 }"/>
+		<FieldSpecification
+			name="prd1TermTemp"
+			setNames="{ prd1 }"
+			objectPath="ElementRegions/Region1"
+			fieldName="temperature"
+			scale="348.15"/>
+		<FieldSpecification
+			name="prd1CO2"
+			setNames="{ prd1 }"
+			objectPath="ElementRegions/Region1"
+			fieldName="globalCompFraction"
+			component="0"
+			scale="0.99999998"/>
+		<FieldSpecification
+			name="prd1C1"
+			setNames="{ prd1 }"
+			objectPath="ElementRegions/Region1"
+			fieldName="globalCompFraction"
+			component="1"
+			scale="1e-08"/>
+		<FieldSpecification
+			name="prd1H2O"
+			setNames="{ prd1 }"
+			objectPath="ElementRegions/Region1"
+			fieldName="globalCompFraction"
+			component="2"
+			scale="1e-08"/>
+	</FieldSpecifications>
+
+</Problem>

pygeos-tools/examples/obl/2ph_comp/main.py

@@ -0,0 +1,117 @@
+# ------------------------------------------------------------------------------------------------------------
+# SPDX-License-Identifier: LGPL-2.1-only
+#
+# Copyright (c) 2016-2024 Lawrence Livermore National Security LLC
+# Copyright (c) 2018-2024 TotalEnergies
+# Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University
+# Copyright (c) 2023-2024 Chevron
+# Copyright (c) 2019-     GEOS/GEOSX Contributors
+# Copyright (c) 2019-     INRIA project-team Makutu
+# All rights reserved
+#
+# See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
+# ------------------------------------------------------------------------------------------------------------
+
+# ---------------------------------------- README ----------------------------------------
+# Requires 'python -m pip install open-darts' and GEOS branch feature/anovikov/adaptive_obl
+# to run this example.
+# This is two-phase three-component model with fluid defined by constant K values.
+
+import numpy as np
+from mpi4py import MPI
+
+from geos.pygeos_tools.input import XML
+from geos.pygeos_tools.solvers import ReservoirSolver
+
+from darts.models.darts_model import DartsModel
+from darts.physics.super.physics import Compositional
+from darts.physics.super.property_container import PropertyContainer
+from darts.physics.properties.flash import ConstantK
+from darts.physics.properties.basic import ConstFunc, PhaseRelPerm
+from darts.physics.properties.density import DensityBasic
+
+
+class Model( DartsModel ):
+
+    def __init__( self, n_points=50 ):
+        # Call base class constructor
+        super().__init__()
+        self.n_obl_points = n_points
+        self.set_physics()
+
+    def set_physics( self ):
+        """Physical properties"""
+        self.zero = 1e-8
+        # Create property containers:
+        components = [ 'CO2', 'C1', 'H2O' ]
+        phases = [ 'gas', 'oil' ]
+        thermal = 0
+        Mw = [ 44.01, 16.04, 18.015 ]
+
+        property_container = PropertyContainer( phases_name=phases,
+                                                components_name=components,
+                                                Mw=Mw,
+                                                min_z=self.zero / 10,
+                                                temperature=1. )
+        """ properties correlations """
+        property_container.flash_ev = ConstantK( len( components ), [ 4, 2, 1e-1 ], self.zero )
+        property_container.density_ev = dict( [ ( 'gas', DensityBasic( compr=1e-3, dens0=200 ) ),
+                                                ( 'oil', DensityBasic( compr=1e-5, dens0=600 ) ) ] )
+        property_container.viscosity_ev = dict( [ ( 'gas', ConstFunc( 0.05 ) ), ( 'oil', ConstFunc( 0.5 ) ) ] )
+        property_container.rel_perm_ev = dict( [ ( 'gas', PhaseRelPerm( "gas" ) ), ( 'oil', PhaseRelPerm( "oil" ) ) ] )
+        """ Activate physics """
+        self.physics = Compositional( components,
+                                      phases,
+                                      self.timer,
+                                      n_points=self.n_obl_points,
+                                      min_p=1,
+                                      max_p=300,
+                                      min_z=self.zero / 10,
+                                      max_z=1 - self.zero / 10 )
+        self.physics.add_property_region( property_container )
+        self.engine = self.physics.init_physics( platform='cpu' )
+        return
+
+
+def run_darts_model( xml_name: str, darts_model=None ):
+    comm = MPI.COMM_WORLD
+    rank = comm.Get_rank()
+
+    xml = XML( xml_name )
+
+    solver = ReservoirSolver( solverType="ReactiveCompositionalMultiphaseOBL" )
+    solver.initialize( rank=rank, xml=xml )
+
+    # connect solver to Python-based operators
+    functions = solver.geosx.get_group( "/Functions" ).groups()
+    for func in functions:
+        if hasattr( func, 'setAxes' ) and darts_model is not None:
+            func.setAxes( darts_model.physics.n_vars, darts_model.physics.n_ops, list( darts_model.physics.axes_min ),
+                          list( darts_model.physics.axes_max ), list( darts_model.physics.n_axes_points ) )
+            func.setEvaluateFunction( darts_model.physics.reservoir_operators[ 0 ].evaluate )
+            print( "Adaptive OBL interpolator is configured." )
+
+    solver.applyInitialConditions()
+    solver.setDtFromTimeVariable( "forceDt" )
+    solver.setMaxTime( solver.getTimeVariables()[ "maxTime" ] )
+
+    time: float = 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.getDt():.4f}, iter = {cycle + 1}" )
+        solver.execute( time )
+        solver.outputVtk( time )
+        time += solver.dt
+        cycle += 1
+    solver.cleanup( time )
+
+
+if __name__ == "__main__":
+    # run adaptive OBL
+    print( "\n" + "=" * 30 + " RUNNING ADAPTIVE OBL " + "=" * 30 + "\n" )
+    darts_model = Model()
+    run_darts_model( xml_name="input_file_adaptive.xml", darts_model=darts_model )