diff --git a/inputFiles/compositionalMultiphaseFlow/deadoil_3ph_staircase_obl_3d.xml b/inputFiles/compositionalMultiphaseFlow/deadoil_3ph_staircase_obl_3d.xml index 0a4e50c5a93..8e8a1b64927 100644 --- a/inputFiles/compositionalMultiphaseFlow/deadoil_3ph_staircase_obl_3d.xml +++ b/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"> @@ -95,7 +94,7 @@ + materialList="{ rock, fluid }"/> + + diff --git a/inputFiles/thermalMultiphaseFlow/co2_thermal_obl_3d.xml b/inputFiles/thermalMultiphaseFlow/co2_thermal_obl_3d.xml index 1fe5b4a1b17..5740331aa6b 100644 --- a/inputFiles/thermalMultiphaseFlow/co2_thermal_obl_3d.xml +++ b/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"> + materialList="{ rock, fluid }"/> @@ -100,7 +99,10 @@ - + diff --git a/src/coreComponents/LvArray b/src/coreComponents/LvArray index e848125162b..3f3dce2502e 160000 --- a/src/coreComponents/LvArray +++ b/src/coreComponents/LvArray @@ -1 +1 @@ -Subproject commit e848125162b5b6af76d2dd40c8cfa7fd6ee18cbe +Subproject commit 3f3dce2502eb930c40f9e52648a972959238294e diff --git a/src/coreComponents/constitutive/CMakeLists.txt b/src/coreComponents/constitutive/CMakeLists.txt index b7ef68bb2c7..72348e8ec1b 100644 --- a/src/coreComponents/constitutive/CMakeLists.txt +++ b/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 @@ -242,6 +244,7 @@ set( constitutive_sources diffusion/DiffusionBase.cpp dispersion/DispersionBase.cpp dispersion/LinearIsotropicDispersion.cpp + fluid/OBLFluid.cpp fluid/multifluid/MultiFluidBase.cpp fluid/multifluid/constant/InvariantImmiscibleFluid.cpp fluid/multifluid/blackOil/BlackOilFluidBase.cpp diff --git a/src/coreComponents/constitutive/fluid/OBLFluid.cpp b/src/coreComponents/constitutive/fluid/OBLFluid.cpp new file mode 100644 index 00000000000..4b7dc754c75 --- /dev/null +++ b/src/coreComponents/constitutive/fluid/OBLFluid.cpp @@ -0,0 +1,131 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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; + } +} + +#if defined(GEOS_USE_PYGEOSX) +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; + } +} +#endif + +OBLFluid::OBLFluid( string const & name, Group * const parent ) + : ConstitutiveBase( name, parent ), + m_OBLOperatorsTable( nullptr ) +#if defined(GEOS_USE_PYGEOSX) + , m_pythonFunction( nullptr ) +#endif +{ + 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 ); + + // currently only multilinear interpolation is supported + 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 ) + { +#if defined(GEOS_USE_PYGEOSX) + m_pythonFunction = makePythonFunction< longIndex >( getName(), FunctionManager::getInstance()); +#else + GEOS_THROW( GEOS_FMT( "{}: adaptive interpolators require Python support, rebuild with ENABLE_PYGEOSX ON", getFullName() ), InputError ); +#endif + } + + // raise intialization flag + m_isInitialized = true; +} + +REGISTER_CATALOG_ENTRY( ConstitutiveBase, OBLFluid, string const &, Group * const ) + +} // namespace constitutive + +} // namespace geos diff --git a/src/coreComponents/constitutive/fluid/OBLFluid.hpp b/src/coreComponents/constitutive/fluid/OBLFluid.hpp new file mode 100644 index 00000000000..bcbd5ddbdba --- /dev/null +++ b/src/coreComponents/constitutive/fluid/OBLFluid.hpp @@ -0,0 +1,153 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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" +#if defined(GEOS_USE_PYGEOSX) + #include "functions/python/PythonFunction.hpp" +#endif + +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; + } +#if defined(GEOS_USE_PYGEOSX) + /** + * @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; + } +#endif + /** + * @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; + +#if defined(GEOS_USE_PYGEOSX) + /// OBL operators with access to Python-base exact evaluator + PythonFunction< longIndex > * m_pythonFunction; +#endif + + /// 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_ diff --git a/src/coreComponents/constitutive/fluid/OBLFluidKernels.hpp b/src/coreComponents/constitutive/fluid/OBLFluidKernels.hpp new file mode 100644 index 00000000000..62ee9cb4bac --- /dev/null +++ b/src/coreComponents/constitutive/fluid/OBLFluidKernels.hpp @@ -0,0 +1,321 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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 OBLFluidKernels.hpp + */ + +#ifndef GEOS_CONSTITUTIVE_FLUID_OBLFLUIDKERNELS_HPP_ +#define GEOS_CONSTITUTIVE_FLUID_OBLFLUIDKERNELS_HPP_ + +#include "mesh/ObjectManagerBase.hpp" +#include "constitutive/fluid/OBLFluid.hpp" +#include "functions/MultilinearInterpolatorStaticKernels.hpp" +#if defined(GEOS_USE_PYGEOSX) + #include "functions/python/MultilinearInterpolatorAdaptiveKernels.hpp" +#endif +#include "physicsSolvers/fluidFlow/CompositionalMultiphaseBaseFields.hpp" +#include "physicsSolvers/fluidFlow/ReactiveCompositionalMultiphaseOBLFields.hpp" + +namespace geos +{ + +namespace oblFluidKernels +{ + +using namespace constitutive; + +// The number of operators in use depends on: +// 1. Number of phases +// 2. Number of components +// 3. Features required in simulation (now its only one - Energy balance) +// This number needs to be used in solver and in kernels (as a template parameter) +// IMHO, this number is too big ( order of 10-100) to be treated via a kernelLaunchSelectorSwitch construct +// Hence, a way to define it once and for all is needed. +// Could be constexpr member of solver, but passing constexpr lambdas require -std=c++17 if I`m not mistaken + +constexpr integer COMPUTE_NUM_OPS ( integer const NP, integer const NC, bool ENERGY ) +{ + auto DOF = NC + ENERGY; + return DOF /*accumulation*/ + + DOF * NP /*flux*/ + + NP /*up_constant*/ + + DOF * NP /*gradient*/ + + DOF /*kinetic rate*/ + + 2 * NP /*gravity and capillarity*/ + + 1 /*rock porosity*/ + + NP /*enthalpies*/ + + 2 /*temperature and pressure*/; +} + + +/******************************** Kernel launch machinery ********************************/ +namespace internal +{ + +template< bool ENABLE_ENERGY, integer NUM_PHASES, typename T, typename LAMBDA > +void kernelLaunchSelectorCompSwitch( T numComps, LAMBDA && lambda ) +{ + static_assert( std::is_integral< T >::value, "kernelLaunchSelectorCompSwitch: type should be integral" ); + + switch( numComps ) + { + case 1: + { lambda( std::integral_constant< T, NUM_PHASES >(), std::integral_constant< T, 1 >(), std::integral_constant< bool, ENABLE_ENERGY >() ); return; } + case 2: + { lambda( std::integral_constant< T, NUM_PHASES >(), std::integral_constant< T, 2 >(), std::integral_constant< bool, ENABLE_ENERGY >() ); return; } + case 3: + { lambda( std::integral_constant< T, NUM_PHASES >(), std::integral_constant< T, 3 >(), std::integral_constant< bool, ENABLE_ENERGY >() ); return; } + case 4: + { lambda( std::integral_constant< T, NUM_PHASES >(), std::integral_constant< T, 4 >(), std::integral_constant< bool, ENABLE_ENERGY >() ); return; } + case 5: + { lambda( std::integral_constant< T, NUM_PHASES >(), std::integral_constant< T, 5 >(), std::integral_constant< bool, ENABLE_ENERGY >()); return; } + default: + { GEOS_ERROR( "Unsupported number of components: " << numComps ); } + } +} + +template< bool ENABLE_ENERGY, typename T, typename LAMBDA > +void kernelLaunchSelectorPhaseSwitch( T numPhases, T numComps, LAMBDA && lambda ) +{ + static_assert( std::is_integral< T >::value, "kernelLaunchSelectorPhaseSwitch: type should be integral" ); + + switch( numPhases ) + { + case 1: + { kernelLaunchSelectorCompSwitch< ENABLE_ENERGY, 1 >( numComps, lambda ); return; } + case 2: + { kernelLaunchSelectorCompSwitch< ENABLE_ENERGY, 2 >( numComps, lambda ); return; } + case 3: + { kernelLaunchSelectorCompSwitch< ENABLE_ENERGY, 3 >( numComps, lambda ); return; } + default: + { GEOS_ERROR( "Unsupported number of phases: " << numPhases ); } + } +} + +template< typename T, typename LAMBDA > +void kernelLaunchSelectorEnergySwitch( T numPhases, T numComps, bool enableEnergyBalance, LAMBDA && lambda ) +{ + if( enableEnergyBalance ) + { + kernelLaunchSelectorPhaseSwitch< true >( numPhases, numComps, lambda ); + } + else + { + kernelLaunchSelectorPhaseSwitch< false >( numPhases, numComps, lambda ); + } +} + +} // namespace internal + +/******************************** OBLOperatorsKernel ********************************/ + +/** + * @class OBLOperatorsKernel + * @tparam NUM_PHASES number of phases + * @tparam NUM_COMPS number of components + * @tparam ENABLE_ENERGY flag if energy balance equation is assembled + * @brief Compute OBL Operators and derivatives + */ +template< integer NUM_PHASES, integer NUM_COMPS, bool ENABLE_ENERGY > +class OBLOperatorsKernel +{ +public: + /// Compile time value for the energy balance switch + static constexpr integer enableEnergyBalance = ENABLE_ENERGY; + /// Compile time value for the number of components + static constexpr integer numComps = NUM_COMPS; + + /// Compile time value for the number of dimensions + static constexpr integer numDofs = numComps + enableEnergyBalance; + // /// Compile time value for the number of operators + static constexpr integer numOps = COMPUTE_NUM_OPS( NUM_PHASES, NUM_COMPS, ENABLE_ENERGY ); + + static constexpr real64 barToPascalMult = 1e5; + static constexpr real64 pascalToBarMult = 1e-5; + + /** + * @brief Performs the kernel launch + * @tparam POLICY the policy used in the RAJA kernels + * @tparam KERNEL_TYPE the kernel type + * @param[in] numElems the number of elements + * @param[inout] kernelComponent the kernel component providing access to the compute function + */ + template< typename POLICY, typename KERNEL_TYPE > + static void + launch( localIndex const numElems, + KERNEL_TYPE const & kernelComponent ) + { + forAll< POLICY >( numElems, [=] GEOS_HOST_DEVICE ( localIndex const ei ) + { + kernelComponent.compute( ei ); + } ); + } + + /** + * @brief Constructor + * @param[in] subRegion the element subregion + * @param[in] OBLOperatorsTable the OBL table function kernel + */ + OBLOperatorsKernel( ObjectManagerBase & subRegion, + MultilinearInterpolatorBaseKernel< numDofs, numOps > const * OBLOperatorsTable ) + : + m_OBLOperatorsTable( OBLOperatorsTable ), + m_pressure( subRegion.getField< fields::flow::pressure >() ), + m_compFrac( subRegion.getField< fields::flow::globalCompFraction >() ), + m_temperature( subRegion.getField< fields::flow::temperature >() ), + m_OBLOperatorValues ( subRegion.getField< fields::flow::OBLOperatorValues >()), + m_OBLOperatorDerivatives ( subRegion.getField< fields::flow::OBLOperatorDerivatives >()) + {} + + /** + * @brief Compute the operator values and derivatives for an element + * @param[in] ei the element index + */ + GEOS_HOST_DEVICE + inline + void compute( localIndex const ei ) const + { + arraySlice1d< real64 const, compflow::USD_COMP - 1 > const compFrac = m_compFrac[ei]; + arraySlice1d< real64, compflow::USD_OBL_VAL - 1 > const & OBLVals = m_OBLOperatorValues[ei]; + arraySlice2d< real64, compflow::USD_OBL_DER - 1 > const & OBLDers = m_OBLOperatorDerivatives[ei]; + real64 state[numDofs]; + + // we need to convert pressure from Pa (internal unit in GEOSX) to bar (internal unit in DARTS) + state[0] = m_pressure[ei] * pascalToBarMult; + + // the last component fraction is not used to define the state + for( integer i = 1; i < numComps; ++i ) + { + state[i] = compFrac[i - 1]; + } + + if( enableEnergyBalance ) + { + state[numDofs - 1] = m_temperature[ei]; + } + + m_OBLOperatorsTable->compute( state, OBLVals, OBLDers ); + + // we do not perform derivatives unit conversion here: + // instead we postpone it till all the derivatives are fully formed, and only then apply the factor only once in 'complete' function + // scaling the whole system might be even better solution (every pressure column needs to be multiplied by pascalToBarMult) + } + +private: + + // inputs + MultilinearInterpolatorBaseKernel< numDofs, numOps > const * m_OBLOperatorsTable; + + // Views on primary variables and their updates + arrayView1d< real64 const > m_pressure; + arrayView2d< real64 const, compflow::USD_COMP > m_compFrac; + arrayView1d< real64 const > m_temperature; + + // outputs + + // Views on OBL operator values and derivatives + arrayView2d< real64, compflow::USD_OBL_VAL > m_OBLOperatorValues; + arrayView3d< real64, compflow::USD_OBL_DER > m_OBLOperatorDerivatives; +}; + +/** + * @class OBLOperatorsKernelFactory + */ +class OBLOperatorsKernelFactory +{ +public: + + /** + * @brief Create a new kernel and launch + * @tparam POLICY the policy used in the RAJA kernel + * @param[in] numPhases the number of phases + * @param[in] numComponents the number of components + * @param[in] enableEnergyBalance flag if energy balance equation is assembled + * @param[in] subRegion the element subregion + * @param[in] function the OBL table function + */ + template< typename POLICY > + static void + createAndLaunch( integer const numPhases, + integer const numComponents, + bool const enableEnergyBalance, + ObjectManagerBase & subRegion, + constitutive::OBLFluid * oblFluid ) + { + internal::kernelLaunchSelectorEnergySwitch( numPhases, numComponents, enableEnergyBalance, [&] ( auto NP, auto NC, auto E ) + { + integer constexpr ENABLE_ENERGY = E(); + integer constexpr NUM_PHASES = NP(); + integer constexpr NUM_COMPS = NC(); + integer constexpr NUM_DIMS = ENABLE_ENERGY + NUM_COMPS; + integer constexpr NUM_OPS = COMPUTE_NUM_OPS( NUM_PHASES, NUM_COMPS, ENABLE_ENERGY ); + + if( oblFluid->getInterpolatorMode() == constitutive::OBLInterpolatorMode::Static ) + { + MultivariableTableFunction const & function = oblFluid->getTable(); + + GEOS_THROW_IF_NE_MSG( NUM_DIMS, function.numDims(), + GEOS_FMT( "The number of degrees of freedom per cell used in the solver has a value of {}, " + "whereas it as a value of {} in the operator table (at {}).", + NUM_DIMS, function.numDims(), + function.getName() ), + InputError ); + + GEOS_THROW_IF_NE_MSG( NUM_OPS, function.numOps(), + GEOS_FMT( "The number of operators per cell used in the solver has a value of {}, " + "whereas it as a value of {} in the operator table (at {}).", + NUM_OPS, function.numOps(), + function.getName() ), + InputError ); + + MultilinearInterpolatorStaticKernel< NUM_DIMS, NUM_OPS > const interpolationKernel( + function.getAxisMinimums(), + function.getAxisMaximums(), + function.getAxisPoints(), + function.getAxisSteps(), + function.getAxisStepInvs(), + function.getAxisHypercubeMults(), + function.getHypercubeData() + ); + + OBLOperatorsKernel< NUM_PHASES, NUM_COMPS, ENABLE_ENERGY > kernel( subRegion, &interpolationKernel ); + OBLOperatorsKernel< NUM_PHASES, NUM_COMPS, ENABLE_ENERGY >::template launch< POLICY >( subRegion.size(), kernel ); + } +#if defined(GEOS_USE_PYGEOSX) + else /* if ( oblFluid->getInterpolatorMode() == constitutive::OBLInterpolatorMode::Adaptive ) */ + { + // Check if Python function was assigned to wrapper + auto pyFunctionPtr = oblFluid->getPythonFunction(); + + GEOS_ERROR_IF( pyFunctionPtr->py_evaluate_func == nullptr, + GEOS_FMT( "{}: py_evaluate_func is not specified", + pyFunctionPtr->getName()) + ); + MultilinearInterpolatorAdaptiveKernel< NUM_DIMS, NUM_OPS > const interpolationKernel( pyFunctionPtr ); + + OBLOperatorsKernel< NUM_PHASES, NUM_COMPS, ENABLE_ENERGY > kernel( subRegion, &interpolationKernel ); + OBLOperatorsKernel< NUM_PHASES, NUM_COMPS, ENABLE_ENERGY >::template launch< POLICY >( subRegion.size(), kernel ); + } +#endif + } ); + } + +}; + +} // namespace oblFluidKernels + +} // namespace geos + +#endif //GEOS_CONSTITUTIVE_FLUID_OBLFLUIDKERNELS_HPP_ diff --git a/src/coreComponents/functions/CMakeLists.txt b/src/coreComponents/functions/CMakeLists.txt index b3869114d73..27c10135f8f 100644 --- a/src/coreComponents/functions/CMakeLists.txt +++ b/src/coreComponents/functions/CMakeLists.txt @@ -53,6 +53,15 @@ if( ENABLE_MATHPRESSO ) list( APPEND dependencyList mathpresso ) endif() +if( ENABLE_PYGEOSX ) + list( APPEND functions_headers + python/PythonFunction.hpp + python/PyPythonFunctionType.hpp ) + list( APPEND functions_sources + python/PyPythonFunction.cpp ) + list( APPEND dependencyList Python3::Python pylvarray ) +endif() + geos_decorate_link_dependencies( LIST decoratedDependencies DEPENDENCIES ${dependencyList} ) diff --git a/src/coreComponents/functions/FunctionManager.cpp b/src/coreComponents/functions/FunctionManager.cpp index 2f8425beb99..f68c0147d52 100644 --- a/src/coreComponents/functions/FunctionManager.cpp +++ b/src/coreComponents/functions/FunctionManager.cpp @@ -18,6 +18,9 @@ */ #include "FunctionManager.hpp" +#if defined(GEOS_USE_PYGEOSX) + #include "python/PythonFunction.hpp" +#endif namespace geos { @@ -54,11 +57,22 @@ Group * FunctionManager::createChild( string const & functionCatalogKey, string const & functionName ) { GEOS_LOG_RANK_0( GEOS_FMT( "{}: adding {} {}", getName(), functionCatalogKey, functionName ) ); - std::unique_ptr< FunctionBase > function = - FunctionBase::CatalogInterface::factory( functionCatalogKey, getDataContext(), functionName, this ); - return &this->registerGroup< FunctionBase >( functionName, std::move( function ) ); -} +#if defined(GEOS_USE_PYGEOSX) + if( functionCatalogKey == PythonFunction< __uint128_t >::catalogName() ) + { + // Create PythonFunction instance + std::unique_ptr< PythonFunction< __uint128_t > > function = std::make_unique< PythonFunction< __uint128_t > >( functionName, this ); + return &this->registerGroup< PythonFunction< __uint128_t > >( functionName, std::move( function )); + } + else +#endif + { + // Create FunctionBase-derived instance + std::unique_ptr< FunctionBase > function = FunctionBase::CatalogInterface::factory( functionCatalogKey, getDataContext(), functionName, this ); + return &this->registerGroup< FunctionBase >( functionName, std::move( function )); + } +} void FunctionManager::expandObjectCatalogs() { @@ -67,6 +81,11 @@ void FunctionManager::expandObjectCatalogs() { createChild( catalogIter.first, catalogIter.first ); } + +#if defined(GEOS_USE_PYGEOSX) + // Register an example of PythonFunction + createChild( "PythonFunction", "PythonFunction" ); +#endif } } // end of namespace geos diff --git a/src/coreComponents/functions/MultivariableTableFunctionKernels.hpp b/src/coreComponents/functions/MultilinearInterpolatorBaseKernels.hpp similarity index 72% rename from src/coreComponents/functions/MultivariableTableFunctionKernels.hpp rename to src/coreComponents/functions/MultilinearInterpolatorBaseKernels.hpp index 99022f23dc0..efdb5ec402b 100644 --- a/src/coreComponents/functions/MultivariableTableFunctionKernels.hpp +++ b/src/coreComponents/functions/MultilinearInterpolatorBaseKernels.hpp @@ -14,34 +14,27 @@ */ /** - * @file MultivariableTableFunctionKernels.hpp + * @file MultilinearInterpolatorBaseKernels.hpp */ -#ifndef GEOS_FUNCTIONS_MULTIVARIABLETABLEFUNCTIONKERNELS_HPP_ -#define GEOS_FUNCTIONS_MULTIVARIABLETABLEFUNCTIONKERNELS_HPP_ +#ifndef GEOS_FUNCTIONS_MULTILINEARINTERPOLATORBASEKERNELS_HPP_ +#define GEOS_FUNCTIONS_MULTILINEARINTERPOLATORBASEKERNELS_HPP_ namespace geos { - - /** - * @class KernelWrapper - * + * @class MultilinearInterpolatorBaseKernel * - */ - -/** - * @class MultivariableTableFunctionStaticKernel - * - * A class for multivariable piecewise interpolation with static storage + * A base class for multilinear piecewise interpolation * All functions are interpolated using the same uniformly discretized space * * @tparam NUM_DIMS number of dimensions (inputs) * @tparam NUM_OPS number of interpolated functions (outputs) + * @tparam INDEX_T datatype used for indexing in multidimensional space */ -template< integer NUM_DIMS, integer NUM_OPS > -class MultivariableTableFunctionStaticKernel +template< integer NUM_DIMS, integer NUM_OPS, typename INDEX_T = __uint128_t > +class MultilinearInterpolatorBaseKernel { public: @@ -54,116 +47,40 @@ class MultivariableTableFunctionStaticKernel /// Compile time value for the number of hypercube vertices static constexpr integer numVerts = 1 << numDims; + /// Datatype used for indexing + typedef INDEX_T longIndex; /** - * @brief Construct a new Multivariable Table Function Static Kernel object + * @brief Construct a new Multilinear Interpolator Base Kernel object * * @param[in] axisMinimums minimum coordinate for each axis * @param[in] axisMaximums maximum coordinate for each axis * @param[in] axisPoints number of discretization points between minimum and maximum for each axis * @param[in] axisSteps axis interval lengths (axes are discretized uniformly) * @param[in] axisStepInvs inversions of axis interval lengths (axes are discretized uniformly) - * @param[in] axisHypercubeMults hypercube index mult factors for each axis - * @param[in] hypercubeData table data stored per hypercube - */ - MultivariableTableFunctionStaticKernel( arrayView1d< real64 const > const & axisMinimums, - arrayView1d< real64 const > const & axisMaximums, - arrayView1d< integer const > const & axisPoints, - arrayView1d< real64 const > const & axisSteps, - arrayView1d< real64 const > const & axisStepInvs, - arrayView1d< globalIndex const > const & axisHypercubeMults, - arrayView1d< real64 const > const & hypercubeData ): - m_axisMinimums ( axisMinimums ), - m_axisMaximums ( axisMaximums ), - m_axisPoints ( axisPoints ), - m_axisSteps ( axisSteps ), - m_axisStepInvs ( axisStepInvs ), - m_axisHypercubeMults ( axisHypercubeMults ), - m_hypercubeData ( hypercubeData ) - {}; - -/** - * @brief interpolate all operators at a given point - * - * @param[in] coordinates point coordinates - * @param[out] values interpolated operator values - */ - template< typename IN_ARRAY, typename OUT_ARRAY > - GEOS_HOST_DEVICE - void - compute( IN_ARRAY const & coordinates, - OUT_ARRAY && values ) const - { - globalIndex hypercubeIndex = 0; - real64 axisLows[numDims]; - real64 axisMults[numDims]; - - for( int i = 0; i < numDims; ++i ) - { - integer const axisIndex = getAxisIntervalIndexLowMult( coordinates[i], - m_axisMinimums[i], m_axisMaximums[i], - m_axisSteps[i], m_axisStepInvs[i], m_axisPoints[i], - axisLows[i], axisMults[i] ); - hypercubeIndex += axisIndex * m_axisHypercubeMults[i]; - } - - interpolatePoint( coordinates, - getHypercubeData( hypercubeIndex ), - &axisLows[0], - &m_axisStepInvs[0], - values ); - } - - /** - * @brief interpolate all operators and compute their derivatives at a given point - * - * @param[in] coordinates point coordinates - * @param[out] values interpolated operator values - * @param[out] derivatives derivatives of interpolated operators + * @param[in] axisHypercubeMults hypercube index mult factors for each axis */ - template< typename IN_ARRAY, typename OUT_ARRAY, typename OUT_2D_ARRAY > - GEOS_HOST_DEVICE - void - compute( IN_ARRAY const & coordinates, - OUT_ARRAY && values, - OUT_2D_ARRAY && derivatives ) const + MultilinearInterpolatorBaseKernel( arrayView1d< real64 const > const & axisMinimums, + arrayView1d< real64 const > const & axisMaximums, + arrayView1d< integer const > const & axisPoints, + arrayView1d< real64 const > const & axisSteps, + arrayView1d< real64 const > const & axisStepInvs, + arrayView1d< longIndex const > const & axisHypercubeMults ): + m_axisMinimums( axisMinimums ), + m_axisMaximums( axisMaximums ), + m_axisPoints( axisPoints ), + m_axisSteps( axisSteps ), + m_axisStepInvs( axisStepInvs ), + m_axisHypercubeMults( axisHypercubeMults ), + m_axisPointMults( numDims ), + m_coordinates( numDims ) { - globalIndex hypercubeIndex = 0; - real64 axisLows[numDims]; - real64 axisMults[numDims]; - - for( int i = 0; i < numDims; ++i ) + // fill remaining properties + m_axisPointMults[numDims - 1] = 1; + for( integer i = numDims - 2; i >= 0; --i ) { - integer const axisIndex = getAxisIntervalIndexLowMult( coordinates[i], - m_axisMinimums[i], m_axisMaximums[i], - m_axisSteps[i], m_axisStepInvs[i], m_axisPoints[i], - axisLows[i], axisMults[i] ); - hypercubeIndex += axisIndex * m_axisHypercubeMults[i]; + m_axisPointMults[i] = m_axisPointMults[i + 1] * m_axisPoints[i + 1]; } - - interpolatePointWithDerivatives( coordinates, - getHypercubeData( hypercubeIndex ), - &axisLows[0], &axisMults[0], - &m_axisStepInvs[0], - values, - derivatives ); - - } - -protected: - - /** - * @brief Get pointer to hypercube data - * - * @param[in] hypercubeIndex - * @return pointer to hypercube data - */ - GEOS_HOST_DEVICE - inline - real64 const * - getHypercubeData( globalIndex const hypercubeIndex ) const - { - return &m_hypercubeData[hypercubeIndex * numVerts * numOps]; } /** @@ -221,7 +138,8 @@ class MultivariableTableFunctionStaticKernel /** * @brief interpolate all operators values at a given point * The algoritm is based on http://dx.doi.org/10.1090/S0025-5718-1988-0917826-0 - * + * @tparam IN_ARRAY type of input array of coordinates + * @tparam OUT_ARRAY type of output array of values * @param[in] axisCoordinates coordinates of a point * @param[in] hypercubeData data of target hypercube * @param[in] axisLows array of left coordinates of target axis intervals @@ -252,7 +170,6 @@ class MultivariableTableFunctionStaticKernel for( integer i = 0; i < numDims; ++i ) { - for( integer j = 0; j < pwr; ++j ) { for( integer op = 0; op < numOps; ++op ) @@ -273,7 +190,9 @@ class MultivariableTableFunctionStaticKernel /** * @brief interpolate all operators values and derivatives at a given point * The algoritm is based on http://dx.doi.org/10.1090/S0025-5718-1988-0917826-0 - * + * @tparam IN_ARRAY type of input array of coordinates + * @tparam OUT_ARRAY type of output array of values + * @tparam OUT_2D_ARRAY type of output array of derivatives * @param[in] axisCoordinates coordinates of a point * @param[in] hypercubeData data of target hypercube * @param[in] axisLows array of left coordinates of target axis intervals @@ -308,7 +227,6 @@ class MultivariableTableFunctionStaticKernel for( integer i = 0; i < numDims; ++i ) { - for( integer j = 0; j < pwr; ++j ) { for( integer op = 0; op < numOps; ++op ) @@ -345,48 +263,146 @@ class MultivariableTableFunctionStaticKernel } } } + /** + * @brief interpolate all operators at a given point + * + * @tparam IN_ARRAY type of input array of coordinates + * @tparam OUT_ARRAY type of output array of values + * @param[in] coordinates point coordinates + * @param[out] values interpolated operator values + */ + template< typename IN_ARRAY, typename OUT_ARRAY > + GEOS_HOST_DEVICE + void + compute( IN_ARRAY const & coordinates, + OUT_ARRAY && values ) const + { + longIndex hypercubeIndex = 0; + real64 axisLows[numDims]; + real64 axisMults[numDims]; + + for( int i = 0; i < numDims; ++i ) + { + integer const axisIndex = getAxisIntervalIndexLowMult( coordinates[i], + m_axisMinimums[i], m_axisMaximums[i], + m_axisSteps[i], m_axisStepInvs[i], m_axisPoints[i], + axisLows[i], axisMults[i] ); + hypercubeIndex += axisIndex * m_axisHypercubeMults[i]; + } + + interpolatePoint( coordinates, + getHypercubeData( hypercubeIndex ), + &axisLows[0], + &m_axisStepInvs[0], + values ); + } + + /** + * @brief interpolate all operators and compute their derivatives at a given point + * + * @tparam IN_ARRAY type of input array of coordinates + * @tparam OUT_ARRAY type of output array of values + * @tparam OUT_2D_ARRAY type of output array of derivatives + * @param[in] coordinates point coordinates + * @param[out] values interpolated operator values + * @param[out] derivatives derivatives of interpolated operators + */ + template< typename IN_ARRAY, typename OUT_ARRAY, typename OUT_2D_ARRAY > + GEOS_HOST_DEVICE + void + compute( IN_ARRAY const & coordinates, + OUT_ARRAY && values, + OUT_2D_ARRAY && derivatives ) const + { + longIndex hypercubeIndex = 0; + real64 axisLows[numDims]; + real64 axisMults[numDims]; + + for( int i = 0; i < numDims; ++i ) + { + integer const axisIndex = this->getAxisIntervalIndexLowMult( coordinates[i], + m_axisMinimums[i], m_axisMaximums[i], + m_axisSteps[i], m_axisStepInvs[i], m_axisPoints[i], + axisLows[i], axisMults[i] ); + hypercubeIndex += axisIndex * m_axisHypercubeMults[i]; + } + + interpolatePointWithDerivatives( coordinates, + getHypercubeData( hypercubeIndex ), + &axisLows[0], &axisMults[0], + &m_axisStepInvs[0], + values, + derivatives ); + } +protected: + /** + * @brief Get pointer to hypercube data + * + * @param[in] hypercubeIndex + * @return pointer to hypercube data + */ + virtual + GEOS_HOST_DEVICE + inline + real64 const * + getHypercubeData( longIndex const hypercubeIndex ) const + { + (void)hypercubeIndex; // Suppress unused parameter warning + return nullptr; + } + /** + * @brief Get coordinates of point given by index + * + * @param[in] pointIndex index of a point + * @param[out] coordinates coordinates of a point + */ + virtual + GEOS_HOST_DEVICE + inline + void + getPointCoordinates( longIndex pointIndex, array1d< real64 > & coordinates ) const + { + longIndex axisIdx, remainderIdx = pointIndex; + for( integer i = 0; i < numDims; ++i ) + { + axisIdx = remainderIdx / m_axisPointMults[i]; + remainderIdx = remainderIdx % m_axisPointMults[i]; + coordinates[i] = m_axisMinimums[i] + m_axisSteps[i] * axisIdx; + } + } + // inputs : table discretization data /// Array [numDims] of axis minimum values - arrayView1d< real64 const > m_axisMinimums; + arrayView1d< real64 const > const m_axisMinimums; /// Array [numDims] of axis maximum values - arrayView1d< real64 const > m_axisMaximums; + arrayView1d< real64 const > const m_axisMaximums; /// Array [numDims] of axis discretization points - arrayView1d< integer const > m_axisPoints; + arrayView1d< integer const > const m_axisPoints; // inputs : service data derived from table discretization data /// Array [numDims] of axis interval lengths (axes are discretized uniformly) - arrayView1d< real64 const > m_axisSteps; + arrayView1d< real64 const > const m_axisSteps; /// Array [numDims] of inversions of axis interval lengths (axes are discretized uniformly) - arrayView1d< real64 const > m_axisStepInvs; + arrayView1d< real64 const > const m_axisStepInvs; /// Array [numDims] of hypercube index mult factors for each axis - arrayView1d< globalIndex const > m_axisHypercubeMults; - - // inputs: operator sample data + arrayView1d< longIndex const > const m_axisHypercubeMults; - /// Main table data stored per hypercube: all values required for interpolation withing give hypercube are stored contiguously - arrayView1d< real64 const > m_hypercubeData; + /// Array [numDims] of point index mult factors for each axis + array1d< longIndex > const m_axisPointMults; // inputs: where to interpolate /// Coordinates in numDims-dimensional space where interpolation is requested - arrayView1d< real64 const > m_coordinates; - - // outputs - - /// Interpolated values - arrayView1d< real64 > m_values; - - /// /// Interpolated derivatives - arrayView1d< real64 > m_derivatives; + mutable array1d< real64 > m_coordinates; }; } /* namespace geos */ -#endif /* GEOS_FUNCTIONS_MULTIVARIABLETABLEFUNCTIONKERNELS_HPP_ */ +#endif /* GEOS_FUNCTIONS_MULTILINEARINTERPOLATORBASEKERNELS_HPP_ */ diff --git a/src/coreComponents/functions/MultilinearInterpolatorStaticKernels.hpp b/src/coreComponents/functions/MultilinearInterpolatorStaticKernels.hpp new file mode 100644 index 00000000000..4fc8adde901 --- /dev/null +++ b/src/coreComponents/functions/MultilinearInterpolatorStaticKernels.hpp @@ -0,0 +1,101 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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 + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file MultilinearInterpolatorStaticKernels.hpp + */ + +#ifndef GEOS_FUNCTIONS_MULTILINEARINTERPOLATORSTATICKERNELS_HPP_ +#define GEOS_FUNCTIONS_MULTILINEARINTERPOLATORSTATICKERNELS_HPP_ + +#include "MultilinearInterpolatorBaseKernels.hpp" + +namespace geos +{ +/** + * @class MultilinearInterpolatorStaticKernel + * + * A class for multilinear piecewise interpolation with static storage + * All functions are interpolated using the same uniformly discretized space + * + * @tparam NUM_DIMS number of dimensions (inputs) + * @tparam NUM_OPS number of interpolated functions (outputs) + * @tparam INDEX_T datatype used for indexing in multidimensional space + */ +template< integer NUM_DIMS, integer NUM_OPS, typename INDEX_T = __uint128_t > +class MultilinearInterpolatorStaticKernel : public MultilinearInterpolatorBaseKernel< NUM_DIMS, NUM_OPS, INDEX_T > +{ +public: + using MultilinearInterpolatorBaseKernel< NUM_DIMS, NUM_OPS, INDEX_T >::numDims; + using MultilinearInterpolatorBaseKernel< NUM_DIMS, NUM_OPS, INDEX_T >::numOps; + using MultilinearInterpolatorBaseKernel< NUM_DIMS, NUM_OPS, INDEX_T >::numVerts; + using MultilinearInterpolatorBaseKernel< NUM_DIMS, NUM_OPS, INDEX_T >::m_axisMinimums; + using MultilinearInterpolatorBaseKernel< NUM_DIMS, NUM_OPS, INDEX_T >::m_axisMaximums; + using MultilinearInterpolatorBaseKernel< NUM_DIMS, NUM_OPS, INDEX_T >::m_axisPoints; + using MultilinearInterpolatorBaseKernel< NUM_DIMS, NUM_OPS, INDEX_T >::m_axisSteps; + using MultilinearInterpolatorBaseKernel< NUM_DIMS, NUM_OPS, INDEX_T >::m_axisStepInvs; + using MultilinearInterpolatorBaseKernel< NUM_DIMS, NUM_OPS, INDEX_T >::m_axisHypercubeMults; + using typename MultilinearInterpolatorBaseKernel< NUM_DIMS, NUM_OPS, INDEX_T >::longIndex; + + /** + * @brief Construct a new Multilinear Interpolator Static Kernel object + * + * @param[in] axisMinimums minimum coordinate for each axis + * @param[in] axisMaximums maximum coordinate for each axis + * @param[in] axisPoints number of discretization points between minimum and maximum for each axis + * @param[in] axisSteps axis interval lengths (axes are discretized uniformly) + * @param[in] axisStepInvs inversions of axis interval lengths (axes are discretized uniformly) + * @param[in] axisHypercubeMults hypercube index mult factors for each axis + * @param[in] hypercubeData table data stored per hypercube + */ + MultilinearInterpolatorStaticKernel( arrayView1d< real64 const > const & axisMinimums, + arrayView1d< real64 const > const & axisMaximums, + arrayView1d< integer const > const & axisPoints, + arrayView1d< real64 const > const & axisSteps, + arrayView1d< real64 const > const & axisStepInvs, + arrayView1d< longIndex const > const & axisHypercubeMults, + arrayView1d< real64 const > const & hypercubeData ): + MultilinearInterpolatorBaseKernel< NUM_DIMS, NUM_OPS >( axisMinimums, + axisMaximums, + axisPoints, + axisSteps, + axisStepInvs, + axisHypercubeMults ), + m_hypercubeData ( hypercubeData ) + {}; +protected: + /** + * @brief Get pointer to hypercube data + * + * @param[in] hypercubeIndex + * @return pointer to hypercube data + */ + GEOS_HOST_DEVICE + inline + real64 const * + getHypercubeData( longIndex const hypercubeIndex ) const override + { + return &m_hypercubeData[hypercubeIndex * numVerts * numOps]; + } + + // inputs: operator sample data + + /// Main table data stored per hypercube: all values required for interpolation withing give hypercube are stored contiguously + arrayView1d< real64 const > const m_hypercubeData; +}; + +} /* namespace geos */ + +#endif /* GEOS_FUNCTIONS_MULTILINEARINTERPOLATORSTATICKERNELS_HPP_ */ diff --git a/src/coreComponents/functions/MultivariableTableFunction.hpp b/src/coreComponents/functions/MultivariableTableFunction.hpp index 927a020d817..1eb61698239 100644 --- a/src/coreComponents/functions/MultivariableTableFunction.hpp +++ b/src/coreComponents/functions/MultivariableTableFunction.hpp @@ -32,7 +32,7 @@ namespace geos * @class MultivariableTableFunction * * An interface class for multivariable table function (function with multiple inputs and outputs) with uniform discretization - * Prepares input data for MultivariableStaticInterpolatorKernel, which performes actual interpolation + * Prepares input data for MultilinearInterpolatorStaticKernel, which performes actual interpolation */ class MultivariableTableFunction : public FunctionBase @@ -157,7 +157,7 @@ class MultivariableTableFunction : public FunctionBase * @brief Get the table axes hypercube index multiplicators * @return a reference to an array of table axes hypercube index multiplicators */ - arrayView1d< globalIndex const > getAxisHypercubeMults() const { return m_axisHypercubeMults.toViewConst(); } + arrayView1d< __uint128_t const > getAxisHypercubeMults() const { return m_axisHypercubeMults.toViewConst(); } /** * @brief Get the table values stored per-hypercube @@ -214,10 +214,10 @@ class MultivariableTableFunction : public FunctionBase real64_array m_axisStepInvs; /// Array [numDims] of point index mult factors for each axis - globalIndex_array m_axisPointMults; + array1d< __uint128_t > m_axisPointMults; /// Array [numDims] of hypercube index mult factors for each axis - globalIndex_array m_axisHypercubeMults; + array1d< __uint128_t > m_axisHypercubeMults; // inputs: operator sample data diff --git a/src/coreComponents/functions/python/MultilinearInterpolatorAdaptiveKernels.hpp b/src/coreComponents/functions/python/MultilinearInterpolatorAdaptiveKernels.hpp new file mode 100644 index 00000000000..1f559965f14 --- /dev/null +++ b/src/coreComponents/functions/python/MultilinearInterpolatorAdaptiveKernels.hpp @@ -0,0 +1,173 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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 + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file MultilinearInterpolatorAdaptiveKernels.hpp + */ + +#ifndef GEOS_FUNCTIONS_MULTILINEARINTERPOLATORADAPTIVEKERNELS_HPP_ +#define GEOS_FUNCTIONS_MULTILINEARINTERPOLATORADAPTIVEKERNELS_HPP_ + +#include "functions/MultilinearInterpolatorBaseKernels.hpp" +#include "functions/python/PythonFunction.hpp" + +namespace geos +{ +/** + * @class KernelWrapper + * + * + */ + +/** + * @class MultilinearInterpolatorAdaptiveKernel + * + * A class for multilinear piecewise interpolation with access to exact function evaluator + * All functions are interpolated using the same uniformly discretized space + * + * @tparam NUM_DIMS number of dimensions (inputs) + * @tparam NUM_OPS number of interpolated functions (outputs) + */ +template< integer NUM_DIMS, integer NUM_OPS, typename INDEX_T = __uint128_t > +class MultilinearInterpolatorAdaptiveKernel : public MultilinearInterpolatorBaseKernel< NUM_DIMS, NUM_OPS, INDEX_T > +{ +public: + using MultilinearInterpolatorBaseKernel< NUM_DIMS, NUM_OPS, INDEX_T >::numDims; + using MultilinearInterpolatorBaseKernel< NUM_DIMS, NUM_OPS, INDEX_T >::numOps; + using MultilinearInterpolatorBaseKernel< NUM_DIMS, NUM_OPS, INDEX_T >::numVerts; + using MultilinearInterpolatorBaseKernel< NUM_DIMS, NUM_OPS, INDEX_T >::m_axisMinimums; + using MultilinearInterpolatorBaseKernel< NUM_DIMS, NUM_OPS, INDEX_T >::m_axisMaximums; + using MultilinearInterpolatorBaseKernel< NUM_DIMS, NUM_OPS, INDEX_T >::m_axisPoints; + using MultilinearInterpolatorBaseKernel< NUM_DIMS, NUM_OPS, INDEX_T >::m_axisSteps; + using MultilinearInterpolatorBaseKernel< NUM_DIMS, NUM_OPS, INDEX_T >::m_axisStepInvs; + using MultilinearInterpolatorBaseKernel< NUM_DIMS, NUM_OPS, INDEX_T >::m_axisHypercubeMults; + using MultilinearInterpolatorBaseKernel< NUM_DIMS, NUM_OPS, INDEX_T >::m_axisPointMults; + using MultilinearInterpolatorBaseKernel< NUM_DIMS, NUM_OPS, INDEX_T >::m_coordinates; + using typename MultilinearInterpolatorBaseKernel< NUM_DIMS, NUM_OPS, INDEX_T >::longIndex; + + /** + * @brief Construct a new Multilinear Interpolator Adaptive Kernel object + * @param[in] evalFunc evaluator of exact function values + */ + MultilinearInterpolatorAdaptiveKernel( PythonFunction< longIndex > * evalFunc ): + MultilinearInterpolatorBaseKernel< NUM_DIMS, NUM_OPS >( evalFunc->getAxisMinimums(), + evalFunc->getAxisMaximums(), + evalFunc->getAxisPoints(), + evalFunc->getAxisSteps(), + evalFunc->getAxisStepInvs(), + evalFunc->getAxisHypercubeMults() ), + m_evalFunc( evalFunc ) + {}; +protected: + /** + * @brief Get values of operators at a given point + * Provide a reference to correct location in the adaptive point storage. + * If the point is not found, compute it first, and then return the reference. + * + * @param[in] pointIndex index of point + * @return operator values at given point + */ + GEOS_HOST_DEVICE + array1d< real64 > const & + getPointData( longIndex const pointIndex ) const + { + auto & m_pointData = m_evalFunc->getPointData(); + auto item = m_pointData.find( pointIndex ); + if( item == m_pointData.end()) + { + array1d< real64 > newPoint ( numOps ); + this->getPointCoordinates( pointIndex, m_coordinates ); + m_evalFunc->evaluate( m_coordinates, newPoint ); + m_pointData[pointIndex] = newPoint; + return m_pointData.at( pointIndex ); + } + else + return item->second; + } + /** + * @brief Get indexes of all vertices for given hypercube + * + * @param[in] hypercubeIndex index of the hyporcube + * @param[out] hypercubePoints indexes of all vertices of hypercube + */ + GEOS_HOST_DEVICE + inline + void + getHypercubePoints( longIndex const hypercubeIndex, + array1d< longIndex > & hypercubePoints ) const + { + longIndex axisIdx, remainderIdx = hypercubeIndex; + integer pwr = numVerts; + + for( integer i = 0; i < numVerts; ++i ) + hypercubePoints[i] = 0; + + for( integer i = 0; i < numDims; ++i ) + { + axisIdx = remainderIdx / m_axisHypercubeMults[i]; + remainderIdx = remainderIdx % m_axisHypercubeMults[i]; + pwr /= 2; + + for( integer j = 0; j < numVerts; ++j ) + { + integer zeroOrOne = (j / pwr) % 2; + hypercubePoints[j] += (axisIdx + zeroOrOne) * m_axisPointMults[i]; + } + } + } + /** + * @brief Get pointer to hypercube data + * + * @param[in] hypercubeIndex + * @return pointer to hypercube data + */ + GEOS_HOST_DEVICE + inline + real64 const * + getHypercubeData( longIndex const hypercubeIndex ) const override + { + auto & m_hypercubeData = m_evalFunc->getHypercubeData(); + auto item = m_hypercubeData.find( hypercubeIndex ); + if( item == m_hypercubeData.end()) + { + array1d< longIndex > points ( numVerts ); + array1d< real64 > newHypercube ( numVerts * numOps ); + + this->getHypercubePoints( hypercubeIndex, points ); + + for( integer i = 0; i < numVerts; ++i ) + { + // obtain point data and copy it to hypercube data + auto const & pData = this->getPointData( points[i] ); + for( integer op = 0; op < numOps; op++ ) + { + newHypercube[i * numOps + op] = pData[op]; + } + } + m_hypercubeData[hypercubeIndex] = newHypercube; + return m_hypercubeData[hypercubeIndex].data(); + } + else + return item->second.data(); + } + /** + * @brief wrapper providing interface to Python function + */ + PythonFunction< longIndex > * m_evalFunc; +}; + +} /* namespace geos */ + +#endif /* GEOS_FUNCTIONS_MULTILINEARINTERPOLATORADAPTIVEKERNELS_HPP_ */ diff --git a/src/coreComponents/functions/python/PyPythonFunction.cpp b/src/coreComponents/functions/python/PyPythonFunction.cpp new file mode 100644 index 00000000000..ffcbc1a56bb --- /dev/null +++ b/src/coreComponents/functions/python/PyPythonFunction.cpp @@ -0,0 +1,206 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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 + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +// Python.h must be included first. +#define PY_SSIZE_T_CLEAN +#include + +// Source includes +#include "LvArray/src/python/pythonHelpers.hpp" +#include "PythonFunction.hpp" +#include "PyPythonFunctionType.hpp" +#include "dataRepository/python/PyGroupType.hpp" + + +#define VERIFY_NON_NULL_SELF( self ) \ + PYTHON_ERROR_IF( self == nullptr, PyExc_RuntimeError, "Passed a nullptr as self.", nullptr ) + +#define VERIFY_INITIALIZED( self ) \ + PYTHON_ERROR_IF( self->group == nullptr, PyExc_RuntimeError, "The PyPythonFunction is not initialized.", nullptr ) + +namespace geos +{ +namespace python +{ + +struct PyPythonFunction +{ + PyObject_HEAD + + static constexpr char const * docString = + "A Python interface to geos::PythonFunction."; + + geos::PythonFunction<> * group; +}; + + +/** + * @brief Allocate a new PyPythonFunction object + */ +static PyObject * PyPythonFunction_new( PyTypeObject *type, PyObject *args, PyObject *kwds ) +{ + GEOS_UNUSED_VAR( args, kwds ); + PyPythonFunction *self; + + // Allocate memory for the Python object + self = (PyPythonFunction *)type->tp_alloc( type, 0 ); + if( self != nullptr ) + { + self->group = nullptr; // Initialize the pyFunc to nullptr (not yet assigned) + } + + return (PyObject *)self; +} +/** + * @brief String representation of PyPythonFunction object + */ +static PyObject * PyPythonFunction_repr( PyObject * const obj ) noexcept +{ + PyPythonFunction const * const pyPythonFunc = reinterpret_cast< PyPythonFunction * >( obj ); + if( pyPythonFunc == nullptr ) + { + return nullptr; + } + + VERIFY_INITIALIZED( pyPythonFunc ); + + std::string const path = pyPythonFunc->group->getPath(); // Assuming PythonFunction has getPath() + std::string const type = LvArray::system::demangle( typeid( *(pyPythonFunc->group) ).name() ); + std::string const repr = path + " ( " + type + " )"; + + return PyUnicode_FromString( repr.c_str() ); +} +/** + * @brief Expose the setAxes function to Python + */ +static PyObject * PyPythonFunction_setAxes( PyPythonFunction * self, PyObject * args ) +{ + VERIFY_NON_NULL_SELF( self ); + VERIFY_INITIALIZED( self ); + + PythonFunction<> * func = self->group; + + if( func == nullptr ) + { + PyErr_SetString( PyExc_RuntimeError, "Not a valid PythonFunction instance." ); + return nullptr; + } + + // Parse arguments + integer numDims, numOps; + PyObject *axisMinList, *axisMaxList, *axisPointsList; + + if( !PyArg_ParseTuple( args, "iiOOO", &numDims, &numOps, &axisMinList, &axisMaxList, &axisPointsList )) + { + return nullptr; + } + + // Check list lengths + if( PyList_Size( axisMinList ) != numDims || PyList_Size( axisMaxList ) != numDims || PyList_Size( axisPointsList ) != numDims ) + { + PyErr_SetString( PyExc_ValueError, "List lengths do not match numDims." ); + return nullptr; + } + + // Convert Python lists to internal array representations + real64_array axisMinimums( numDims ); + real64_array axisMaximums( numDims ); + integer_array axisPoints( numDims ); + + for( integer i = 0; i < numDims; ++i ) + { + axisMinimums[i] = PyFloat_AsDouble( PyList_GetItem( axisMinList, i )); + axisMaximums[i] = PyFloat_AsDouble( PyList_GetItem( axisMaxList, i )); + axisPoints[i] = PyLong_AsLong( PyList_GetItem( axisPointsList, i )); + } + + // Call the C++ function to set axes + func->setAxes( numDims, numOps, axisMinimums, axisMaximums, axisPoints ); + + Py_RETURN_NONE; +} + +/** + * @brief Expose the setEvaluateFunction function to Python + */ +static PyObject * PyPythonFunction_setEvaluateFunction( PyPythonFunction * self, PyObject * args ) +{ + VERIFY_NON_NULL_SELF( self ); + VERIFY_INITIALIZED( self ); + + PythonFunction<> * func = self->group; + + if( func == nullptr ) + { + PyErr_SetString( PyExc_RuntimeError, "Not a valid PythonFunction instance." ); + return nullptr; + } + + PyObject * pyFuncObj; + + if( !PyArg_ParseTuple( args, "O", &pyFuncObj )) + { + return nullptr; + } + + // Ensure pyFuncObj is callable + if( !PyCallable_Check( pyFuncObj )) + { + PyErr_SetString( PyExc_TypeError, "Argument must be callable." ); + return nullptr; + } + + // Call the C++ function to set the evaluate function + func->setEvaluateFunction( pyFuncObj ); + + Py_RETURN_NONE; +} + +static PyMethodDef PyPythonFunction_methods[] = { + { "setAxes", (PyCFunction) PyPythonFunction_setAxes, METH_VARARGS, "Set axes for the Python function." }, + { "setEvaluateFunction", (PyCFunction) PyPythonFunction_setEvaluateFunction, METH_VARARGS, "Set Python evaluate function and optionally a derivative function." }, + { nullptr, nullptr, 0, nullptr } // Sentinel +}; + +/** + * Initialize the module object for Python with the exported functions + */ + +BEGIN_ALLOW_DESIGNATED_INITIALIZERS + +static PyTypeObject PyPythonFunctionType = { + PyVarObject_HEAD_INIT( nullptr, 0 ) + .tp_name = "pygeosx.PythonFunction", + .tp_basicsize = sizeof(PyPythonFunction), + .tp_repr = PyPythonFunction_repr, + .tp_flags = Py_TPFLAGS_DEFAULT | Py_TPFLAGS_BASETYPE, + .tp_doc = PyPythonFunction::docString, + .tp_methods = PyPythonFunction_methods, + .tp_base = getPyGroupType(), + .tp_new = PyPythonFunction_new, +}; + +END_ALLOW_DESIGNATED_INITIALIZERS + +/** + * Return the PyTypeObject for PyPythonFunction + */ +PyTypeObject * getPyPythonFunctionType() +{ + return &PyPythonFunctionType; +} + +} // namespace python +} // namespace geos diff --git a/src/coreComponents/functions/python/PyPythonFunctionType.hpp b/src/coreComponents/functions/python/PyPythonFunctionType.hpp new file mode 100644 index 00000000000..95b37f277f1 --- /dev/null +++ b/src/coreComponents/functions/python/PyPythonFunctionType.hpp @@ -0,0 +1,31 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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 + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +#ifndef GEOS_PYTHON_PYPYTHONFUNCTIONTYPE_HPP_ +#define GEOS_PYTHON_PYPYTHONFUNCTIONTYPE_HPP_ + +#include "LvArray/src/python/pythonForwardDeclarations.hpp" + +namespace geos +{ +namespace python +{ + +PyTypeObject * getPyPythonFunctionType(); + +} +} + +#endif diff --git a/src/coreComponents/functions/python/PythonFunction.hpp b/src/coreComponents/functions/python/PythonFunction.hpp new file mode 100644 index 00000000000..d32a1c52cee --- /dev/null +++ b/src/coreComponents/functions/python/PythonFunction.hpp @@ -0,0 +1,330 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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 + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file PythonFunction.hpp + */ + +#ifndef GEOS_FUNCTIONS_PYTHONFUNCTION_HPP_ +#define GEOS_FUNCTIONS_PYTHONFUNCTION_HPP_ + +#include +#include "common/DataTypes.hpp" +#include "dataRepository/Group.hpp" +#include "LvArray/src/python/PyFunc.hpp" +#include "PyPythonFunctionType.hpp" + +namespace std +{ +template<> +struct hash< __uint128_t > +{ + size_t operator()( const __uint128_t & x ) const noexcept + { + size_t h1 = std::hash< uint64_t >{} (static_cast< uint64_t >(x)); + size_t h2 = std::hash< uint64_t >{} (static_cast< uint64_t >(x >> 64)); + return h1 ^ (h2 * 0x9e3779b97f4a7c15 + 0x7f4a7c15); // Use a large prime multiplier and a random offset + } +}; +}; + +namespace geos +{ +/** + * @class PythonFunction + * @brief A wrapper class to interface Python functions for use in C++ computations. + * + * @tparam IN_ARRAY type of input array of coordinates + * @tparam OUT_ARRAY type of output array of values + * @tparam OUT_2D_ARRAY type of output array of derivatives + * @tparam INDEX_T datatype used for indexing in multidimensional space + */ +template< typename INDEX_T = __uint128_t > +class PythonFunction : public dataRepository::Group +{ +public: + /// Number of table dimensions (inputs) + integer numDims; + + /// Number of operators (interpolated functions, outputs) + integer numOps; + + /// Number of hypercube vertices + integer numVerts; + + /// Datatype used for indexing + typedef INDEX_T longIndex; + + typedef LvArray::python::PythonFunction< array1d< real64 > const &, array1d< real64 > & > PyEvaluateFunction; + typedef LvArray::python::PythonFunction< array1d< real64 > const &, array2d< real64 > & > PyDerivativeFunction; + + /// Python function reference to be called in evaluate + // PyObject* py_evaluate_func = nullptr; + mutable std::unique_ptr< PyEvaluateFunction > py_evaluate_func; + + /// Python function reference to be called for evaluation of derivatives + // PyObject* py_derivative_func = nullptr; + mutable std::unique_ptr< PyDerivativeFunction > py_derivative_func; + + /** + * @brief Construct a new wrapper for python function + * @param[in] name the name of this object manager + * @param[in] parent the parent Group + */ + PythonFunction( const string & name, Group * const parent ): + dataRepository::Group( name, parent ), + py_evaluate_func( nullptr ), + py_derivative_func( nullptr ) + { + if( !Py_IsInitialized()) + Py_Initialize(); + } + /** + * @brief The catalog name interface + * @return name of the PythonFunction in the FunctionBase catalog + */ + static string catalogName() { return "PythonFunction"; } + /** + * @brief Set Python functions for values and (optionally) for their derivatives + * + * @param[in] pyFunc pointer to the function + * @param[in] pyDerivativeFunc pointer to the function evaluating derivatives + */ + void setEvaluateFunction( PyObject * pyFunc, PyObject * pyDerivativeFunc = nullptr ) + { + // Set the evaluation function (required) + if( PyCallable_Check( pyFunc )) + py_evaluate_func = std::make_unique< PyEvaluateFunction >( pyFunc ); + else + GEOS_THROW( GEOS_FMT( "{}: Provided Python function is not callable.", getName()), InputError ); + + // Set the derivative function (optional) + if( pyDerivativeFunc ) + { + if( PyCallable_Check( pyDerivativeFunc )) + py_derivative_func = std::make_unique< PyDerivativeFunction >( pyDerivativeFunc ); + else + GEOS_THROW( GEOS_FMT( "{}: Provided Python derivative function is not callable.", getName()), InputError ); + } + } + /** + * @brief Set dimensions of input and output mapping spaces + * + * @param[in] numberOfDimesions dimension of input states + * @param[in] numberOfOperators dimension of output operators + */ + void setDimensions( integer numberOfDimesions, + integer numberOfOperators ) + { + numDims = numberOfDimesions; + numOps = numberOfOperators; + numVerts = 1 << numDims; + } + /** + * @brief Set the properties of parametrization + * + * @param[in] numberOfDimesions dimension of input states + * @param[in] numberOfOperators dimension of output operators + * @param[in] axisMinimums array with state axes' minimums + * @param[in] axisMaximums array with state axes' maximums + * @param[in] axisPoints array defining number of points along each state axis + */ + void setAxes( integer numberOfDimesions, + integer numberOfOperators, + real64_array axisMinimums, + real64_array axisMaximums, + integer_array axisPoints ) + { + numDims = numberOfDimesions; + numOps = numberOfOperators; + numVerts = 1 << numDims; + + m_axisMinimums = axisMinimums; + m_axisMaximums = axisMaximums; + m_axisPoints = axisPoints; + + m_axisSteps.resize( numDims ); + m_axisStepInvs.resize( numDims ); + m_axisHypercubeMults.resize( numDims ); + + for( integer dim = 0; dim < numDims; dim++ ) + { + m_axisSteps[dim] = (m_axisMaximums[dim] - m_axisMinimums[dim]) / (m_axisPoints[dim] - 1); + m_axisStepInvs[dim] = 1 / m_axisSteps[dim]; + } + + m_axisHypercubeMults[numDims - 1] = 1; + for( integer dim = numDims - 2; dim >= 0; --dim ) + { + m_axisHypercubeMults[dim] = m_axisHypercubeMults[dim + 1] * (m_axisPoints[dim + 1] - 1); + } + } + /** + * @brief interpolate all operators values at a given point + * + * @param[in] state function arguments + * @param[out] values function values + */ + GEOS_HOST_DEVICE + inline + void + evaluate( array1d< real64 > const & state, + array1d< real64 > & values ) const + { + if( py_evaluate_func ) + (*py_evaluate_func)( state, values ); + else + GEOS_THROW( GEOS_FMT( "{}: Evaluation function not set.", getName()), InputError ); + } + + /** + * @brief interpolate all operators values at a given point + * + * @param[in] state function arguments + * @param[out] values function values + * @param[out] derivatives derivatives of function values w.r.t. arguments + */ + GEOS_HOST_DEVICE + inline + void + evaluate( array1d< real64 > const & state, + array1d< real64 > & values, + array2d< real64 > & derivatives ) const + { + if( py_evaluate_func ) + (*py_evaluate_func)( state, values ); + else + GEOS_THROW( GEOS_FMT( "{}: Evaluation function not set.", getName()), InputError ); + + if( py_derivative_func ) + (*py_derivative_func)( state, derivatives ); + else + GEOS_THROW( GEOS_FMT( "{}: Derivative function not set.", getName()), InputError ); + } + /** + * @brief Retrieves constant point data. + * @return const unordered_map>& + * A constant reference to the point data mapping. + */ + GEOS_HOST_DEVICE + inline + unordered_map< longIndex, array1d< real64 > > const & + getPointData() const { return m_pointData; } + /** + * @brief Retrieves modifiable point data. + * @return unordered_map>& + * A reference to the point data mapping, which can be modified by the caller. + */ + GEOS_HOST_DEVICE + inline + unordered_map< longIndex, array1d< real64 > > & + getPointData() { return m_pointData; } + /** + * @brief Retrieves constant hypercube data. + * @return const unordered_map>& + * A constant reference to the hypercube data mapping. + */ + GEOS_HOST_DEVICE + inline + unordered_map< longIndex, array1d< real64 > > const & + getHypercubeData() const { return m_hypercubeData; } + /** + * @brief Retrieves modifiable hypercube data. + * @return unordered_map>& + * A reference to the hypercube data mapping, which can be modified by the caller. + */ + GEOS_HOST_DEVICE + inline + unordered_map< longIndex, array1d< real64 > > & + getHypercubeData() { return m_hypercubeData; } + /** + * @brief Get the axes minimums + * @return a reference to an array of axes minimums + */ + arrayView1d< real64 const > getAxisMinimums() const { return m_axisMinimums.toViewConst(); } + /** + * @brief Get the axes maximums + * @return a reference to an array of axes maximums + */ + arrayView1d< real64 const > getAxisMaximums() const { return m_axisMaximums.toViewConst(); } + /** + * @brief Get the axes discretization points numbers + * @return a reference to an array of axes discretization points numbers + */ + arrayView1d< integer const > getAxisPoints() const { return m_axisPoints.toViewConst(); } + /** + * @brief Get the axes step sizes + * @return a reference to an array of axes step sizes + */ + arrayView1d< real64 const > getAxisSteps() const { return m_axisSteps.toViewConst(); } + /** + * @brief Get the axes step sizes inversions + * @return a reference to an array of axes step sizes inversions + */ + arrayView1d< real64 const > getAxisStepInvs() const { return m_axisStepInvs.toViewConst(); } + /** + * @brief Get the table axes hypercube index multiplicators + * @return a reference to an array of table axes hypercube index multiplicators + */ + arrayView1d< __uint128_t const > getAxisHypercubeMults() const { return m_axisHypercubeMults.toViewConst(); } + /** + * @brief Return PySolver type. + * @return Return PySolver type. + */ + virtual PyTypeObject * getPythonType() const override + { + return python::getPyPythonFunctionType(); + } +private: + /// Array [numDims] of axis minimum values + array1d< real64 > m_axisMinimums; + + /// Array [numDims] of axis maximum values + array1d< real64 > m_axisMaximums; + + /// Array [numDims] of axis discretization points numbers + array1d< integer > m_axisPoints; + + /// Array [numDims] of axis interval lengths (axes are discretized uniformly) + array1d< real64 > m_axisSteps; + + /// Array [numDims] of inversions of axis interval lengths (axes are discretized uniformly) + array1d< real64 > m_axisStepInvs; + + /// Array [numDims] of hypercube index mult factors for each axis + array1d< __uint128_t > m_axisHypercubeMults; + + /** + * @brief adaptive point storage: the values of operators at requested supporting points + * Storage is grown dynamically in the process of simulation. + * Only supporting points that are required for interpolation are computed and added + */ + unordered_map< longIndex, array1d< real64 > > m_pointData; + /** + * @brief adaptive hypercube storage: the values of operators at every vertex of reqested hypercubes + * Storage is grown dynamically in the process of simulation + * Only hypercubes that are required for interpolation are computed and added + * + * In fact it is an excess storage used to reduce memory accesses during interpolation. + * Here all values of all vertexes of requested hypercube are stored consecutevely and are accessed via a single index + * Usage of point_data for interpolation directly would require N_VERTS memory accesses (>1000 accesses for 10-dimensional space) + * * + */ + unordered_map< longIndex, array1d< real64 > > m_hypercubeData; +}; + +} /* namespace geos */ + +#endif /* GEOS_FUNCTIONS_PYTHONFUNCTION_HPP_ */ diff --git a/src/coreComponents/functions/unitTests/testFunctions.cpp b/src/coreComponents/functions/unitTests/testFunctions.cpp index 04ca26acca6..fbce32c39fc 100644 --- a/src/coreComponents/functions/unitTests/testFunctions.cpp +++ b/src/coreComponents/functions/unitTests/testFunctions.cpp @@ -20,9 +20,16 @@ #include "functions/FunctionBase.hpp" #include "functions/TableFunction.hpp" #include "functions/MultivariableTableFunction.hpp" -#include "functions/MultivariableTableFunctionKernels.hpp" +#include "functions/MultilinearInterpolatorStaticKernels.hpp" +#include "common/logger/Logger.hpp" //#include "mainInterface/GeosxState.hpp" +#if defined(GEOS_USE_PYGEOSX) + #include + #include "functions/python/PythonFunction.hpp" + #include "functions/python/MultilinearInterpolatorAdaptiveKernels.hpp" +#endif + #ifdef GEOS_USE_MATHPRESSO #include "functions/SymbolicFunction.hpp" #endif @@ -643,14 +650,15 @@ void testMutivariableFunction( MultivariableTableFunction & function, arrayView2d< real64 > evaluatedDerivativesView = evaluatedDerivatives.toView(); - MultivariableTableFunctionStaticKernel< NUM_DIMS, NUM_OPS > kernel( function.getAxisMinimums(), - function.getAxisMaximums(), - function.getAxisPoints(), - function.getAxisSteps(), - function.getAxisStepInvs(), - function.getAxisHypercubeMults(), - function.getHypercubeData() - ); + MultilinearInterpolatorStaticKernel< NUM_DIMS, NUM_OPS > kernel( function.getAxisMinimums(), + function.getAxisMaximums(), + function.getAxisPoints(), + function.getAxisSteps(), + function.getAxisStepInvs(), + function.getAxisHypercubeMults(), + function.getHypercubeData() + ); + // Test values evaluation first forAll< geos::parallelDevicePolicy< > >( numElems, [=] GEOS_HOST_DEVICE ( localIndex const elemIndex ) @@ -886,6 +894,216 @@ TEST( FunctionTests, MultivariableTableFromFile ) testMutivariableFunction< nDims, nOps >( table_h, testCoordinates, testExpectedValues, testExpectedDerivatives ); } +#if defined(GEOS_USE_PYGEOSX) + +std::pair< PyObject *, PyObject * > definePythonFunction() +{ + const char * pythonCode = + R"( +import math +import pylvarray + +def python_evaluate(state, values): + x1, x2, x3, x4, x5, x6 = state.to_numpy() + values = values.to_numpy() + + values[0] = (math.sin(x1) * math.cos(x2) * math.sin(x3) * + math.cos(x4) * math.sin(x5) * math.cos(x6)) + + values[1] = (math.sin(x1) + math.cos(x2) + math.sin(x3) + + math.cos(x4) + math.sin(x5) + math.cos(x6)) + +def python_derivatives(state, derivatives): + x1, x2, x3, x4, x5, x6 = state.to_numpy() + derivatives = derivatives.to_numpy() + + derivatives[0][0] = math.cos(x1) * math.cos(x2) * math.sin(x3) * math.cos(x4) * math.sin(x5) * math.cos(x6) + derivatives[0][1] = -math.sin(x1) * math.sin(x2) * math.sin(x3) * math.cos(x4) * math.sin(x5) * math.cos(x6) + derivatives[0][2] = math.sin(x1) * math.cos(x2) * math.cos(x3) * math.cos(x4) * math.sin(x5) * math.cos(x6) + derivatives[0][3] = -math.sin(x1) * math.cos(x2) * math.sin(x3) * math.sin(x4) * math.sin(x5) * math.cos(x6) + derivatives[0][4] = math.sin(x1) * math.cos(x2) * math.sin(x3) * math.cos(x4) * math.cos(x5) * math.cos(x6) + derivatives[0][5] = -math.sin(x1) * math.cos(x2) * math.sin(x3) * math.cos(x4) * math.sin(x5) * math.sin(x6) + + derivatives[1][0] = math.cos(x1) + derivatives[1][1] = -math.sin(x2) + derivatives[1][2] = math.cos(x3) + derivatives[1][3] = -math.sin(x4) + derivatives[1][4] = math.cos(x5) + derivatives[1][5] = -math.sin(x6) +)"; + + PyRun_SimpleString( pythonCode ); + + PyObject * mainModule = PyImport_AddModule( "__main__" ); + PyObject * func = PyObject_GetAttrString( mainModule, "python_evaluate" ); + PyObject * deriv_func = PyObject_GetAttrString( mainModule, "python_derivatives" ); + + if( func == nullptr || + deriv_func == nullptr || + !PyCallable_Check( func ) || + !PyCallable_Check( deriv_func )) + { + PyErr_Print(); + GEOS_THROW( "Failed to define Python function or derivatives.", InputError ); + } + + return {func, deriv_func}; +} + +template< integer NUM_DIMS, integer NUM_OPS > +void getMultilinearAdaptiveInterpolation( integer numAxisPts, + real64 lowerBound, + real64 upperBound, + std::vector< array1d< real64 > > const & coordinates, + arrayView2d< real64, compflow::USD_OBL_VAL > const & values, + arrayView3d< real64, compflow::USD_OBL_DER > const & derivatives ) +{ + // Set parameters of state space grid + array1d< real64 > const axesMinimums ( NUM_DIMS ); + array1d< real64 > const axesMaximums ( NUM_DIMS ); + array1d< integer > const axesPoints ( NUM_DIMS ); + for( integer dim = 0; dim < NUM_DIMS; dim++ ) + { + axesMinimums[dim] = lowerBound; + axesMaximums[dim] = upperBound; + axesPoints[dim] = numAxisPts; + } + + // Create interface to Python functions + FunctionManager * functionManager = &FunctionManager::getInstance(); + string const pythonFunctionName = "Evaluator_" + std::to_string( numAxisPts ); + PythonFunction<> & func = dynamicCast< PythonFunction<> & >( *functionManager->createChild( "PythonFunction", pythonFunctionName ) ); + auto [pyFunc, pyDerivFunc] = definePythonFunction(); + func.setEvaluateFunction( pyFunc, pyDerivFunc ); + func.setAxes( NUM_DIMS, NUM_OPS, axesMinimums, axesMaximums, axesPoints ); + + // Create interpolation kernel + MultilinearInterpolatorAdaptiveKernel< NUM_DIMS, NUM_OPS > kernel( &func ); + + // Evaluation + const integer numPts = coordinates.size(); + forAll< geos::parallelDevicePolicy< > >( numPts, + [=] GEOS_HOST_DEVICE ( localIndex const ptIndex ) + { + arrayView1d< real64 const > const coordinatesView = coordinates[ptIndex].toViewConst(); + arraySlice1d< real64, compflow::USD_OBL_VAL - 1 > const & ptValues = values[ptIndex]; + arraySlice2d< real64, compflow::USD_OBL_DER - 1 > const & ptDers = derivatives[ptIndex]; + + kernel.compute( coordinatesView, ptValues, ptDers ); + } ); +} + +TEST( FunctionTests, MultilinearInterpolatorAdaptiveKernels ) +{ + // Test the convergence rate of multilinear adaptive interpolation + // Calculates exact values and derivatives defined by Python functions in random points + // Then interpolates using two resolutions and estimate convergence rate + + constexpr integer numDims = 6; + constexpr integer numOps = 2; + constexpr real64 lowerBound = -M_PI; + constexpr real64 upperBound = M_PI; + constexpr integer numPts = 100; + constexpr integer numResolutions = 2; + constexpr std::array< integer, numResolutions > numAxesPts = {128, 512}; + array1d< real64 > residuals ( 2 * numResolutions ); + + // Fill array of coordinates + std::vector< array1d< real64 > > states ( numPts, array1d< real64 >( numDims ) ); + std::default_random_engine generator; + std::uniform_real_distribution< double > distribution( lowerBound, upperBound ); + for( integer i = 0; i < numPts; ++i ) + { + arrayView1d< real64 > const stateView = states[i].toView(); + for( integer j = 0; j < numDims; ++j ) + stateView[j] = distribution( generator ); + } + + // Initialize output arrays + using array2dLayoutOBLOpVals = array2d< real64, compflow::LAYOUT_OBL_OPERATOR_VALUES >; + using array3dLayoutOBLOpDers = array3d< real64, compflow::LAYOUT_OBL_OPERATOR_DERIVATIVES >; + array2dLayoutOBLOpVals evaluatedValues( numPts, numOps ); + array3dLayoutOBLOpDers evaluatedDerivatives ( numPts, numOps, numDims ); + + // Other types for exact values due to restriction on type of arguments of Python function + std::vector< array1d< real64 > > trueValues( numPts, array1d< real64 >( numOps )); + std::vector< array2d< real64 > > trueDerivatives( numPts, array2d< real64 >( numOps, numDims )); + + // Create interface to Python functions for evaluation of exact values + FunctionManager * functionManager = &FunctionManager::getInstance(); + PythonFunction<> & func = dynamicCast< PythonFunction<> & >( *functionManager->createChild( "PythonFunction", "ExactEvaluator" ) ); + auto [pyFunc, pyDerivFunc] = definePythonFunction(); + func.setEvaluateFunction( pyFunc, pyDerivFunc ); + func.setDimensions( numDims, numOps ); + + // Python section + Py_Initialize(); + + // Do exact evaluations + forAll< geos::parallelDevicePolicy< > >( numPts, + [&] GEOS_HOST_DEVICE ( localIndex const ptIndex ) + { + func.evaluate( states[ptIndex], trueValues[ptIndex], trueDerivatives[ptIndex] ); + } ); + + // Do multilinear interpolation + arrayView2d< real64, compflow::USD_OBL_VAL > const m_evaluatedValues ( evaluatedValues ); + arrayView3d< real64, compflow::USD_OBL_DER > const m_evaluatedDerivatives ( evaluatedDerivatives ); + for( integer i = 0; i < numResolutions; ++i ) + { + getMultilinearAdaptiveInterpolation< numDims, numOps >( numAxesPts[i], + lowerBound, + upperBound, + states, + m_evaluatedValues, + m_evaluatedDerivatives ); + + residuals[2 * i] = residuals[2 * i + 1] = 0.0; + for( integer ptIndex = 0; ptIndex < numPts; ++ptIndex ) + { + arrayView1d< real64 > const ptTrueValues = trueValues[ptIndex].toView(); + arrayView2d< real64 > const ptTrueDers = trueDerivatives[ptIndex].toView(); + + arraySlice1d< real64, compflow::USD_OBL_VAL - 1 > const & ptEvalValues = m_evaluatedValues[ptIndex]; + arraySlice2d< real64, compflow::USD_OBL_DER - 1 > const & ptEvalDers = m_evaluatedDerivatives[ptIndex]; + + for( integer opIndex = 0; opIndex < numOps; ++opIndex ) + { + real64 diff = ptEvalValues[opIndex] - ptTrueValues[opIndex]; + residuals[2 * i] += diff * diff; + + for( integer axIndex = 0; axIndex < numDims; ++axIndex ) + { + real64 d_diff = ptEvalDers[opIndex][axIndex] - ptTrueDers[opIndex][axIndex]; + residuals[2 * i + 1] += d_diff * d_diff; + } + } + } + residuals[2 * i] /= numPts * numOps; + residuals[2 * i + 1] /= numPts * numOps * numDims; + + residuals[2 * i] = sqrt( residuals[2 * i] ); + residuals[2 * i + 1] = sqrt( residuals[2 * i + 1] ); + + // printf("Residuals with %d points per axes: %f \t %f \n", numAxesPts[i], residuals[2 * i], residuals[2 * i + 1]); + } + + Py_Finalize(); + + // Calculate convergence rate + const real64 dLogDx = log((upperBound - lowerBound) / numAxesPts[numResolutions - 1] ) - + log((upperBound - lowerBound) / numAxesPts[numResolutions - 2] ); + const real64 convOrderVals = (log( residuals[2 * (numResolutions - 1)] ) - + log( residuals[2 * (numResolutions - 2)] )) / dLogDx; + const real64 convOrderDers = (log( residuals[2 * (numResolutions - 1) + 1] ) - + log( residuals[2 * (numResolutions - 2) + 1] )) / dLogDx; + // printf("Covergence order: values %f, \t derivatives %f \n", convOrderVals, convOrderDers); + + ASSERT_GT( convOrderVals, 1.8 ) << "interpolation convergence rate must be greater than 1.8"; + ASSERT_GT( convOrderDers, 0.8 ) << "derivatives convergence rate must be greater than 0.8"; +} + +#endif // The `ENUM_STRING` implementation relies on consistency between the order of the `enum`, // and the order of the `string` array provided. Since this consistency is not enforced, it can be corrupted anytime. diff --git a/src/coreComponents/linearAlgebra/interfaces/hypre/mgrStrategies/ReactiveCompositionalMultiphaseOBL.hpp b/src/coreComponents/linearAlgebra/interfaces/hypre/mgrStrategies/ReactiveCompositionalMultiphaseOBL.hpp index ee2561ad491..77ace7566f9 100644 --- a/src/coreComponents/linearAlgebra/interfaces/hypre/mgrStrategies/ReactiveCompositionalMultiphaseOBL.hpp +++ b/src/coreComponents/linearAlgebra/interfaces/hypre/mgrStrategies/ReactiveCompositionalMultiphaseOBL.hpp @@ -68,8 +68,8 @@ class ReactiveCompositionalMultiphaseOBL : public MGRStrategyBase< 1 > m_levelFRelaxType[0] = MGRFRelaxationType::none; m_levelInterpType[0] = MGRInterpolationType::injection; - m_levelRestrictType[0] = MGRRestrictionType::injection; - m_levelCoarseGridMethod[0] = MGRCoarseGridMethod::cprLikeBlockDiag; + m_levelRestrictType[0] = MGRRestrictionType::blockColLumped; + m_levelCoarseGridMethod[0] = MGRCoarseGridMethod::galerkin; m_levelGlobalSmootherType[0] = MGRGlobalSmootherType::ilu0; m_levelGlobalSmootherIters[0] = 1; } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/ReactiveCompositionalMultiphaseOBL.cpp b/src/coreComponents/physicsSolvers/fluidFlow/ReactiveCompositionalMultiphaseOBL.cpp index af7b3b8e51b..715c2da387a 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/ReactiveCompositionalMultiphaseOBL.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/ReactiveCompositionalMultiphaseOBL.cpp @@ -20,6 +20,7 @@ #include "ReactiveCompositionalMultiphaseOBL.hpp" #include "constitutive/solid/CoupledSolidBase.hpp" +#include "constitutive/fluid/OBLFluid.hpp" #include "dataRepository/Group.hpp" #include "discretizationMethods/NumericalMethodsManager.hpp" #include "fieldSpecification/FieldSpecificationManager.hpp" @@ -44,26 +45,6 @@ using namespace constitutive; using namespace fields; using namespace reactiveCompositionalMultiphaseOBLKernels; -namespace -{ - -MultivariableTableFunction const * makeOBLOperatorsTable( string const & OBLOperatorsTableFile, - FunctionManager & functionManager ) -{ - string const tableName = "OBL_operators_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; - } -} -} - ReactiveCompositionalMultiphaseOBL::ReactiveCompositionalMultiphaseOBL( const string & name, Group * const parent ) : @@ -78,6 +59,11 @@ ReactiveCompositionalMultiphaseOBL::ReactiveCompositionalMultiphaseOBL( const st setInputFlag( InputFlags::REQUIRED ). setDescription( "Number of components" ); + this->registerWrapper( viewKeyStruct::numSolidComponentsString(), &m_numSolidComponents ). + setApplyDefaultValue( 0 ). + setInputFlag( InputFlags::OPTIONAL ). + setDescription( "Number of solid components" ); + this->registerWrapper( viewKeyStruct::numPhasesString(), &m_numPhases ). setInputFlag( InputFlags::REQUIRED ). setDescription( "Number of phases" ); @@ -86,11 +72,6 @@ ReactiveCompositionalMultiphaseOBL::ReactiveCompositionalMultiphaseOBL( const st setInputFlag( InputFlags::REQUIRED ). setDescription( "Enable energy balance calculation and temperature degree of freedom" ); - this->registerWrapper( viewKeyStruct::OBLOperatorsTableFileString(), &m_OBLOperatorsTableFile ). - setInputFlag( InputFlags::REQUIRED ). - setRestartFlags( RestartFlags::NO_WRITE ). - setDescription( "File containing OBL operator values" ); - this->registerWrapper( viewKeyStruct::maxCompFracChangeString(), &m_maxCompFracChange ). setApplyDefaultValue( 1.0 ). setInputFlag( InputFlags::OPTIONAL ). @@ -187,6 +168,15 @@ void ReactiveCompositionalMultiphaseOBL::postInputInitialization() // need to override to skip the check for fluidModel, which is enabled in FlowSolverBase PhysicsSolverBase::postInputInitialization(); + GEOS_THROW_IF_GT_MSG( m_numSolidComponents, m_numComponents, + GEOS_FMT( "{}: The number of solid component is set to {}, while it must not be greater than the total number of components {}", + getWrapperDataContext( viewKeyStruct::numSolidComponentsString() ), m_numSolidComponents, m_numComponents ), + InputError ); + GEOS_THROW_IF_LT_MSG( m_numSolidComponents, 0, + GEOS_FMT( "{}: The number of solid component is set to {}, while it must non-negative", + getWrapperDataContext( viewKeyStruct::numSolidComponentsString() ), m_numSolidComponents ), + InputError ); + GEOS_THROW_IF_GT_MSG( m_maxCompFracChange, 1.0, GEOS_FMT( "{}: The maximum absolute change in component fraction is set to {}, while it must not be greater than 1.0", getWrapperDataContext( viewKeyStruct::maxCompFracChangeString() ), m_maxCompFracChange ), @@ -196,30 +186,11 @@ void ReactiveCompositionalMultiphaseOBL::postInputInitialization() getWrapperDataContext( viewKeyStruct::maxCompFracChangeString() ), m_maxCompFracChange ), InputError ); - m_OBLOperatorsTable = makeOBLOperatorsTable( m_OBLOperatorsTableFile, FunctionManager::getInstance()); - // Equations: [NC] Molar mass balance, ([1] energy balance if enabled) // Primary variables: [1] pressure, [NC-1] global component fractions, ([1] temperature) m_numDofPerCell = m_numComponents + m_enableEnergyBalance; m_numOBLOperators = COMPUTE_NUM_OPS( m_numPhases, m_numComponents, m_enableEnergyBalance ); - - GEOS_THROW_IF_NE_MSG( m_numDofPerCell, m_OBLOperatorsTable->numDims(), - GEOS_FMT( "The number of degrees of freedom per cell used in the solver (at {}) has a value of {}, " - "whereas it as a value of {} in the operator table (at {}).", - getWrapperDataContext( viewKeyStruct::elemDofFieldString() ), - m_numDofPerCell, m_OBLOperatorsTable->numDims(), - m_OBLOperatorsTableFile ), - InputError ); - - GEOS_THROW_IF_NE_MSG( m_numOBLOperators, m_OBLOperatorsTable->numOps(), - GEOS_FMT( "The number of operators per cell used in the solver (at {}) has a value of {}, " - "whereas it as a value of {} in the operator table (at {}).", - getWrapperDataContext( viewKeyStruct::elemDofFieldString() ), - m_numDofPerCell, m_OBLOperatorsTable->numDims(), - m_OBLOperatorsTableFile ), - InputError ); - } void ReactiveCompositionalMultiphaseOBL::registerDataOnMesh( Group & meshBodies ) @@ -1228,6 +1199,7 @@ void ReactiveCompositionalMultiphaseOBL::chopPrimaryVariables( DomainPartition & { real64 const eps = LvArray::NumericLimits< real64 >::epsilon; integer const numComp = m_numComponents; + integer const numSolidComp = m_numSolidComponents; forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, MeshLevel & mesh, @@ -1248,7 +1220,7 @@ void ReactiveCompositionalMultiphaseOBL::chopPrimaryVariables( DomainPartition & // Step 1: chop the component fractions between 0 and 1 real64 sum = 0.0; - for( integer ic = 0; ic < numComp-1; ++ic ) + for( integer ic = numSolidComp; ic < numComp-1; ++ic ) { if( compFrac[ei][ic] > 1.0-eps ) { @@ -1275,7 +1247,7 @@ void ReactiveCompositionalMultiphaseOBL::chopPrimaryVariables( DomainPartition & // Step 3: rescale the component fractions so they sum to 1, if needed if( isScalingRequired ) { - for( integer ic = 0; ic < numComp-1; ++ic ) + for( integer ic = numSolidComp; ic < numComp-1; ++ic ) { compFrac[ei][ic] /= sum; } @@ -1336,12 +1308,26 @@ void ReactiveCompositionalMultiphaseOBL::updateOBLOperators( ObjectManagerBase & { GEOS_MARK_FUNCTION; - OBLOperatorsKernelFactory:: - createAndLaunch< parallelDevicePolicy<> >( m_numPhases, - m_numComponents, - m_enableEnergyBalance, - dataGroup, - *m_OBLOperatorsTable ); + auto const constitutiveModels = dataGroup.getGroupPointer( ElementSubRegionBase::groupKeyStruct::constitutiveModelsString() ); + + if( constitutiveModels ) + { + for( auto & subGroup : constitutiveModels->getSubGroups() ) + { + if( typeid(*subGroup.second) == typeid(constitutive::OBLFluid) ) + { + auto * oblFluid = dynamic_cast< constitutive::OBLFluid * >( subGroup.second ); + oblFluid->initialize(); + + OBLOperatorsKernelFactory:: + createAndLaunch< parallelDevicePolicy<> >( m_numPhases, + m_numComponents, + m_enableEnergyBalance, + dataGroup, + oblFluid ); + } + } + } } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/ReactiveCompositionalMultiphaseOBL.hpp b/src/coreComponents/physicsSolvers/fluidFlow/ReactiveCompositionalMultiphaseOBL.hpp index c50d111d2e0..0dc9b53c9aa 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/ReactiveCompositionalMultiphaseOBL.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/ReactiveCompositionalMultiphaseOBL.hpp @@ -179,10 +179,16 @@ class ReactiveCompositionalMultiphaseOBL : public FlowSolverBase void updateOBLOperators( ObjectManagerBase & dataGroup ) const; /** - * @brief Get the number of fluid components (species) + * @brief Get the total number of components (species) + * @return the total number of components + */ + localIndex numComponents() const { return m_numComponents; } + + /** + * @brief Get the number of solid components (species) * @return the number of components */ - localIndex numFluidComponents() const { return m_numComponents; } + localIndex numSolidComponents() const { return m_numSolidComponents; } /** * @brief Get the number of fluid phases @@ -235,14 +241,14 @@ class ReactiveCompositionalMultiphaseOBL : public FlowSolverBase static constexpr char const * numComponentsString() { return "numComponents"; } + static constexpr char const * numSolidComponentsString() { return "numSolidComponents"; } + static constexpr char const * enableEnergyBalanceString() { return "enableEnergyBalance"; } static constexpr char const * componentNamesString() { return "componentNames"; } static constexpr char const * phaseNamesString() { return "phaseNames"; } - static constexpr char const * OBLOperatorsTableFileString() { return "OBLOperatorsTableFile"; } - static constexpr char const * transMultExpString() { return "transMultExp"; } static constexpr char const * maxCompFracChangeString() { return "maxCompFractionChange"; } @@ -316,9 +322,12 @@ class ReactiveCompositionalMultiphaseOBL : public FlowSolverBase /// the max number of fluid phases integer m_numPhases; - /// the number of fluid components + /// the number of all components integer m_numComponents; + /// the number of solid components + integer m_numSolidComponents; + /// list of component names string_array m_componentNames; @@ -328,12 +337,6 @@ class ReactiveCompositionalMultiphaseOBL : public FlowSolverBase /// the number of OBL operators integer m_numOBLOperators; - /// 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 = nullptr; - /// flag indicating whether energy balance will be enabled or not integer m_enableEnergyBalance; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/ReactiveCompositionalMultiphaseOBLKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/ReactiveCompositionalMultiphaseOBLKernels.hpp index cb2a651ccd0..bc97492e978 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/ReactiveCompositionalMultiphaseOBLKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/ReactiveCompositionalMultiphaseOBLKernels.hpp @@ -25,8 +25,8 @@ #include "common/GEOS_RAJA_Interface.hpp" #include "constitutive/permeability/PermeabilityBase.hpp" #include "constitutive/permeability/PermeabilityFields.hpp" -#include "functions/MultivariableTableFunction.hpp" -#include "functions/MultivariableTableFunctionKernels.hpp" +#include "functions/MultilinearInterpolatorStaticKernels.hpp" +#include "constitutive/fluid/OBLFluidKernels.hpp" #include "mesh/ElementSubRegionBase.hpp" #include "mesh/ObjectManagerBase.hpp" #include "mesh/utilities/MeshMapUtilities.hpp" @@ -42,251 +42,10 @@ namespace geos namespace reactiveCompositionalMultiphaseOBLKernels { -static constexpr real64 minValueForDivision = 1e-10; - - -// The number of operators in use depends on: -// 1. Number of phases -// 2. Number of components -// 3. Features required in simulation (now its only one - Energy balance) -// This number needs to be used in solver and in kernels (as a template parameter) -// IMHO, this number is too big ( order of 10-100) to be treated via a kernelLaunchSelectorSwitch construct -// Hence, a way to define it once and for all is needed. -// Could be constexpr member of solver, but passing constexpr lambdas require -std=c++17 if I`m not mistaken - -constexpr integer COMPUTE_NUM_OPS ( integer const NP, integer const NC, bool ENERGY ) -{ - auto DOF = NC + ENERGY; - return DOF /*accumulation*/ + - DOF * NP /*flux*/ + - NP /*up_constant*/ + - DOF * NP /*gradient*/ + - DOF /*kinetic rate*/ + - 2 /*rock internal energy and conduction*/ + - 2 * NP /*gravity and capillarity*/ + - 1 /*rock porosity*/ + - 1; -} - - -/******************************** Kernel launch machinery ********************************/ -namespace internal -{ - -template< bool ENABLE_ENERGY, integer NUM_PHASES, typename T, typename LAMBDA > -void kernelLaunchSelectorCompSwitch( T numComps, LAMBDA && lambda ) -{ - static_assert( std::is_integral< T >::value, "kernelLaunchSelectorCompSwitch: type should be integral" ); - - switch( numComps ) - { - case 1: - { lambda( std::integral_constant< T, NUM_PHASES >(), std::integral_constant< T, 1 >(), std::integral_constant< bool, ENABLE_ENERGY >() ); return; } - case 2: - { lambda( std::integral_constant< T, NUM_PHASES >(), std::integral_constant< T, 2 >(), std::integral_constant< bool, ENABLE_ENERGY >() ); return; } - case 3: - { lambda( std::integral_constant< T, NUM_PHASES >(), std::integral_constant< T, 3 >(), std::integral_constant< bool, ENABLE_ENERGY >() ); return; } - case 4: - { lambda( std::integral_constant< T, NUM_PHASES >(), std::integral_constant< T, 4 >(), std::integral_constant< bool, ENABLE_ENERGY >() ); return; } - case 5: - { lambda( std::integral_constant< T, NUM_PHASES >(), std::integral_constant< T, 5 >(), std::integral_constant< bool, ENABLE_ENERGY >()); return; } - default: - { GEOS_ERROR( "Unsupported number of components: " << numComps ); } - } -} - -template< bool ENABLE_ENERGY, typename T, typename LAMBDA > -void kernelLaunchSelectorPhaseSwitch( T numPhases, T numComps, LAMBDA && lambda ) -{ - static_assert( std::is_integral< T >::value, "kernelLaunchSelectorPhaseSwitch: type should be integral" ); - - switch( numPhases ) - { - case 1: - { kernelLaunchSelectorCompSwitch< ENABLE_ENERGY, 1 >( numComps, lambda ); return; } - case 2: - { kernelLaunchSelectorCompSwitch< ENABLE_ENERGY, 2 >( numComps, lambda ); return; } - case 3: - { kernelLaunchSelectorCompSwitch< ENABLE_ENERGY, 3 >( numComps, lambda ); return; } - default: - { GEOS_ERROR( "Unsupported number of phases: " << numPhases ); } - } -} - -template< typename T, typename LAMBDA > -void kernelLaunchSelectorEnergySwitch( T numPhases, T numComps, bool enableEnergyBalance, LAMBDA && lambda ) -{ - if( enableEnergyBalance ) - { - kernelLaunchSelectorPhaseSwitch< true >( numPhases, numComps, lambda ); - } - else - { - kernelLaunchSelectorPhaseSwitch< false >( numPhases, numComps, lambda ); - } -} +using namespace oblFluidKernels; -} // namespace internal - - -/******************************** OBLOperatorsKernel ********************************/ - -/** - * @class OBLOperatorsKernel - * @tparam NUM_PHASES number of phases - * @tparam NUM_COMPS number of components - * @tparam ENABLE_ENERGY flag if energy balance equation is assembled - * @brief Compute OBL Operators and derivatives - */ -template< integer NUM_PHASES, integer NUM_COMPS, bool ENABLE_ENERGY > -class OBLOperatorsKernel -{ -public: - /// Compile time value for the energy balance switch - static constexpr integer enableEnergyBalance = ENABLE_ENERGY; - /// Compile time value for the number of components - static constexpr integer numComps = NUM_COMPS; - - /// Compile time value for the number of dimensions - static constexpr integer numDofs = numComps + enableEnergyBalance; - // /// Compile time value for the number of operators - static constexpr integer numOps = COMPUTE_NUM_OPS( NUM_PHASES, NUM_COMPS, ENABLE_ENERGY ); - - static constexpr real64 barToPascalMult = 1e5; - static constexpr real64 pascalToBarMult = 1.0 / 1e5; - - /** - * @brief Performs the kernel launch - * @tparam POLICY the policy used in the RAJA kernels - * @tparam KERNEL_TYPE the kernel type - * @param[in] numElems the number of elements - * @param[inout] kernelComponent the kernel component providing access to the compute function - */ - template< typename POLICY, typename KERNEL_TYPE > - static void - launch( localIndex const numElems, - KERNEL_TYPE const & kernelComponent ) - { - forAll< POLICY >( numElems, [=] GEOS_HOST_DEVICE ( localIndex const ei ) - { - kernelComponent.compute( ei ); - } ); - } - - /** - * @brief Constructor - * @param[in] subRegion the element subregion - * @param[in] OBLOperatorsTable the OBL table function kernel - */ - OBLOperatorsKernel( ObjectManagerBase & subRegion, - MultivariableTableFunctionStaticKernel< numDofs, numOps > OBLOperatorsTable ) - : - m_OBLOperatorsTable( OBLOperatorsTable ), - m_pressure( subRegion.getField< fields::flow::pressure >() ), - m_compFrac( subRegion.getField< fields::flow::globalCompFraction >() ), - m_temperature( subRegion.getField< fields::flow::temperature >() ), - m_OBLOperatorValues ( subRegion.getField< fields::flow::OBLOperatorValues >()), - m_OBLOperatorDerivatives ( subRegion.getField< fields::flow::OBLOperatorDerivatives >()) - {} - - /** - * @brief Compute the operator values and derivatives for an element - * @param[in] ei the element index - */ - GEOS_HOST_DEVICE - inline - void compute( localIndex const ei ) const - { - arraySlice1d< real64 const, compflow::USD_COMP - 1 > const compFrac = m_compFrac[ei]; - arraySlice1d< real64, compflow::USD_OBL_VAL - 1 > const & OBLVals = m_OBLOperatorValues[ei]; - arraySlice2d< real64, compflow::USD_OBL_DER - 1 > const & OBLDers = m_OBLOperatorDerivatives[ei]; - real64 state[numDofs]; - - // we need to convert pressure from Pa (internal unit in GEOSX) to bar (internal unit in DARTS) - state[0] = m_pressure[ei] * pascalToBarMult; - - // the last component fraction is not used to define the state - for( integer i = 1; i < numComps; ++i ) - { - state[i] = compFrac[i - 1]; - } - - if( enableEnergyBalance ) - { - state[numDofs - 1] = m_temperature[ei]; - } - - m_OBLOperatorsTable.compute( state, OBLVals, OBLDers ); - - // we do not perform derivatives unit conversion here: - // instead we postpone it till all the derivatives are fully formed, and only then apply the factor only once in 'complete' function - // scaling the whole system might be even better solution (every pressure column needs to be multiplied by pascalToBarMult) - } - -private: - - // inputs - MultivariableTableFunctionStaticKernel< numDofs, numOps > m_OBLOperatorsTable; - - // Views on primary variables and their updates - arrayView1d< real64 const > m_pressure; - arrayView2d< real64 const, compflow::USD_COMP > m_compFrac; - arrayView1d< real64 const > m_temperature; - - // outputs - - // Views on OBL operator values and derivatives - arrayView2d< real64, compflow::USD_OBL_VAL > m_OBLOperatorValues; - arrayView3d< real64, compflow::USD_OBL_DER > m_OBLOperatorDerivatives; -}; - -/** - * @class OBLOperatorsKernelFactory - */ -class OBLOperatorsKernelFactory -{ -public: - - /** - * @brief Create a new kernel and launch - * @tparam POLICY the policy used in the RAJA kernel - * @param[in] numPhases the number of phases - * @param[in] numComponents the number of components - * @param[in] enableEnergyBalance flag if energy balance equation is assembled - * @param[in] subRegion the element subregion - * @param[in] function the OBL table function - */ - template< typename POLICY > - static void - createAndLaunch( integer const numPhases, - integer const numComponents, - bool const enableEnergyBalance, - ObjectManagerBase & subRegion, - MultivariableTableFunction const & function ) - { - internal::kernelLaunchSelectorEnergySwitch( numPhases, numComponents, enableEnergyBalance, [&] ( auto NP, auto NC, auto E ) - { - integer constexpr ENABLE_ENERGY = E(); - integer constexpr NUM_PHASES = NP(); - integer constexpr NUM_COMPS = NC(); - integer constexpr NUM_DIMS = ENABLE_ENERGY + NUM_COMPS; - integer constexpr NUM_OPS = COMPUTE_NUM_OPS( NUM_PHASES, NUM_COMPS, ENABLE_ENERGY ); - - OBLOperatorsKernel< NUM_PHASES, NUM_COMPS, ENABLE_ENERGY > - kernel( subRegion, - MultivariableTableFunctionStaticKernel< NUM_DIMS, NUM_OPS >( function.getAxisMinimums(), - function.getAxisMaximums(), - function.getAxisPoints(), - function.getAxisSteps(), - function.getAxisStepInvs(), - function.getAxisHypercubeMults(), - function.getHypercubeData() - ) ); - OBLOperatorsKernel< NUM_PHASES, NUM_COMPS, ENABLE_ENERGY >::template launch< POLICY >( subRegion.size(), kernel ); - } ); - } +static constexpr real64 minValueForDivision = 1e-10; -}; /******************************** ElementBasedAssemblyKernel ********************************/ @@ -318,8 +77,8 @@ class ElementBasedAssemblyKernel static constexpr integer ACC_OP = 0; // kinetic reaction static constexpr integer KIN_OP = numDofs + numDofs * numPhases + numPhases + numDofs * numPhases; - // rock internal energy - static constexpr integer RE_INTER_OP = numDofs + numDofs * numPhases + numPhases + numDofs * numPhases + numDofs; + // temperature + static constexpr integer TEMP_OP = KIN_OP + numDofs + 2 * numPhases + 1 + numPhases; static constexpr real64 secondsToDaysMult = 1.0 / (60 * 60 * 24); static constexpr real64 pascalToBarMult = 1.0 / 1e5; @@ -443,11 +202,11 @@ class ElementBasedAssemblyKernel // + rock energy if( enableEnergyBalance ) { - stack.localResidual[numDofs-1] += m_referenceRockVolume[ei] * (OBLVals[RE_INTER_OP] - OBLVals_n[RE_INTER_OP]) * m_rockVolumetricHeatCapacity[ei]; + stack.localResidual[numDofs-1] += m_referenceRockVolume[ei] * (OBLVals[TEMP_OP] - OBLVals_n[TEMP_OP]) * m_rockVolumetricHeatCapacity[ei]; for( integer v = 0; v < numDofs; ++v ) { - stack.localJacobian[numDofs-1][v] += m_referenceRockVolume[ei] * OBLDers[RE_INTER_OP][v] * m_rockVolumetricHeatCapacity[ei]; + stack.localJacobian[numDofs-1][v] += m_referenceRockVolume[ei] * OBLDers[TEMP_OP][v] * m_rockVolumetricHeatCapacity[ei]; } // end of fill offdiagonal part + contribute to diagonal } } @@ -571,7 +330,7 @@ class ElementBasedAssemblyKernelFactory CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) { - internal::kernelLaunchSelectorEnergySwitch( numPhases, numComps, enableEnergyBalance, [&] ( auto NP, auto NC, auto E ) + oblFluidKernels::internal::kernelLaunchSelectorEnergySwitch( numPhases, numComps, enableEnergyBalance, [&] ( auto NP, auto NC, auto E ) { integer constexpr ENABLE_ENERGY = E(); integer constexpr NUM_PHASES = NP(); @@ -622,6 +381,9 @@ class FluxComputeKernelBase static constexpr real64 transUnitMult = 8639693349945.239; static constexpr real64 transDUnitMult = 0.00852671467191601; + // tolerance of phase existence check in the calculation of diffusion flux + static constexpr real64 phaseExistenceTolerance = 1.e-6; + /** * @brief The type for element-based data. Consists entirely of ArrayView's. * @@ -758,18 +520,16 @@ class FluxComputeKernel : public FluxComputeKernelBase static constexpr integer FLUX_OP = numDofs; // diffusion static constexpr integer UPSAT_OP = numDofs + numDofs * numPhases; + // diffusion static constexpr integer GRAD_OP = numDofs + numDofs * numPhases + numPhases; - - // temperature - static constexpr integer RE_TEMP_OP = numDofs + numDofs * numPhases + numPhases + numDofs * numPhases + numDofs + 1; - // rock conduction - static constexpr integer ROCK_COND_OP = numDofs + numDofs * numPhases + numPhases + numDofs * numPhases + numDofs + 2; // gravity (density) - static constexpr integer GRAV_OP = numDofs + numDofs * numPhases + numPhases + numDofs * numPhases + numDofs + 3; + static constexpr integer GRAV_OP = numDofs + numDofs * numPhases + numPhases + numDofs * numPhases + numDofs; // capillary pressure - static constexpr integer PC_OP = numDofs + numDofs * numPhases + numPhases + numDofs * numPhases + numDofs + 3 + numPhases; + static constexpr integer PC_OP = numDofs + numDofs * numPhases + numPhases + numDofs * numPhases + numDofs + numPhases; // porosity - static constexpr integer PORO_OP = numDofs + numDofs * numPhases + numPhases + numDofs * numPhases + numDofs + 3 + 2 * numPhases; + static constexpr integer PORO_OP = numDofs + numDofs * numPhases + numPhases + numDofs * numPhases + numDofs + 2 * numPhases; + // temperature + static constexpr integer TEMP_OP = numDofs + numDofs * numPhases + numPhases + numDofs * numPhases + numDofs + 3 * numPhases + 1; /// Maximum number of elements at the face static constexpr localIndex maxNumElems = STENCILWRAPPER::maxNumPointsInFlux; @@ -935,6 +695,7 @@ class FluxComputeKernel : public FluxComputeKernelBase real64 transMult = 1; real64 transMultD = 0; + real64 phasePresenceMult = 1.0; if( m_transMultExp > 0 ) { // Calculate transmissibility multiplier: @@ -1005,36 +766,35 @@ class FluxComputeKernel : public FluxComputeKernelBase } // end of loop over number of phases for convective operator with gravity and capillarity - // [3] Additional diffusion code here: (phi_p * S_p) * (rho_p * D_cp * Delta_x_cp) or (phi_p * S_p) * (kappa_p * Delta_T) - - real64 const poroAverage = (m_referencePorosity[erI][esrI][eiI] + m_referencePorosity[erJ][esrJ][eiJ]) * 0.5; // diffusion term - // depends on total - // porosity! - - - // Add diffusion term to the residual: + // [3] Diffusion code here + // TODO: make sure this works only for matrix-matrix connections, avoiding wells for( integer c = 0; c < numDofs; ++c ) { for( integer p = 0; p < numPhases; ++p ) { - real64 phaseCompGrad = OBLValsJ[GRAD_OP + c * numPhases + p] - OBLValsI[GRAD_OP + c * numPhases + p]; + // mass concentration gradient + real64 phaseCompGrad = OBLValsJ[GRAD_OP + p * numDofs + c] - OBLValsI[GRAD_OP + p * numDofs + c]; - arraySlice1d< real64 const, compflow::USD_OBL_VAL - 1 > const & OBLValsUp = (phaseCompGrad < 0) ? OBLValsI : OBLValsJ; - arraySlice2d< real64 const, compflow::USD_OBL_DER - 1 > const & OBLDersUp = (phaseCompGrad < 0) ? OBLDersI : OBLDersJ; - integer const upOffset = (phaseCompGrad < 0) ? 0 : numDofs; + // phase existence check + if( OBLValsI[UPSAT_OP + p] * OBLValsJ[UPSAT_OP + p] > phaseExistenceTolerance ) + phasePresenceMult = 1.0; + else + phasePresenceMult = 0.0; - // use upstream quantity from cell i for compressibility and saturation (mass or energy): - real64 phaseGammaDDiff = m_dt * transD * poroAverage * OBLValsUp[UPSAT_OP + p]; + // use arithmetic mean for multiplier + real64 phaseGammaDDiff = m_dt * phasePresenceMult * transD * (m_referencePorosity[erI][esrI][eiI] * OBLValsI[UPSAT_OP + p] + + m_referencePorosity[erJ][esrJ][eiJ] * OBLValsJ[UPSAT_OP + p]) * 0.5; stack.localFlux[c] -= phaseGammaDDiff * phaseCompGrad; // diffusion term // Add diffusion terms to Jacobian: for( integer v = 0; v < numDofs; ++v ) { - stack.localFluxJacobian[c][v] += phaseGammaDDiff * OBLDersI[GRAD_OP + c * numPhases + p][v]; - stack.localFluxJacobian[c][numDofs + v] -= phaseGammaDDiff * OBLDersJ[GRAD_OP + c * numPhases + p][v]; + stack.localFluxJacobian[c][v] += phaseGammaDDiff * OBLDersI[GRAD_OP + p * numDofs + c][v]; + stack.localFluxJacobian[c][numDofs + v] -= phaseGammaDDiff * OBLDersJ[GRAD_OP + p * numDofs + c][v]; - stack.localFluxJacobian[c][v + upOffset] -= phaseCompGrad * m_dt * transD * poroAverage * OBLDersUp[UPSAT_OP + p][v]; + stack.localFluxJacobian[c][v] -= phaseCompGrad * m_dt * phasePresenceMult * transD * m_referencePorosity[erI][esrI][eiI] * OBLDersI[UPSAT_OP + p][v] * 0.5; + stack.localFluxJacobian[c][numDofs + v] -= phaseCompGrad * m_dt * phasePresenceMult * transD * m_referencePorosity[erJ][esrJ][eiJ] * OBLDersJ[UPSAT_OP + p][v] * 0.5; } } @@ -1043,31 +803,16 @@ class FluxComputeKernel : public FluxComputeKernelBase // [4] add rock conduction if( enableEnergyBalance ) { - real64 const tDiff = OBLValsJ[RE_TEMP_OP] - OBLValsI[RE_TEMP_OP]; - real64 const gammaTDiff = transD * m_dt * tDiff; - - arraySlice1d< real64 const, compflow::USD_OBL_VAL - 1 > const & OBLValsUp = (tDiff < 0) ? OBLValsI : OBLValsJ; - arraySlice2d< real64 const, compflow::USD_OBL_DER - 1 > const & OBLDersUp = (tDiff < 0) ? OBLDersI : OBLDersJ; - integer const upOffset = (tDiff < 0) ? 0 : numDofs; - localIndex const erUp = (tDiff < 0) ? erI : erJ; - localIndex const esrUp = (tDiff < 0) ? esrI : esrJ; - localIndex const eiUp = (tDiff < 0) ? eiI : eiJ; - - real64 const poroUp = m_referencePorosity[erUp][esrUp][eiUp]; - real64 const rockCondUp = m_rockThermalConductivity[erUp][esrUp][eiUp]; + real64 const tDiff = OBLValsJ[TEMP_OP] - OBLValsI[TEMP_OP]; + real64 const gammaTi = transD * m_dt * (1 - m_referencePorosity[erI][esrI][eiI]) * m_rockThermalConductivity[erI][esrI][eiI]; + real64 const gammaTj = transD * m_dt * (1 - m_referencePorosity[erJ][esrJ][eiJ]) * m_rockThermalConductivity[erJ][esrJ][eiJ]; // rock heat transfers flows from cell i to j - stack.localFlux[numComps] -= gammaTDiff * OBLValsUp[ROCK_COND_OP] * (1 - poroUp) * rockCondUp; + stack.localFlux[numComps] -= tDiff * (gammaTi + gammaTj) / 2; for( integer v = 0; v < numDofs; ++v ) { - stack.localFluxJacobian[numComps][v + upOffset] -= gammaTDiff * OBLDersUp[ROCK_COND_OP][v] * (1 - poroUp) * rockCondUp; - // the last variable - T - if( v == numDofs - 1 ) - { - real64 const tFlux = transD * m_dt * OBLValsUp[ROCK_COND_OP] * (1 - poroUp) * rockCondUp; - stack.localFluxJacobian[numComps][v] += tFlux; - stack.localFluxJacobian[numComps][numDofs + v] -= tFlux; - } + stack.localFluxJacobian[numComps][v] += OBLDersI[TEMP_OP][v] * (gammaTi + gammaTj) / 2; + stack.localFluxJacobian[numComps][numDofs + v] -= OBLDersJ[TEMP_OP][v] * (gammaTi + gammaTj) / 2; } } } @@ -1208,7 +953,7 @@ class FluxComputeKernelFactory CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) { - internal::kernelLaunchSelectorEnergySwitch( numPhases, numComps, enableEnergyBalance, [&] ( auto NP, auto NC, auto E ) + oblFluidKernels::internal::kernelLaunchSelectorEnergySwitch( numPhases, numComps, enableEnergyBalance, [&] ( auto NP, auto NC, auto E ) { integer constexpr ENABLE_ENERGY = E(); integer constexpr NUM_PHASES = NP(); diff --git a/src/coreComponents/schema/schema.xsd b/src/coreComponents/schema/schema.xsd index 64de5523946..1f845342039 100644 --- a/src/coreComponents/schema/schema.xsd +++ b/src/coreComponents/schema/schema.xsd @@ -807,6 +807,10 @@ + + + + @@ -4345,8 +4349,6 @@ When set to 2 output convergence information to a csv--> - - @@ -4389,6 +4391,8 @@ Information output from lower logLevels is added with the desired log level + + @@ -5918,6 +5922,7 @@ Information output from lower logLevels is added with the desired log level + @@ -7059,6 +7064,16 @@ If you want to do a three-phase simulation, please use instead wettingIntermedia + + + + + + + + + + diff --git a/src/coreComponents/schema/schema.xsd.other b/src/coreComponents/schema/schema.xsd.other index e3d4e22da8a..cbd7798327e 100644 --- a/src/coreComponents/schema/schema.xsd.other +++ b/src/coreComponents/schema/schema.xsd.other @@ -343,12 +343,14 @@ + + @@ -1534,6 +1536,7 @@ + @@ -2668,6 +2671,7 @@ + diff --git a/src/coreComponents/unitTests/fluidFlowTests/testReactiveCompositionalMultiphaseOBL.cpp b/src/coreComponents/unitTests/fluidFlowTests/testReactiveCompositionalMultiphaseOBL.cpp index 326ffc92140..3eef251df49 100644 --- a/src/coreComponents/unitTests/fluidFlowTests/testReactiveCompositionalMultiphaseOBL.cpp +++ b/src/coreComponents/unitTests/fluidFlowTests/testReactiveCompositionalMultiphaseOBL.cpp @@ -194,7 +194,7 @@ void testOperatorsNumericalDerivatives( ReactiveCompositionalMultiphaseOBL & sol real64 const perturbParameter, real64 const relTol ) { - localIndex const NC = solver.numFluidComponents(); + localIndex const NC = solver.numComponents(); localIndex const NOPS = solver.numOBLOperators(); string_array operators( NOPS ); @@ -322,7 +322,7 @@ void testNumericalJacobian( ReactiveCompositionalMultiphaseOBL & solver, real64 const relTol, LAMBDA assembleFunction ) { - localIndex const NC = solver.numFluidComponents(); + localIndex const NC = solver.numComponents(); CRSMatrix< real64, globalIndex > const & jacobian = solver.getLocalMatrix(); array1d< real64 > residual( jacobian.numRows() ); diff --git a/src/pygeosx/pygeosx.cpp b/src/pygeosx/pygeosx.cpp index fff9d7aea9c..3ae79fdd5be 100644 --- a/src/pygeosx/pygeosx.cpp +++ b/src/pygeosx/pygeosx.cpp @@ -25,6 +25,7 @@ #include "fileIO/python/PyHistoryCollectionType.hpp" #include "fileIO/python/PyHistoryOutputType.hpp" #include "fileIO/python/PyVTKOutputType.hpp" +#include "functions/python/PyPythonFunctionType.hpp" #include "mainInterface/initialization.hpp" #include "LvArray/src/python/PyArray.hpp" @@ -430,6 +431,11 @@ PyInit_pygeosx() return nullptr; } + if( !LvArray::python::addTypeToModule( module, geos::python::getPyPythonFunctionType(), "PythonFunction" ) ) + { + return nullptr; + } + // Add the LvArray submodule. if( !LvArray::python::addPyLvArrayModule( module ) ) {