Bugfixes for CUDA build (multiphase model kernel wrappers, HypreVector::extract)

src/CMakeLists.txt

@@ -96,6 +96,11 @@ blt_add_executable(NAME geosx
                               ${externalComponentsLinkList}
                   )
 
+# Seems to be required on some CMake versions (e.g. 3.16) to get enforce device linking
+if( ENABLE_HYPRE_CUDA )
+  set_target_properties( geosx PROPERTIES CUDA_RESOLVE_DEVICE_SYMBOLS TRUE )
+endif()
+
 # Removing all transitive link dependencies from geosx_core target to circumvent
 # the BLT behavior which imposes all dependencies to be public
 #set_target_properties(geosx_core PROPERTIES INTERFACE_LINK_LIBRARIES "")

src/coreComponents/constitutive/fluid/PVTFunctions/BrineCO2Density.cpp

@@ -123,8 +123,7 @@ void BrineCO2Density::calculateBrineDensity( PTTableCoordinates const & tableCoo
 
 BrineCO2Density::KernelWrapper BrineCO2Density::createKernelWrapper()
 {
-  return KernelWrapper( m_componentNames,
-                        m_componentMolarWeight,
+  return KernelWrapper( m_componentMolarWeight,
                         m_brineDensityTable,
                         m_CO2Index,
                         m_waterIndex );

src/coreComponents/constitutive/fluid/PVTFunctions/BrineCO2Density.hpp

@@ -37,13 +37,11 @@ class BrineCO2DensityUpdate final : public PVTFunctionBaseUpdate
 {
 public:
 
-  BrineCO2DensityUpdate( arrayView1d< string const > const & componentNames,
-                         arrayView1d< real64 const > const & componentMolarWeight,
+  BrineCO2DensityUpdate( arrayView1d< real64 const > const & componentMolarWeight,
                          TableFunction * brineDensityTable,
                          localIndex const CO2Index,
                          localIndex const waterIndex )
-    : PVTFunctionBaseUpdate( componentNames,
-                             componentMolarWeight ),
+    : PVTFunctionBaseUpdate( componentMolarWeight ),
     m_brineDensityTable( brineDensityTable->createKernelWrapper() ),
     m_CO2Index( CO2Index ),
     m_waterIndex( waterIndex )
@@ -73,12 +71,18 @@ class BrineCO2DensityUpdate final : public PVTFunctionBaseUpdate
                         real64 & dValue_dPressure,
                         real64 & dValue_dTemperature,
                         arraySlice1d< real64 > const & dValue_dGlobalCompFraction,
-                        bool useMass = 0 ) const override;
+                        bool useMass ) const override;
+
+  virtual void move( LvArray::MemorySpace const space, bool const touch ) override
+  {
+    PVTFunctionBaseUpdate::move( space, touch );
+    m_brineDensityTable.move( space, touch );
+  }
 
 protected:
 
   /// Table with brine density tabulated as a function (P,T,sal)
-  TableFunction::KernelWrapper const m_brineDensityTable;
+  TableFunction::KernelWrapper m_brineDensityTable;
 
   /// Index of the CO2 component
   localIndex const m_CO2Index;

src/coreComponents/constitutive/fluid/PVTFunctions/BrineViscosity.cpp

@@ -68,8 +68,7 @@ void BrineViscosity::makeCoefficients( string_array const & inputPara )
 
 BrineViscosity::KernelWrapper BrineViscosity::createKernelWrapper()
 {
-  return KernelWrapper( m_componentNames,
-                        m_componentMolarWeight,
+  return KernelWrapper( m_componentMolarWeight,
                         m_coef0,
                         m_coef1 );
 }

src/coreComponents/constitutive/fluid/PVTFunctions/BrineViscosity.hpp

@@ -34,12 +34,10 @@ class BrineViscosityUpdate final : public PVTFunctionBaseUpdate
 {
 public:
 
-  BrineViscosityUpdate( arrayView1d< string const > const & componentNames,
-                        arrayView1d< real64 const > const & componentMolarWeight,
+  BrineViscosityUpdate( arrayView1d< real64 const > const & componentMolarWeight,
                         real64 const & coef0,
                         real64 const & coef1 )
-    : PVTFunctionBaseUpdate( componentNames,
-                             componentMolarWeight ),
+    : PVTFunctionBaseUpdate( componentMolarWeight ),
     m_coef0( coef0 ),
     m_coef1( coef1 )
   {}
@@ -68,7 +66,7 @@ class BrineViscosityUpdate final : public PVTFunctionBaseUpdate
                         real64 & dValue_dPressure,
                         real64 & dValue_dTemperature,
                         arraySlice1d< real64 > const & dValue_dGlobalCompFraction,
-                        bool useMass = 0 ) const override;
+                        bool useMass ) const override;
 
 protected:
 

src/coreComponents/constitutive/fluid/PVTFunctions/CO2Solubility.cpp

@@ -291,8 +291,7 @@ void CO2Solubility::CO2SolubilityFunction( real64 const & tolerance,
 
 CO2Solubility::KernelWrapper CO2Solubility::createKernelWrapper()
 {
-  return KernelWrapper( m_componentNames,
-                        m_componentMolarWeight,
+  return KernelWrapper( m_componentMolarWeight,
                         m_CO2SolubilityTable,
                         m_CO2Index,
                         m_waterIndex,

src/coreComponents/constitutive/fluid/PVTFunctions/CO2Solubility.hpp

@@ -39,15 +39,13 @@ class CO2SolubilityUpdate final : public FlashModelBaseUpdate
 {
 public:
 
-  CO2SolubilityUpdate( arrayView1d< string const > const & componentNames,
-                       arrayView1d< real64 const > const & componentMolarWeight,
+  CO2SolubilityUpdate( arrayView1d< real64 const > const & componentMolarWeight,
                        TableFunction * CO2SolubilityTable,
                        localIndex const CO2Index,
                        localIndex const waterIndex,
                        localIndex const phaseGasIndex,
                        localIndex const phaseLiquidIndex )
-    : FlashModelBaseUpdate( componentNames,
-                            componentMolarWeight ),
+    : FlashModelBaseUpdate( componentMolarWeight ),
     m_CO2SolubilityTable( CO2SolubilityTable->createKernelWrapper() ),
     m_CO2Index( CO2Index ),
     m_waterIndex( waterIndex ),
@@ -79,12 +77,18 @@ class CO2SolubilityUpdate final : public FlashModelBaseUpdate
                         arraySlice2d< real64 > const & phaseCompFraction,
                         arraySlice2d< real64 > const & dPhaseCompFraction_dPressure,
                         arraySlice2d< real64 > const & dPhaseCompFraction_dTemperature,
-                        arraySlice3d< real64 > const & dPhaseCompFraction_dCompFraction ) const;
+                        arraySlice3d< real64 > const & dPhaseCompFraction_dCompFraction ) const override;
+
+  virtual void move( LvArray::MemorySpace const space, bool const touch ) override
+  {
+    FlashModelBaseUpdate::move( space, touch );
+    m_CO2SolubilityTable.move( space, touch );
+  }
 
 protected:
 
   /// Table with CO2 solubility tabulated as a function (P,T)
-  TableFunction::KernelWrapper const m_CO2SolubilityTable;
+  TableFunction::KernelWrapper m_CO2SolubilityTable;
 
   /// Index of the CO2 phase
   localIndex const m_CO2Index;

src/coreComponents/constitutive/fluid/PVTFunctions/FenghourCO2Viscosity.cpp

@@ -131,8 +131,7 @@ void FenghourCO2Viscosity::calculateCO2Viscosity( PTTableCoordinates const & tab
 
 FenghourCO2Viscosity::KernelWrapper FenghourCO2Viscosity::createKernelWrapper()
 {
-  return KernelWrapper( m_componentNames,
-                        m_componentMolarWeight,
+  return KernelWrapper( m_componentMolarWeight,
                         m_CO2ViscosityTable );
 }
 

src/coreComponents/constitutive/fluid/PVTFunctions/FenghourCO2Viscosity.hpp

@@ -37,11 +37,9 @@ class FenghourCO2ViscosityUpdate final : public PVTFunctionBaseUpdate
 {
 public:
 
-  FenghourCO2ViscosityUpdate( arrayView1d< string const > const & componentNames,
-                              arrayView1d< real64 const > const & componentMolarWeight,
+  FenghourCO2ViscosityUpdate( arrayView1d< real64 const > const & componentMolarWeight,
                               TableFunction * CO2ViscosityTable )
-    : PVTFunctionBaseUpdate( componentNames,
-                             componentMolarWeight ),
+    : PVTFunctionBaseUpdate( componentMolarWeight ),
     m_CO2ViscosityTable( CO2ViscosityTable->createKernelWrapper() )
   {}
 
@@ -69,12 +67,18 @@ class FenghourCO2ViscosityUpdate final : public PVTFunctionBaseUpdate
                         real64 & dValue_dPressure,
                         real64 & dValue_dTemperature,
                         arraySlice1d< real64 > const & dValue_dGlobalCompFraction,
-                        bool useMass = 0 ) const override;
+                        bool useMass ) const override;
+
+  virtual void move( LvArray::MemorySpace const space, bool const touch ) override
+  {
+    PVTFunctionBaseUpdate::move( space, touch );
+    m_CO2ViscosityTable.move( space, touch );
+  }
 
 protected:
 
   /// Table with viscosity tabulated as a function (P,T)
-  TableFunction::KernelWrapper const m_CO2ViscosityTable;
+  TableFunction::KernelWrapper m_CO2ViscosityTable;
 
 };
 

src/coreComponents/constitutive/fluid/PVTFunctions/FlashModelBase.hpp

@@ -35,10 +35,8 @@ class FlashModelBaseUpdate
 {
 public:
 
-  FlashModelBaseUpdate( arrayView1d< string const > const & componentNames,
-                        arrayView1d< real64 const > const & componentMolarWeight )
+  FlashModelBaseUpdate( arrayView1d< real64 const > const & componentMolarWeight )
     :
-    m_componentNames( componentNames ),
     m_componentMolarWeight( componentMolarWeight )
   {}
 
@@ -70,10 +68,19 @@ class FlashModelBaseUpdate
                         arraySlice2d< real64 > const & dPhaseCompFraction_dTemperature,
                         arraySlice3d< real64 > const & dPhaseCompFraction_dCompFraction ) const = 0;
 
-protected:
+  /**
+   * @brief Move the KernelWrapper to the given execution space, optionally touching it.
+   * @param space the space to move the KernelWrapper to
+   * @param touch whether the KernelWrapper should be touched in the new space or not
+   * @note This function exists to enable holding KernelWrapper objects in an ArrayView
+   *       and have their contents properly moved between memory spaces.
+   */
+  virtual void move( LvArray::MemorySpace const space, bool const touch )
+  {
+    m_componentMolarWeight.move( space, touch );
+  }
 
-  /// Array storing the name of the components
-  arrayView1d< string const > m_componentNames;
+protected:
 
   /// Array storing the component molar weights
   arrayView1d< real64 const > m_componentMolarWeight;

src/coreComponents/constitutive/fluid/PVTFunctions/PVTFunctionBase.hpp

@@ -36,10 +36,8 @@ enum class PVTFunctionType { UNKNOWN, DENSITY, VISCOSITY };
 class PVTFunctionBaseUpdate
 {
 public:
-  PVTFunctionBaseUpdate( arrayView1d< string const > const & componentNames,
-                         arrayView1d< real64 const > const & componentMolarWeight )
+  PVTFunctionBaseUpdate( arrayView1d< real64 const > const & componentMolarWeight )
     :
-    m_componentNames( componentNames ),
     m_componentMolarWeight( componentMolarWeight )
   {}
 
@@ -69,13 +67,22 @@ class PVTFunctionBaseUpdate
                         real64 & dValue_dPressure,
                         real64 & dValue_dTemperature,
                         arraySlice1d< real64 > const & dValue_dGlobalCompFraction,
-                        bool useMass = 0 ) const = 0;
+                        bool useMass ) const = 0;
+
+  /**
+   * @brief Move the KernelWrapper to the given execution space, optionally touching it.
+   * @param space the space to move the KernelWrapper to
+   * @param touch whether the KernelWrapper should be touched in the new space or not
+   * @note This function exists to enable holding KernelWrapper objects in an ArrayView
+   *       and have their contents properly moved between memory spaces.
+   */
+  virtual void move( LvArray::MemorySpace const space, bool const touch )
+  {
+    m_componentMolarWeight.move( space, touch );
+  }
 
 protected:
 
-  /// Array storing the name of the components
-  arrayView1d< string const > m_componentNames;
-
   /// Array storing the component molar weights
   arrayView1d< real64 const > m_componentMolarWeight;
 

src/coreComponents/constitutive/fluid/PVTFunctions/SpanWagnerCO2Density.cpp

@@ -278,8 +278,7 @@ void SpanWagnerCO2Density::spanWagnerCO2DensityFunction( real64 const & toleranc
 
 SpanWagnerCO2Density::KernelWrapper SpanWagnerCO2Density::createKernelWrapper()
 {
-  return KernelWrapper( m_componentNames,
-                        m_componentMolarWeight,
+  return KernelWrapper( m_componentMolarWeight,
                         m_CO2DensityTable,
                         m_CO2Index );
 }

src/coreComponents/constitutive/fluid/PVTFunctions/SpanWagnerCO2Density.hpp

@@ -37,12 +37,10 @@ class SpanWagnerCO2DensityUpdate final : public PVTFunctionBaseUpdate
 {
 public:
 
-  SpanWagnerCO2DensityUpdate( arrayView1d< string const > const & componentNames,
-                              arrayView1d< real64 const > const & componentMolarWeight,
+  SpanWagnerCO2DensityUpdate( arrayView1d< real64 const > const & componentMolarWeight,
                               TableFunction * CO2DensityTable,
                               localIndex const CO2Index )
-    : PVTFunctionBaseUpdate( componentNames,
-                             componentMolarWeight ),
+    : PVTFunctionBaseUpdate( componentMolarWeight ),
     m_CO2DensityTable( CO2DensityTable->createKernelWrapper() ),
     m_CO2Index( CO2Index )
   {}
@@ -71,12 +69,18 @@ class SpanWagnerCO2DensityUpdate final : public PVTFunctionBaseUpdate
                         real64 & dValue_dPressure,
                         real64 & dValue_dTemperature,
                         arraySlice1d< real64 > const & dValue_dGlobalCompFraction,
-                        bool useMass = 0 ) const override;
+                        bool useMass ) const override;
+
+  virtual void move( LvArray::MemorySpace const space, bool const touch ) override
+  {
+    PVTFunctionBaseUpdate::move( space, touch );
+    m_CO2DensityTable.move( space, touch );
+  }
 
 protected:
 
   /// Table with brine density tabulated as a function (P,T,sal)
-  TableFunction::KernelWrapper const m_CO2DensityTable;
+  TableFunction::KernelWrapper m_CO2DensityTable;
 
   /// Index of the CO2 component
   localIndex const m_CO2Index;

src/coreComponents/functions/TableFunction.hpp

@@ -126,6 +126,21 @@ class TableFunction : public FunctionBase
     void
     compute( IN_ARRAY const & input, real64 & value, OUT_ARRAY && derivatives ) const;
 
+    /**
+     * @brief Move the KernelWrapper to the given execution space, optionally touching it.
+     * @param space the space to move the KernelWrapper to
+     * @param touch whether the KernelWrapper should be touched in the new space or not
+     * @note This function exists to enable holding KernelWrapper objects in an ArrayView
+     *       and have their contents properly moved between memory spaces.
+     */
+    void move( LvArray::MemorySpace const space, bool const touch )
+    {
+      m_coordinates.move( space, touch );
+      m_values.move( space, touch );
+      m_size.move( space, touch );
+      m_indexIncrement.move( space, touch );
+    }
+
 private:
 
     /// Table interpolation method

src/coreComponents/linearAlgebra/interfaces/hypre/HypreMatrix.cpp

@@ -787,7 +787,7 @@ void HypreMatrix::scale( real64 const scalingFactor )
   HYPRE_Int const diag_nnz = hypre_CSRMatrixNumNonzeros( prt_diag_CSR );
   HYPRE_Real * const ptr_diag_data = hypre_CSRMatrixData( prt_diag_CSR );
 
-  forAll< execPolicy >( diag_nnz, [=] GEOSX_HOST_DEVICE ( HYPRE_Int const i )
+  forAll< hypre::execPolicy >( diag_nnz, [=] GEOSX_HYPRE_HOST_DEVICE ( HYPRE_Int const i )
   {
     ptr_diag_data[i] *= scalingFactor;
   } );
@@ -796,7 +796,7 @@ void HypreMatrix::scale( real64 const scalingFactor )
   HYPRE_Int const offdiag_nnz = hypre_CSRMatrixNumNonzeros( prt_offdiag_CSR );
   HYPRE_Real * const ptr_offdiag_data = hypre_CSRMatrixData( prt_offdiag_CSR );
 
-  forAll< execPolicy >( offdiag_nnz, [=] GEOSX_HOST_DEVICE ( HYPRE_Int const i )
+  forAll< hypre::execPolicy >( offdiag_nnz, [=] GEOSX_HYPRE_HOST_DEVICE ( HYPRE_Int const i )
   {
     ptr_offdiag_data[i] *= scalingFactor;
   } );
@@ -819,7 +819,7 @@ void HypreMatrix::leftScale( HypreVector const & vec )
   HYPRE_Real * const va_diag = hypre_CSRMatrixData( csr_diag );
   HYPRE_Real * const va_offdiag = hypre_CSRMatrixData( csr_offdiag );
 
-  forAll< execPolicy >( numLocalRows(), [=] GEOSX_HOST_DEVICE ( localIndex const i )
+  forAll< hypre::execPolicy >( numLocalRows(), [=] GEOSX_HYPRE_HOST_DEVICE ( localIndex const i )
   {
     for( HYPRE_Int j = ia_diag[i]; j < ia_diag[i + 1]; ++j )
     {
@@ -971,7 +971,7 @@ void HypreMatrix::extractDiagonal( HypreVector & dst ) const
   HYPRE_Int const * const ja        = hypre_CSRMatrixJ( csr );
   HYPRE_Real const * const va       = hypre_CSRMatrixData( csr );
 
-  forAll< execPolicy >( numLocalRows(), [=] GEOSX_HOST_DEVICE ( localIndex const localRow )
+  forAll< hypre::execPolicy >( numLocalRows(), [=] GEOSX_HYPRE_HOST_DEVICE ( localIndex const localRow )
   {
     for( HYPRE_Int j = ia[localRow]; j < ia[localRow + 1]; ++j )
     {
@@ -1173,7 +1173,7 @@ void HypreMatrix::print( std::ostream & os ) const
         {
 
           sprintf( str,
-#if defined(GEOSX_USE_CUDA) && defined(GEOSX_LA_INTERFACE_HYPRE)
+#ifdef GEOSX_USE_HYPRE_CUDA
                    "%i%20i%20i%24.10e\n",
 #else
                    "%i%20lli%20lli%24.10e\n",
@@ -1187,7 +1187,7 @@ void HypreMatrix::print( std::ostream & os ) const
         for( HYPRE_Int j = offdiag_IA[i]; j < offdiag_IA[i + 1]; ++j )
         {
           sprintf( str,
-#if defined(GEOSX_USE_CUDA) && defined(GEOSX_LA_INTERFACE_HYPRE)
+#ifdef GEOSX_USE_HYPRE_CUDA
                    "%i%20i%20i%24.10e\n",
 #else
                    "%i%20lli%20lli%24.10e\n",
@@ -1300,8 +1300,8 @@ real64 HypreMatrix::normInf() const
   HYPRE_Real const * const va_diag = hypre_CSRMatrixData( csr_diag );
   HYPRE_Real const * const va_offdiag = hypre_CSRMatrixData( csr_offdiag );
 
-  RAJA::ReduceMax< ReducePolicy< execPolicy >, HYPRE_Real > maxRowAbsSum( 0.0 );
-  forAll< execPolicy >( numLocalRows(), [=] GEOSX_HOST_DEVICE ( localIndex const i )
+  RAJA::ReduceMax< ReducePolicy< hypre::execPolicy >, HYPRE_Real > maxRowAbsSum( 0.0 );
+  forAll< hypre::execPolicy >( numLocalRows(), [=] GEOSX_HYPRE_HOST_DEVICE ( localIndex const i )
   {
     HYPRE_Real rowAbsSum = 0.0;
     for( HYPRE_Int j = ia_diag[i]; j < ia_diag[i + 1]; ++j )