Implement CO2Brine model on GPU

src/coreComponents/constitutive/CMakeLists.txt

@@ -15,14 +15,14 @@ set( constitutive_headers
      contact/Coulomb.hpp
      fluid/DeadOilFluid.hpp
      fluid/MultiPhaseMultiComponentFluid.hpp     
-     fluid/PVTFunctions/BrineCO2DensityFunction.hpp
-     fluid/PVTFunctions/BrineViscosityFunction.hpp
-     fluid/PVTFunctions/CO2SolubilityFunction.hpp
-     fluid/PVTFunctions/FenghourCO2ViscosityFunction.hpp
+     fluid/PVTFunctions/BrineCO2Density.hpp
+     fluid/PVTFunctions/BrineViscosity.hpp
+     fluid/PVTFunctions/CO2Solubility.hpp
+     fluid/PVTFunctions/FenghourCO2Viscosity.hpp
      fluid/PVTFunctions/FlashModelBase.hpp
      fluid/PVTFunctions/PVTFunctionBase.hpp
-     fluid/PVTFunctions/SpanWagnerCO2DensityFunction.hpp
-     fluid/PVTFunctions/UtilityFunctions.hpp
+     fluid/PVTFunctions/SpanWagnerCO2Density.hpp
+     fluid/PVTFunctions/PVTFunctionHelpers.hpp
      fluid/SingleFluidBase.hpp
      fluid/singleFluidSelector.hpp
      fluid/MultiFluidBase.hpp
@@ -75,13 +75,12 @@ set( constitutive_sources
      contact/Coulomb.cpp
      fluid/DeadOilFluid.cpp
      fluid/CompressibleSinglePhaseFluid.cpp
-     fluid/MultiPhaseMultiComponentFluid.cpp     
-     fluid/PVTFunctions/BrineCO2DensityFunction.cpp
-     fluid/PVTFunctions/BrineViscosityFunction.cpp 
-     fluid/PVTFunctions/CO2SolubilityFunction.cpp
-     fluid/PVTFunctions/FenghourCO2ViscosityFunction.cpp
-     fluid/PVTFunctions/SpanWagnerCO2DensityFunction.cpp    
-     fluid/PVTFunctions/UtilityFunctions.cpp
+     fluid/MultiPhaseMultiComponentFluid.cpp
+     fluid/PVTFunctions/BrineCO2Density.cpp
+     fluid/PVTFunctions/BrineViscosity.cpp
+     fluid/PVTFunctions/CO2Solubility.cpp
+     fluid/PVTFunctions/FenghourCO2Viscosity.cpp
+     fluid/PVTFunctions/SpanWagnerCO2Density.cpp
      fluid/SingleFluidBase.cpp
      fluid/MultiFluidBase.cpp
      fluid/CompressibleSinglePhaseFluid.cpp

src/coreComponents/constitutive/docs/CO2Brine.rst

@@ -189,11 +189,11 @@ where :math:`m` is the user-defined salinity.
 Parameters
 =========================
 
-The model is represented by ``<MultiPhaseMultiComponentFluid>`` node in the input.
+The model is represented by ``<CO2BrineFluid>`` node in the input.
 
 The following attributes are supported:
 
-.. include:: ../../../coreComponents/fileIO/schema/docs/MultiPhaseMultiComponentFluid.rst
+.. include:: ../../../coreComponents/fileIO/schema/docs/CO2BrineFluid.rst
 
 Supported phase names are:
 
@@ -219,7 +219,7 @@ Example
 .. code-block:: xml
 
     <Constitutive>
-        <MultiPhaseMultiComponentFluid
+        <CO2BrineFluid
           name="fluid"
           phaseNames="{ gas, water }"
           componentNames="{ co2, water }"

src/coreComponents/constitutive/fluid/BlackOilFluid.hpp

@@ -32,6 +32,8 @@ class BlackOilFluid : public MultiFluidPVTPackageWrapper
 {
 public:
 
+  using exec_policy = serialPolicy;
+
   BlackOilFluid( string const & name, Group * const parent );
 
   virtual ~BlackOilFluid() override;

src/coreComponents/constitutive/fluid/CompositionalMultiphaseFluid.hpp

@@ -30,6 +30,8 @@ class CompositionalMultiphaseFluid : public MultiFluidPVTPackageWrapper
 {
 public:
 
+  using exec_policy = serialPolicy;
+
   CompositionalMultiphaseFluid( string const & name, Group * const parent );
 
   virtual ~CompositionalMultiphaseFluid() override;

src/coreComponents/constitutive/fluid/DeadOilFluid.hpp

@@ -30,7 +30,7 @@ namespace constitutive
 
 /**
  * @brief Kernel wrapper class for DeadOilFluild
- *        This kernel can be called on the GPU (after adding GEOSX_HOST_DEVICE flags in this class and in MultiFluidBase)
+ *        This kernel can be called on the GPU
  */
 class DeadOilFluidUpdate final : public MultiFluidBaseUpdate
 {
@@ -122,7 +122,7 @@ class DeadOilFluidUpdate final : public MultiFluidBaseUpdate
   /// Deleted move assignment operator
   DeadOilFluidUpdate & operator=( DeadOilFluidUpdate && ) = delete;
 
-  //GEOSX_HOST_DEVICE
+  GEOSX_HOST_DEVICE
   GEOSX_FORCE_INLINE
   virtual void compute( real64 const pressure,
                         real64 const temperature,
@@ -134,7 +134,7 @@ class DeadOilFluidUpdate final : public MultiFluidBaseUpdate
                         arraySlice2d< real64 > const & phaseCompFraction,
                         real64 & totalDensity ) const override;
 
-  //GEOSX_HOST_DEVICE
+  GEOSX_HOST_DEVICE
   GEOSX_FORCE_INLINE
   virtual void compute( real64 const pressure,
                         real64 const temperature,
@@ -164,7 +164,7 @@ class DeadOilFluidUpdate final : public MultiFluidBaseUpdate
                         real64 & dTotalDensity_dTemperature,
                         arraySlice1d< real64 > const & dTotalDensity_dGlobalCompFraction ) const override;
 
-  //GEOSX_HOST_DEVICE
+  GEOSX_HOST_DEVICE
   GEOSX_FORCE_INLINE
   virtual void update( localIndex const k,
                        localIndex const q,
@@ -203,24 +203,24 @@ class DeadOilFluidUpdate final : public MultiFluidBaseUpdate
 
 private:
 
-  //GEOSX_HOST_DEVICE
+  GEOSX_HOST_DEVICE
   GEOSX_FORCE_INLINE
   void computeDensities( real64 pressure,
                          arraySlice1d< real64 > const & phaseMassDens ) const;
 
-  //GEOSX_HOST_DEVICE
+  GEOSX_HOST_DEVICE
   GEOSX_FORCE_INLINE
   void computeDensities( real64 pressure,
                          arraySlice1d< real64 > const & phaseMassDens,
                          arraySlice1d< real64 > const & dPhaseMassDens_dPres,
                          arraySlice2d< real64 > const & dPhaseMassDens_dGlobalCompFrac ) const;
 
-  //GEOSX_HOST_DEVICE
+  GEOSX_HOST_DEVICE
   GEOSX_FORCE_INLINE
   void computeViscosities( real64 pressure,
                            arraySlice1d< real64 > const & phaseVisc ) const;
 
-  //GEOSX_HOST_DEVICE
+  GEOSX_HOST_DEVICE
   GEOSX_FORCE_INLINE
   void computeViscosities( real64 pressure,
                            arraySlice1d< real64 > const & phaseVisc,
@@ -260,6 +260,8 @@ class DeadOilFluid : public MultiFluidBase
 {
 public:
 
+  using exec_policy = parallelDevicePolicy<>;
+
   struct PhaseType
   {
     static constexpr integer OIL            = 0;
@@ -432,7 +434,7 @@ class DeadOilFluid : public MultiFluidBase
 
 };
 
-//GEOSX_HOST_DEVICE
+GEOSX_HOST_DEVICE
 GEOSX_FORCE_INLINE
 void DeadOilFluidUpdate::computeDensities( real64 pressure,
                                            arraySlice1d< real64 > const & phaseMassDens ) const
@@ -467,7 +469,7 @@ void DeadOilFluidUpdate::computeDensities( real64 pressure,
   }
 }
 
-//GEOSX_HOST_DEVICE
+GEOSX_HOST_DEVICE
 GEOSX_FORCE_INLINE
 void DeadOilFluidUpdate::computeDensities( real64 pressure,
                                            arraySlice1d< real64 > const & phaseMassDens,
@@ -516,7 +518,7 @@ void DeadOilFluidUpdate::computeDensities( real64 pressure,
   }
 }
 
-//GEOSX_HOST_DEVICE
+GEOSX_HOST_DEVICE
 GEOSX_FORCE_INLINE
 void DeadOilFluidUpdate::computeViscosities( real64 pressure,
                                              arraySlice1d< real64, 0 > const & phaseVisc ) const
@@ -547,7 +549,7 @@ void DeadOilFluidUpdate::computeViscosities( real64 pressure,
   }
 }
 
-//GEOSX_HOST_DEVICE
+GEOSX_HOST_DEVICE
 GEOSX_FORCE_INLINE
 void DeadOilFluidUpdate::computeViscosities( real64 pressure,
                                              arraySlice1d< real64 > const & phaseVisc,
@@ -583,7 +585,7 @@ void DeadOilFluidUpdate::computeViscosities( real64 pressure,
   }
 }
 
-//GEOSX_HOST_DEVICE
+GEOSX_HOST_DEVICE
 GEOSX_FORCE_INLINE
 void DeadOilFluidUpdate::compute( real64 pressure,
                                   real64 temperature,
@@ -632,7 +634,7 @@ void DeadOilFluidUpdate::compute( real64 pressure,
   totalDens = 1.0 / totalDens;
 }
 
-//GEOSX_HOST_DEVICE
+GEOSX_HOST_DEVICE
 GEOSX_FORCE_INLINE
 void DeadOilFluidUpdate::compute( real64 pressure,
                                   real64 temperature,
@@ -660,7 +662,7 @@ void DeadOilFluidUpdate::compute( real64 pressure,
                                   real64 & totalDensity,
                                   real64 & dTotalDensity_dPressure,
                                   real64 & dTotalDensity_dTemperature,
-                                  arraySlice1d< real64, 0 > const & dTotalDensity_dGlobalCompFraction ) const
+                                  arraySlice1d< real64 > const & dTotalDensity_dGlobalCompFraction ) const
 {
   GEOSX_UNUSED_VAR( temperature );
   GEOSX_UNUSED_VAR( dPhaseFraction_dTemperature );
@@ -702,7 +704,7 @@ void DeadOilFluidUpdate::compute( real64 pressure,
     dPhaseFraction_dPressure[ip] = 0.0;
     for( localIndex ic = 0; ic < nComps; ++ic )
     {
-      dPhaseFraction_dGlobalCompFraction[ip][ic] = (ip == ic) ? 1-composition[ip] : -composition[ip];
+      dPhaseFraction_dGlobalCompFraction[ip][ic] = (ip == ic) ? 1.0 : 0.0;
 
       phaseCompFraction[ip][ic] = (ip == ic) ? 1.0 : 0.0;
       dPhaseCompFraction_dPressure[ip][ic] = 0.0;

src/coreComponents/constitutive/fluid/MultiFluidBase.hpp

@@ -163,6 +163,7 @@ class MultiFluidBaseUpdate
 
 private:
 
+  GEOSX_HOST_DEVICE
   virtual void compute( real64 const pressure,
                         real64 const temperature,
                         arraySlice1d< real64 const > const & composition,
@@ -173,6 +174,7 @@ class MultiFluidBaseUpdate
                         arraySlice2d< real64 > const & phaseCompFraction,
                         real64 & totalDensity ) const = 0;
 
+  GEOSX_HOST_DEVICE
   virtual void compute( real64 const pressure,
                         real64 const temperature,
                         arraySlice1d< real64 const > const & composition,
@@ -201,12 +203,12 @@ class MultiFluidBaseUpdate
                         real64 & dTotalDensity_dTemperature,
                         arraySlice1d< real64 > const & dTotalDensity_dGlobalCompFraction ) const = 0;
 
+  GEOSX_HOST_DEVICE
   virtual void update( localIndex const k,
                        localIndex const q,
                        real64 const pressure,
                        real64 const temperature,
                        arraySlice1d< real64 const > const & composition ) const = 0;
-
 };
 
 class MultiFluidBase : public ConstitutiveBase

src/coreComponents/constitutive/fluid/MultiFluidPVTPackageWrapper.cpp

@@ -81,6 +81,7 @@ MultiFluidPVTPackageWrapper::deliverClone( string const & name,
   return clone;
 }
 
+GEOSX_HOST_DEVICE
 void MultiFluidPVTPackageWrapperUpdate::compute( real64 pressure,
                                                  real64 temperature,
                                                  arraySlice1d< real64 const > const & composition,
@@ -91,6 +92,9 @@ void MultiFluidPVTPackageWrapperUpdate::compute( real64 pressure,
                                                  arraySlice2d< real64 > const & phaseCompFrac,
                                                  real64 & totalDens ) const
 {
+#if defined(__CUDA_ARCH__)
+  GEOSX_ERROR( "This function cannot be used on GPU" );
+#else
   localIndex const NC = m_componentMolarWeight.size();
   localIndex const NP = m_phaseTypes.size();
 
@@ -195,8 +199,10 @@ void MultiFluidPVTPackageWrapperUpdate::compute( real64 pressure,
     // 5.2. Invert the previous quantity to get actual density
     totalDens = 1.0 / totalDens;
   }
+#endif
 }
 
+GEOSX_HOST_DEVICE
 void MultiFluidPVTPackageWrapperUpdate::compute( real64 pressure,
                                                  real64 temperature,
                                                  arraySlice1d< real64 const > const & composition,
@@ -225,6 +231,9 @@ void MultiFluidPVTPackageWrapperUpdate::compute( real64 pressure,
                                                  real64 & dTotalDensity_dTemperature,
                                                  arraySlice1d< real64 > const & dTotalDensity_dGlobalCompFraction ) const
 {
+#if defined(__CUDA_ARCH__)
+  GEOSX_ERROR( "This function cannot be used on GPU" );
+#else
 // 0. make shortcut structs to avoid long names (TODO maybe remove)
   CompositionalVarContainer< 1 > phaseFrac{
     phaseFraction,
@@ -268,11 +277,6 @@ void MultiFluidPVTPackageWrapperUpdate::compute( real64 pressure,
     dTotalDensity_dGlobalCompFraction
   };
 
-#if defined(__CUDACC__)
-  // For some reason nvcc thinks these aren't used.
-  GEOSX_UNUSED_VAR( phaseFrac, phaseDens, phaseMassDens, phaseVisc, phaseCompFrac, totalDens );
-#endif
-
   localIndex constexpr maxNumComp = MultiFluidBase::MAX_NUM_COMPONENTS;
   localIndex const NC = numComponents();
   localIndex const NP = numPhases();
@@ -489,6 +493,7 @@ void MultiFluidPVTPackageWrapperUpdate::compute( real64 pressure,
       totalDens.dComp[jc] *= minusDens2;
     }
   }
+#endif
 }
 
 } //namespace constitutive

src/coreComponents/constitutive/fluid/MultiFluidPVTPackageWrapper.hpp

@@ -117,6 +117,7 @@ class MultiFluidPVTPackageWrapperUpdate final : public MultiFluidBaseUpdate
   /// Deleted move assignment operator
   MultiFluidPVTPackageWrapperUpdate & operator=( MultiFluidPVTPackageWrapperUpdate && ) = delete;
 
+  GEOSX_HOST_DEVICE
   virtual void compute( real64 const pressure,
                         real64 const temperature,
                         arraySlice1d< real64 const > const & composition,
@@ -127,6 +128,7 @@ class MultiFluidPVTPackageWrapperUpdate final : public MultiFluidBaseUpdate
                         arraySlice2d< real64 > const & phaseCompFraction,
                         real64 & totalDensity ) const override;
 
+  GEOSX_HOST_DEVICE
   virtual void compute( real64 const pressure,
                         real64 const temperature,
                         arraySlice1d< real64 const > const & composition,
@@ -155,6 +157,7 @@ class MultiFluidPVTPackageWrapperUpdate final : public MultiFluidBaseUpdate
                         real64 & dTotalDensity_dTemperature,
                         arraySlice1d< real64 > const & dTotalDensity_dGlobalCompFraction ) const override;
 
+  GEOSX_HOST_DEVICE
   GEOSX_FORCE_INLINE
   virtual void update( localIndex const k,
                        localIndex const q,