Link compositional fluid model to static function implementation

src/coreComponents/constitutive/CMakeLists.txt

@@ -62,6 +62,7 @@ set( constitutive_headers
      fluid/multifluid/CO2Brine/functions/SpanWagnerCO2Density.hpp
      fluid/multifluid/CO2Brine/functions/WaterDensity.hpp
      fluid/multifluid/compositional/functions/CompositionalProperties.hpp     
+     fluid/multifluid/compositional/functions/CompositionalPropertiesImpl.hpp     
      fluid/multifluid/compositional/functions/CubicEOSPhaseModel.hpp     
      fluid/multifluid/compositional/functions/KValueInitialization.hpp
      fluid/multifluid/compositional/functions/NegativeTwoPhaseFlash.hpp
@@ -70,6 +71,8 @@ set( constitutive_headers
      fluid/multifluid/compositional/models/CompositionalDensity.hpp
      fluid/multifluid/compositional/models/ConstantViscosity.hpp
      fluid/multifluid/compositional/models/FunctionBase.hpp
+     fluid/multifluid/compositional/models/LohrenzBrayClarkViscosity.hpp
+     fluid/multifluid/compositional/models/LohrenzBrayClarkViscosityImpl.hpp
      fluid/multifluid/compositional/models/NegativeTwoPhaseFlashModel.hpp
      fluid/multifluid/compositional/models/NullModel.hpp
      fluid/multifluid/compositional/models/PhaseModel.hpp
@@ -214,13 +217,16 @@ set( constitutive_sources
      fluid/multifluid/CO2Brine/functions/CO2EOSSolver.cpp
      fluid/multifluid/CO2Brine/functions/PureWaterProperties.cpp
      fluid/multifluid/CO2Brine/functions/WaterDensity.cpp
-     fluid/multifluid/compositional/functions/CompositionalProperties.cpp
      fluid/multifluid/compositional/models/CompositionalDensity.cpp
      fluid/multifluid/compositional/models/ConstantViscosity.cpp
+     fluid/multifluid/compositional/models/LohrenzBrayClarkViscosity.cpp
      fluid/multifluid/compositional/models/NegativeTwoPhaseFlashModel.cpp
      fluid/multifluid/compositional/CompositionalMultiphaseFluid.cpp
      fluid/multifluid/compositional/CompositionalMultiphaseFluidUpdates.cpp
-     fluid/multifluid/compositional/PVTDriverRunTestCompositionalMultiphaseFluid.cpp
+     fluid/multifluid/compositional/PVTDriverRunTestCompositionalPR.cpp
+     fluid/multifluid/compositional/PVTDriverRunTestCompositionalPRLBC.cpp
+     fluid/multifluid/compositional/PVTDriverRunTestCompositionalSRK.cpp
+     fluid/multifluid/compositional/PVTDriverRunTestCompositionalSRKLBC.cpp
      fluid/multifluid/reactive/ReactiveBrineFluid.cpp
      fluid/multifluid/reactive/ReactiveMultiFluid.cpp
      fluid/multifluid/reactive/ReactiveFluidDriver.cpp
@@ -253,10 +259,10 @@ set( constitutive_sources
      relativePermeability/VanGenuchtenStone2RelativePermeability.cpp
      relativePermeability/RelpermDriver.cpp
      relativePermeability/RelpermDriverBrooksCoreyBakerRunTest.cpp
-        relativePermeability/RelpermDriverBrooksCoreyStone2RunTest.cpp
+     relativePermeability/RelpermDriverBrooksCoreyStone2RunTest.cpp
      relativePermeability/RelpermDriverBrooksCoreyRunTest.cpp
      relativePermeability/RelpermDriverVanGenuchtenBakerRunTest.cpp
-        relativePermeability/RelpermDriverVanGenuchtenStone2RunTest.cpp
+     relativePermeability/RelpermDriverVanGenuchtenStone2RunTest.cpp
      relativePermeability/RelpermDriverTableRelativeRunTest.cpp
      relativePermeability/RelpermDriverTableRelativeHysteresisRunTest.cpp
      solid/CompressibleSolid.cpp
@@ -300,7 +306,8 @@ if( ENABLE_PVTPackage )
 
     set( constitutive_sources
          ${constitutive_sources}
-         fluid/multifluid/compositional/CompositionalMultiphaseFluidPVTPackage.cpp )
+         fluid/multifluid/compositional/CompositionalMultiphaseFluidPVTPackage.cpp
+         fluid/multifluid/compositional/PVTDriverRunTestCompositionalMultiphaseFluid.cpp )
 
     add_subdirectory( PVTPackage )
 

src/coreComponents/constitutive/docs/PVTDriver.rst

@@ -107,6 +107,7 @@ In this case, we have a two-phase, two-component mixture.
 The total density is reported in column 4, while phase fractions, phase densities, and phase viscosities are reported in subsequent columns.
 If the ``outputCompressibility`` flag is activated, an extra column will be added for the total fluid compressibility after the density.
 This is defined as :math:`c_t=\frac{1}{\rho_t}\left(\partial{\rho_t}/\partial P\right)` where :math:`\rho_t` is the total density.
+If the ``outputMassDensity`` flag is activated, extra columns will be added for the mass density of each phase.
 The number of columns will also depend on whether the ``outputPhaseComposition`` flag is activated or not. If it is activated, there will be an extra column for the mole fraction of each component in each phase.
 The phase order will match the one defined in the input XML (here, the co2-rich phase followed by the water-rich phase).
 This file can be readily plotted using any number of plotting tools.  Each row corresponds to one timestep of the driver, starting from initial conditions in the first row.

src/coreComponents/constitutive/fluid/multifluid/MultiFluidConstants.hpp

@@ -49,7 +49,7 @@ struct MultiFluidConstants
   /**
    * @brief Max number of SSI iterations
    */
-  static constexpr integer maxSSIIterations = 200;
+  static constexpr integer maxSSIIterations = 1000;
 
   /**
    * @brief Max number of Newton iterations

src/coreComponents/constitutive/fluid/multifluid/MultiFluidSelector.hpp

@@ -43,6 +43,8 @@ void constitutiveUpdatePassThru( MultiFluidBase const & fluid,
                                BlackOilFluid,
                                CompositionalTwoPhasePengRobinsonConstantViscosity,
                                CompositionalTwoPhaseSoaveRedlichKwongConstantViscosity,
+//                               CompositionalTwoPhasePengRobinsonLBCViscosity,
+//                               CompositionalTwoPhaseSoaveRedlichKwongLBCViscosity,
 #ifdef GEOSX_USE_PVTPackage
                                CompositionalMultiphaseFluidPVTPackage,
 #endif
@@ -62,6 +64,8 @@ void constitutiveUpdatePassThru( MultiFluidBase & fluid,
                                BlackOilFluid,
                                CompositionalTwoPhasePengRobinsonConstantViscosity,
                                CompositionalTwoPhaseSoaveRedlichKwongConstantViscosity,
+//                               CompositionalTwoPhasePengRobinsonLBCViscosity,
+//                               CompositionalTwoPhaseSoaveRedlichKwongLBCViscosity,
 #ifdef GEOSX_USE_PVTPackage
                                CompositionalMultiphaseFluidPVTPackage,
 #endif

src/coreComponents/constitutive/fluid/multifluid/MultiFluidUtils.hpp

@@ -63,6 +63,8 @@ struct MultiFluidVarSlice
 
   internal::ArraySliceOrRef< T, DIM, USD > value;        /// variable value
   internal::ArraySliceOrRef< T, DIM + 1, USD_DC > derivs; /// derivative w.r.t. pressure, temperature, compositions
+
+  using ValueType = internal::ArraySliceOrRef< T, DIM, USD >;
 };
 
 /**
@@ -87,13 +89,15 @@ struct MultiFluidVarView
                       ArrayView< T, NDIM + 1, USD_DC > const & derivsSrc ):
     value( valueSrc ),
     derivs( derivsSrc )
-  {};
+  {}
 
   ArrayView< T, NDIM, USD > value;        ///< View into property values
   ArrayView< T, NDIM + 1, USD_DC > derivs; ///< View into property derivatives w.r.t. pressure, temperature, compositions
 
   using SliceType = MultiFluidVarSlice< T, NDIM - 2, USD - 2, USD_DC - 2 >;
 
+  using ValueType = ArrayView< T, NDIM, USD >;
+
   GEOS_HOST_DEVICE
   SliceType operator()( localIndex const k, localIndex const q ) const
   {
@@ -119,6 +123,12 @@ struct MultiFluidVar
   using SliceType = typename ViewType::SliceType;
   using SliceTypeConst = typename ViewTypeConst::SliceType;
 
+  using ValueType = Array< real64, NDIM, PERM >;
+  template< int MAXSIZE >
+  using StackValueType = StackArray< real64, NDIM, MAXSIZE, PERM >;
+  using ViewValueType = typename ViewType::ValueType;
+  using SliceValueType = typename SliceType::ValueType;
+
   ViewType toView()
   {
     return { value.toView(), derivs.toView() };

src/coreComponents/constitutive/fluid/multifluid/PVTDriver.cpp

@@ -61,6 +61,11 @@
     setInputFlag( InputFlags::REQUIRED ).
     setDescription( "Function controlling temperature time history" );
 
+  registerWrapper( viewKeyStruct::outputMassDensityString(), &m_outputMassDensity ).
+    setInputFlag( InputFlags::OPTIONAL ).
+    setApplyDefaultValue( 0 ).
+    setDescription( "Flag to indicate that the mass density of each phase should be output" );
+
   registerWrapper( viewKeyStruct::outputCompressibilityString(), &m_outputCompressibility ).
     setInputFlag( InputFlags::OPTIONAL ).
     setApplyDefaultValue( 0 ).
@@ -90,6 +95,10 @@
 void PVTDriver::postProcessInput()
 {
   // Validate some inputs
+  GEOS_ERROR_IF( m_outputMassDensity != 0 && m_outputMassDensity != 1,
+                 getWrapperDataContext( viewKeyStruct::outputMassDensityString() ) <<
+                 ": option can be either 0 (false) or 1 (true)" );
+
   GEOS_ERROR_IF( m_outputCompressibility != 0 && m_outputCompressibility != 1,
                  getWrapperDataContext( viewKeyStruct::outputCompressibilityString() ) <<
                  ": option can be either 0 (false) or 1 (true)" );
@@ -112,6 +121,12 @@
   // Default column order = time, pressure, temp, totalDensity, phaseFraction_{1:NP}, phaseDensity_{1:NP}, phaseViscosity_{1:NP}
   integer numCols = 3*m_numPhases+4;
 
+  // If the mass density is requested then add NP columns
+  if( m_outputMassDensity != 0 )
+  {
+    numCols += m_numPhases;
+  }
+
   // If the total compressibility is requested then add a column
   if( m_outputCompressibility != 0 )
   {
@@ -185,6 +200,7 @@
     GEOS_LOG_RANK_0( "  Steps .................. " << m_numSteps );
     GEOS_LOG_RANK_0( "  Output ................. " << m_outputFile );
     GEOS_LOG_RANK_0( "  Baseline ............... " << m_baselineFile );
+    GEOS_LOG_RANK_0( "  Output Mass Density .... " << m_outputMassDensity );
     GEOS_LOG_RANK_0( "  Output Compressibility . " << m_outputCompressibility );
     GEOS_LOG_RANK_0( "  Output Phase Comp. ..... " << m_outputPhaseComposition );
   }
@@ -248,6 +264,11 @@
   columnIndex += m_numPhases;
   fprintf( fp, "# columns %d-%d = phase densities\n", columnIndex+1, columnIndex + m_numPhases );
   columnIndex += m_numPhases;
+  if( m_outputMassDensity != 0 )
+  {
+    fprintf( fp, "# columns %d-%d = phase mass densities\n", columnIndex+1, columnIndex + m_numPhases );
+    columnIndex += m_numPhases;
+  }
   fprintf( fp, "# columns %d-%d = phase viscosities\n", columnIndex+1, columnIndex + m_numPhases );
   columnIndex += m_numPhases;
 

src/coreComponents/constitutive/fluid/multifluid/PVTDriver.hpp

@@ -90,6 +90,7 @@ class PVTDriver : public TaskBase
     constexpr static char const * outputString() { return "output"; }
     constexpr static char const * baselineString() { return "baseline"; }
     constexpr static char const * feedString() { return "feedComposition"; }
+    constexpr static char const * outputMassDensityString() { return "outputMassDensity"; }
     constexpr static char const * outputCompressibilityString() { return "outputCompressibility"; }
     constexpr static char const * outputPhaseCompositionString() { return "outputPhaseComposition"; }
   };
@@ -102,6 +103,7 @@ class PVTDriver : public TaskBase
   string m_pressureFunctionName;        ///< Time-dependent function controlling pressure
   string m_temperatureFunctionName;     ///< Time-dependent function controlling temperature
   string m_outputFile;                  ///< Output file (optional, no output if not specified)
+  integer m_outputMassDensity{0};       ///< Flag to indicate that the mass density of each phase should be output
   integer m_outputCompressibility{0};   ///< Flag to indicate that the total compressibility should be output
   integer m_outputPhaseComposition{0};  ///< Flag to indicate that phase compositions should be output
 

src/coreComponents/constitutive/fluid/multifluid/PVTDriverRunTest.hpp

@@ -63,17 +63,23 @@
   // perform fluid update using table (P,T) and save resulting total density, etc.
   // note: column indexing should be kept consistent with output file header below.
 
+  bool const outputMassDensity = (m_outputMassDensity != 0);
+  bool const outputCompressibility = (m_outputCompressibility != 0);
+  bool const outputPhaseComposition = (m_outputPhaseComposition != 0);
+
   integer const numSteps = m_numSteps;
-  integer const outputCompressibility = m_outputCompressibility;
-  integer const outputPhaseComposition = m_outputPhaseComposition;
   using ExecPolicy = typename FLUID_TYPE::exec_policy;
   forAll< ExecPolicy >( composition.size( 0 ),
-                        [ outputCompressibility, outputPhaseComposition, numPhases, numComponents, numSteps, kernelWrapper,
-                          table, composition]
+                        [outputMassDensity, outputCompressibility, outputPhaseComposition,
+                         numPhases, numComponents, numSteps, kernelWrapper,
+                         table, composition]
                         GEOS_HOST_DEVICE ( localIndex const i )
   {
     // Index for start of phase properties
-    integer const PHASE = outputCompressibility != 0 ? TEMP + 3 : TEMP + 2;
+    integer const PHASE_FRACTION = outputCompressibility ? TEMP + 3 : TEMP + 2;
+    integer const PHASE_DENSITY = PHASE_FRACTION + numPhases;
+    integer const PHASE_VISCOSITY = outputMassDensity ? PHASE_DENSITY + 2*numPhases : PHASE_DENSITY + numPhases;
+    integer const PHASE_COMP = PHASE_VISCOSITY + numPhases;
 
     // Temporary space for phase mole fractions
     stackArray1d< real64, constitutive::MultiFluidBase::MAX_NUM_COMPONENTS > phaseComposition( numComponents );
@@ -83,22 +89,29 @@
       kernelWrapper.update( i, 0, table( n, PRES ), table( n, TEMP ), composition[i] );
       table( n, TEMP + 1 ) = kernelWrapper.totalDensity()( i, 0 );
 
-      if( outputCompressibility != 0 )
+      if( outputCompressibility )
       {
         table( n, TEMP + 2 ) = kernelWrapper.totalCompressibility( i, 0 );
       }
 
       for( integer p = 0; p < numPhases; ++p )
       {
-        table( n, PHASE + p ) = kernelWrapper.phaseFraction()( i, 0, p );
-        table( n, PHASE + p + numPhases ) = kernelWrapper.phaseDensity()( i, 0, p );
-        table( n, PHASE + p + 2 * numPhases ) = kernelWrapper.phaseViscosity()( i, 0, p );
+        table( n, PHASE_FRACTION + p ) = kernelWrapper.phaseFraction()( i, 0, p );
+        table( n, PHASE_DENSITY + p ) = kernelWrapper.phaseDensity()( i, 0, p );
+        table( n, PHASE_VISCOSITY + p ) = kernelWrapper.phaseViscosity()( i, 0, p );
+      }
+      if( outputMassDensity )
+      {
+        for( integer p = 0; p < numPhases; ++p )
+        {
+          table( n, PHASE_DENSITY + numPhases + p ) = kernelWrapper.phaseMassDensity()( i, 0, p );
+        }
       }
-      if( outputPhaseComposition != 0 )
+      if( outputPhaseComposition )
       {
         for( integer p = 0; p < numPhases; ++p )
         {
-          integer const compStartIndex = PHASE + 3 * numPhases + p * numComponents;
+          integer const compStartIndex = PHASE_COMP + p * numComponents;
 
           kernelWrapper.phaseCompMoleFraction( i, 0, p, phaseComposition );
           for( integer ic = 0; ic < numComponents; ++ic )