feat: Adaptive multilinear interpolation for OBL solver

inputFiles/compositionalMultiphaseFlow/deadoil_3ph_staircase_obl_3d.xml

@@ -12,8 +12,7 @@
       enableEnergyBalance="0"
       useDARTSL2Norm="1"
       numComponents="3"
-      numPhases="3"
-      OBLOperatorsTableFile="obl_do_static.txt">
+      numPhases="3">
       <NonlinearSolverParameters
         newtonTol="1.0e-10"
         newtonMaxIter="10"/>
@@ -95,7 +94,7 @@
     <CellElementRegion
       name="Channel"
       cellBlocks="{ cb-1_0_0, cb-0_0_0, cb-0_0_1, cb-0_1_1, cb-0_1_2, cb-1_1_2, cb-1_1_3, cb-1_0_3 }"
-      materialList="{ rock }"/>
+      materialList="{ rock, fluid }"/>
 
     <CellElementRegion
       name="Barrier"
@@ -134,6 +133,11 @@
     <ConstantPermeability
       name="rockPerm"
       permeabilityComponents="{ 1.0e-16, 1.0e-16, 1.0e-16 }"/>
+
+    <OBLFluid
+			name="fluid"
+			interpolatorMode="static"
+			oblOperatorsTableFile="obl_do_static.txt"/>
   </Constitutive>
 
   <!-- END_SPHINX_INCLUDE_CONST_BLOCK -->

inputFiles/thermalMultiphaseFlow/co2_thermal_obl_3d.xml

@@ -10,8 +10,7 @@
       enableEnergyBalance="1"
       maxCompFractionChange="1"
       numComponents="3"
-      numPhases="2"
-      OBLOperatorsTableFile="ccs_thermal_10p_static.txt">
+      numPhases="2">
       <NonlinearSolverParameters
 	timeStepCutFactor="0.5"
 	newtonTol="1.0e-4"
@@ -79,7 +78,7 @@
     <CellElementRegion
       name="Region1"
       cellBlocks="{ * }"
-      materialList="{ rock }"/>
+      materialList="{ rock, fluid }"/>
   </ElementRegions>
 
   <Constitutive>
@@ -100,7 +99,10 @@
     <ConstantPermeability
       name="rockPerm"
       permeabilityComponents="{ 9.869233e-14, 9.869233e-14, 9.869233e-14 }"/>
-
+    <OBLFluid
+			name="fluid"
+			interpolatorMode="static"
+			oblOperatorsTableFile="ccs_thermal_10p_static.txt"/>
   </Constitutive>
 
   <FieldSpecifications>

src/coreComponents/constitutive/CMakeLists.txt

@@ -53,6 +53,8 @@ set( constitutive_headers
      dispersion/DispersionFields.hpp
      dispersion/DispersionSelector.hpp
      dispersion/LinearIsotropicDispersion.hpp
+     fluid/OBLFluid.hpp
+     fluid/OBLFluidKernels.hpp   
      fluid/multifluid/Layouts.hpp
      fluid/multifluid/LogLevelsInfo.hpp
      fluid/multifluid/MultiFluidSelector.hpp
@@ -230,6 +232,7 @@ set( constitutive_sources
      diffusion/DiffusionBase.cpp
      dispersion/DispersionBase.cpp
      dispersion/LinearIsotropicDispersion.cpp
+     fluid/OBLFluid.cpp
      fluid/multifluid/MultiFluidBase.cpp
      fluid/multifluid/blackOil/BlackOilFluidBase.cpp
      fluid/multifluid/blackOil/BlackOilFluid.cpp

src/coreComponents/constitutive/fluid/OBLFluid.cpp

@@ -0,0 +1,139 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * SPDX-License-Identifier: LGPL-2.1-only
+ *
+ * Copyright (c) 2018-2020 Lawrence Livermore National Security LLC
+ * Copyright (c) 2018-2020 The Board of Trustees of the Leland Stanford Junior University
+ * Copyright (c) 2018-2020 TotalEnergies
+ * Copyright (c) 2019-     GEOSX Contributors
+ * All rights reserved
+ *
+ * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
+ * ------------------------------------------------------------------------------------------------------------
+ */
+
+/**
+ * @file OBLFluid.cpp
+ */
+
+#include "OBLFluid.hpp"
+#include "functions/FunctionManager.hpp"
+#include "dataRepository/InputFlags.hpp"
+#include "dataRepository/RestartFlags.hpp"
+
+namespace geos
+{
+
+using namespace dataRepository;
+
+namespace constitutive
+{
+
+MultivariableTableFunction const *
+makeOBLOperatorsTable( string const & OBLOperatorsTableFile,
+                       string const & OBLFluidName,
+                       FunctionManager & functionManager )
+{
+  string const tableName = OBLFluidName + "OBL_table";
+  if( functionManager.hasGroup< MultivariableTableFunction >( tableName ) )
+  {
+    return functionManager.getGroupPointer< MultivariableTableFunction >( tableName );
+  }
+  else
+  {
+    MultivariableTableFunction * const table = dynamicCast< MultivariableTableFunction * >( functionManager.createChild( "MultivariableTableFunction", tableName ) );
+    table->initializeFunctionFromFile ( OBLOperatorsTableFile );
+    return table;
+  }
+}
+
+template< typename INDEX_T = __uint128_t >
+PythonFunction< INDEX_T > *
+makePythonFunction( string const & OBLFluidName, FunctionManager & functionManager )
+{
+  string const pythonFunctionName = OBLFluidName + "PythonFunction";
+  if( functionManager.hasGroup< PythonFunction< INDEX_T > >( pythonFunctionName ))
+  {
+    return functionManager.getGroupPointer< PythonFunction< INDEX_T > >( pythonFunctionName );
+  }
+  else
+  {
+    PythonFunction< INDEX_T > * function = dynamicCast< PythonFunction< INDEX_T > * >( functionManager.createChild( "PythonFunction", pythonFunctionName ));
+    return function;
+  }
+}
+
+OBLFluid::OBLFluid( string const & name, Group * const parent )
+  : ConstitutiveBase( name, parent ),
+  m_OBLOperatorsTable( nullptr ),
+  m_pythonFunction( nullptr )
+{
+  this->registerWrapper( viewKeyStruct::interpolatorModeString(), &m_interpolatorModeString ).
+    setInputFlag( InputFlags::REQUIRED ).
+    setDescription( "OBL interpolator mode: static or adaptive" );
+
+  this->registerWrapper( viewKeyStruct::interpolatorTypeString(), &m_interpolatorTypeString ).
+    setInputFlag( InputFlags::OPTIONAL ).
+    setDescription( "OBL interpolator type: multilinear or linear" );
+
+  this->registerWrapper( viewKeyStruct::oblOperatorsTableFileString(), &m_OBLOperatorsTableFile ).
+    setInputFlag( InputFlags::OPTIONAL ).
+    setDescription( "File containing OBL operator values for static mode interpolation" );
+}
+
+void OBLFluid::postInputInitialization()
+{
+  ConstitutiveBase::postInputInitialization();
+
+  // set interpolator mode
+  GEOS_THROW_IF( m_interpolatorModeString.empty(), 
+                 GEOS_FMT( "{}: Interpolator mode string is empty", 
+                            getFullName() ), 
+                 InputError );
+  m_interpolatorMode = EnumStrings<OBLInterpolatorMode>::fromString( m_interpolatorModeString );
+
+  // set interpolator type
+  /*GEOS_WARNING_IF( m_interpolatorTypeString.empty() ||
+                  (m_interpolatorTypeString != "multilinear" &&
+                   m_interpolatorTypeString != "linear"),
+                  GEOS_FMT( "{}: Invalid interpolator type, using multilinear interpolator",
+                            getFullName() ) );*/
+
+  GEOS_WARNING_IF( m_interpolatorTypeString == "linear",
+                   GEOS_FMT( "{}: Linear interpolator type is not supported yet, using multilinear interpolator",
+                             getFullName()) );
+
+  if( !m_interpolatorTypeString.empty())
+  {
+    if( m_interpolatorTypeString == "multilinear" )
+      m_interpolatorType = OBLInterpolatorType::Multilinear;
+    else if( m_interpolatorTypeString == "linear" )
+      m_interpolatorType = OBLInterpolatorType::Multilinear;
+    else
+      m_interpolatorType = OBLInterpolatorType::Multilinear;
+  }
+  else
+    m_interpolatorType = OBLInterpolatorType::Multilinear;
+
+  // set table file
+  GEOS_THROW_IF( m_OBLOperatorsTableFile.empty() &&
+                 m_interpolatorMode == OBLInterpolatorMode::Static,
+                 GEOS_FMT( "{}: Invalid operator table file",
+                           getFullName() ),
+                 InputError );
+
+
+  if( m_interpolatorMode == OBLInterpolatorMode::Static )
+    m_OBLOperatorsTable = makeOBLOperatorsTable( m_OBLOperatorsTableFile, getName(), FunctionManager::getInstance());
+  else if( m_interpolatorMode == OBLInterpolatorMode::Adaptive )
+    m_pythonFunction = makePythonFunction< longIndex >( getName(), FunctionManager::getInstance());
+
+  // raise intialization flag
+  m_isInitialized = true;
+}
+
+REGISTER_CATALOG_ENTRY( ConstitutiveBase, OBLFluid, string const &, Group * const )
+
+} // namespace constitutive
+
+} // namespace geos

src/coreComponents/constitutive/fluid/OBLFluid.hpp

@@ -0,0 +1,147 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * SPDX-License-Identifier: LGPL-2.1-only
+ *
+ * Copyright (c) 2018-2020 Lawrence Livermore National Security LLC
+ * Copyright (c) 2018-2020 The Board of Trustees of the Leland Stanford Junior University
+ * Copyright (c) 2018-2020 TotalEnergies
+ * Copyright (c) 2019-     GEOSX Contributors
+ * All rights reserved
+ *
+ * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
+ * ------------------------------------------------------------------------------------------------------------
+ */
+
+/**
+ * @file OBLFluid.hpp
+ */
+
+#ifndef GEOS_CONSTITUTIVE_FLUID_OBLFLUID_HPP_
+#define GEOS_CONSTITUTIVE_FLUID_OBLFLUID_HPP_
+
+#include "constitutive/ConstitutiveBase.hpp"
+#include "functions/MultivariableTableFunction.hpp"
+#include "functions/PythonFunction.hpp"
+
+namespace geos
+{
+namespace constitutive
+{
+
+enum OBLInterpolatorMode : integer
+{
+  Static = 0,         ///< static interpolation from a given tabled data
+  Adaptive = 1        ///< adaptive interpolation from a given interface to function
+};
+
+ENUM_STRINGS( OBLInterpolatorMode,
+              "static",
+              "adaptive" );
+
+enum OBLInterpolatorType : integer
+{
+  Multilinear = 0,    ///< multilinear interpolation
+  Linear = 1          ///< linear interpolation
+};
+
+class OBLFluid : public ConstitutiveBase
+{
+public:
+  using longIndex = __uint128_t;
+
+  OBLFluid( string const & name, Group * const parent );
+  /**
+   * @brief name of the constitutive in the object catalog
+   * @return string that contains the catalog name to generate a new object through the object catalog.
+   */
+  static string catalogName() { return "OBLFluid"; }
+  /**
+   * @copydoc ConstitutiveBase::getCatalogName()
+   */
+  virtual string getCatalogName() const override { return catalogName(); }
+  /**
+   * @struct viewKeyStruct holds char strings and viewKeys for fast lookup
+   */
+  struct viewKeyStruct : ConstitutiveBase::viewKeyStruct
+  {
+    // input
+    static constexpr char const * interpolatorModeString() { return "interpolatorMode"; }
+    static constexpr char const * interpolatorTypeString() { return "interpolatorType"; }
+    static constexpr char const * oblOperatorsTableFileString() { return "oblOperatorsTableFile"; }
+  };
+  /**
+   * @brief getter to the pointer to OBL operators table
+   * @return pointer to OBL operators table.
+   */
+  MultivariableTableFunction const & getTable() const
+  {
+    GEOS_ERROR_IF( m_OBLOperatorsTable == nullptr, "m_OBLOperatorsTable is not initialized" );
+    return *m_OBLOperatorsTable;
+  }
+  /**
+   * @brief getter to the Python-based evaluator
+   * @return pointer to the Python-based evaluator.
+   */
+  template< typename INDEX_T = longIndex >
+  PythonFunction< INDEX_T > * getPythonFunction()
+  {
+    GEOS_ERROR_IF( m_pythonFunction == nullptr, "m_pythonFunction is not initialized" );
+    return m_pythonFunction;
+  }
+  /**
+   * @brief initialize input
+   */
+  void initialize()
+  {
+    if( !m_isInitialized )
+    {
+      postInputInitialization();
+    }
+  }
+  /**
+   * @brief Retrieves the current OBL interpolator mode.
+   * @return OBLInterpolatorMode The current interpolation mode (Static or Adaptive).
+   */
+  OBLInterpolatorMode getInterpolatorMode() const { return m_interpolatorMode; };
+  /**
+   * @brief Retrieves the current OBL interpolator type.
+   * @return OBLInterpolatorType The current interpolation type (Multilinear or Linear).
+   */
+  OBLInterpolatorType getInterpolatorType() const { return m_interpolatorType; };
+private:
+  /// OBL interpolator mode
+  OBLInterpolatorMode m_interpolatorMode;
+
+  /// corresponding input string
+  string m_interpolatorModeString;
+
+  /// OBL interpolator type
+  OBLInterpolatorType m_interpolatorType;
+
+  /// corresponding input string
+  string m_interpolatorTypeString;
+
+  /// OBL operators table file (if OBL physics becomes consitutive, multiple regions will be supported )
+  Path m_OBLOperatorsTableFile;
+
+  /// OBL operators table function tabulated vs all primary variables
+  MultivariableTableFunction const * m_OBLOperatorsTable;
+
+  /// OBL operators with access to Python-base exact evaluator
+  PythonFunction< longIndex > * m_pythonFunction;
+
+  /// Flag to check if contitutive is initialized or not
+  bool m_isInitialized = false;
+
+  /**
+   * @copydoc dataRepository::Group::postInputInitialization()
+   */
+  virtual void postInputInitialization() override;
+};
+
+} //namespace constitutive
+
+} //namespace geos
+
+
+#endif //GEOS_CONSTITUTIVE_FLUID_OBLFLUID_HPP_