feat: Adaptive multilinear interpolation for OBL solver

inputFiles/compositionalMultiphaseFlow/2ph3c_square/2ph3c_square_obl.py

@@ -0,0 +1,121 @@
+# ------------------------------------------------------------------------------------------------------------
+# SPDX-License-Identifier: LGPL-2.1-only
+#
+# Copyright (c) 2016-2024 Lawrence Livermore National Security LLC
+# Copyright (c) 2018-2024 Total, S.A
+# 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.
+# ------------------------------------------------------------------------------------------------------------
+
+import numpy as np
+from mpi4py import MPI
+
+from geos.pygeos_tools.utilities.input import XML
+from geos.pygeos_tools.utilities.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()
+    solver.initialize(rank=rank, xml=xml)
+    solver.initialDt = solver.geosx.get_wrapper("Events/solver_1/forceDt").value()[0]
+    solver.dt = solver.initialDt
+
+    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.")
+
+    if rank == 0:
+       print()
+       print("===============================================================")
+       print(" Running pygeosx with solver:", solver.solver)
+       print(" simulation duration: ", solver.maxTime)
+       print(" initial timestep: ", solver.initialDt)
+       print("===============================================================")
+    solver.applyInitialConditions()
+
+    time = 0
+    cycle = 0
+
+    solver.outputVtk(time)
+    while time < solver.maxTime:
+        if rank == 0:
+          if solver.initialDt is not None:
+              print(f"time = {time:.3f}s, dt = {solver.initialDt:.4f}, iter = {cycle+1}")
+          else:
+              print(f"time = {time:.3f}s, dt is None, iter = {cycle+1}")
+        solver.execute(time)
+        solver.updateTimeStep()
+        solver.outputVtk(time)
+        time += solver.initialDt
+        cycle += 1
+    solver.cleanup(time)
+
+    comm.Barrier()
+
+
+# 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)
+
+# run static OBL
+print("\n" + "="*30 + " RUNNING STATIC OBL " + "="*30 + "\n")
+run_darts_model(xml_name="input_file_static.xml")

inputFiles/compositionalMultiphaseFlow/2ph3c_square/input_file_adaptive.xml

@@ -0,0 +1,216 @@
+<?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>
+
+	<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>
+
+	<NumericalMethods>
+		<FiniteVolume>
+			<TwoPointFluxApproximation name="fluidTPFA"/>
+		</FiniteVolume>
+	</NumericalMethods>
+
+	<ElementRegions>
+		<CellElementRegion
+			name="Region1"
+			cellBlocks="{block1}"
+			materialList="{rock, fluid}"/>
+	</ElementRegions>
+
+	<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>
+
+	<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>
+
+	<Outputs>
+		<VTK
+			name="vtkOutput"
+			plotFileRoot="2ph_square_adaptive"/>
+	</Outputs>
+</Problem>