Add phase composition to PVTDriver output

src/coreComponents/constitutive/docs/PVTDriver.rst

@@ -102,9 +102,12 @@ The file is a simple ASCII format with a brief header followed by test data:
   4.0816e-02 3.0000e+06 3.5000e+02 4.9901e+01 1.0000e+00 4.1563e-11 4.9901e+01 1.0066e+03 1.7778e-05 9.9525e-04
   ...
 
-Note that the number of columns will depend on how may phases are present. 
+Note that the number of columns will depend on how many phases and components are present.
 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.
+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/MultiFluidBase.hpp

@@ -463,6 +463,46 @@ class MultiFluidBase : public ConstitutiveBase
     PhaseComp::ViewType m_phaseCompFraction;
     FluidProp::ViewType m_totalDensity;
 
+public:
+    /**
+     * @brief Calculate the total fluid compressibility
+     * @param i Element index
+     * @param q Quadrature node index
+     * @return The total fluid compressibility
+     */
+    GEOS_HOST_DEVICE
+    real64 totalCompressibility( integer const i, integer const q ) const
+    {
+      real64 const totalFluidDensity = totalDensity()( i, q );
+      real64 const dTotalFluidDensity_dP = m_totalDensity.derivs( i, q, multifluid::DerivativeOffset::dP );
+      return 0.0 < totalFluidDensity ? dTotalFluidDensity_dP / totalFluidDensity : 0.0;
+    }
+
+    /**
+     * @brief Extract the phase mole fractions for a phase
+     * @param i Element index
+     * @param q Quadrature node index
+     * @param p Phase index
+     * @param[out] moleFractions The calculated mole fractions
+     */
+    template< typename OUT_ARRAY >
+    GEOS_HOST_DEVICE
+    void phaseCompMoleFraction( integer const i,
+                                integer const q,
+                                integer const p,
+                                OUT_ARRAY && moleFractions ) const
+    {
+      integer const numComponents = m_componentMolarWeight.size( 0 );
+      for( integer ic = 0; ic < numComponents; ++ic )
+      {
+        moleFractions[ic] = m_phaseCompFraction.value( i, q, p, ic );
+      }
+      detail::convertToMoleFractions( numComponents,
+                                      m_componentMolarWeight,
+                                      moleFractions,
+                                      moleFractions );
+    }
+
 private:
 
     /**

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

@@ -130,6 +130,42 @@ struct MultiFluidVar
   }
 };
 
+namespace detail
+{
+/**
+ * @brief Utility function to convert mass fractions to mole fractions
+ * @tparam ARRAY1 the type of array storing the component molar weights
+ * @tparam ARRAY2 the type of array storing the component mass fractions
+ * @tparam ARRAY3 the type of array storing the component mole fractions
+ * @param[in] componentMolarWeight the component molar weights
+ * @param[in] massFractions the component mass fractions
+ * @param[out] moleFractions the newly converted component mole fractions
+ */
+template< typename ARRAY1, typename ARRAY2, typename ARRAY3 >
+GEOS_HOST_DEVICE
+void convertToMoleFractions( integer numComponents,
+                             ARRAY1 const & componentMolarWeight,
+                             ARRAY2 const & massFractions,
+                             ARRAY3 & moleFractions )
+{
+  real64 totalMolality = 0.0;
+  for( integer ic = 0; ic < numComponents; ++ic )
+  {
+    moleFractions[ic] = massFractions[ic] / componentMolarWeight[ic];
+    totalMolality += moleFractions[ic];
+  }
+
+  constexpr real64 epsilon = LvArray::NumericLimits< real64 >::epsilon;
+
+  real64 const totalMolalityInv = epsilon < totalMolality ? 1.0 / totalMolality : 0.0;
+  for( integer ic = 0; ic < numComponents; ++ic )
+  {
+    moleFractions[ic] *= totalMolalityInv;
+  }
+}
+
+} // namespace detail
+
 } // namespace constitutive
 
 } // namespace geos

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

@@ -26,6 +26,7 @@
 #include "fileIO/Outputs/OutputBase.hpp"
 #include "functions/FunctionManager.hpp"
 #include "functions/TableFunction.hpp"
+#include "codingUtilities/StringUtilities.hpp"
 
 #include <fstream>
 
@@ -57,6 +58,15 @@ PVTDriver::PVTDriver( const string & name,
     setInputFlag( InputFlags::REQUIRED ).
     setDescription( "Function controlling temperature time history" );
 
+  registerWrapper( viewKeyStruct::outputCompressibilityString(), &m_outputCompressibility ).
+    setInputFlag( InputFlags::OPTIONAL ).
+    setApplyDefaultValue( 0 ).
+    setDescription( "Flag to indicate that the total compressibility should be output" );
+
+  registerWrapper( viewKeyStruct::outputPhaseCompositionString(), &m_outputPhaseComposition ).
+    setInputFlag( InputFlags::OPTIONAL ).
+    setApplyDefaultValue( 0 ).
+    setDescription( "Flag to indicate that phase compositions should be output" );
 
   //todo refactor in mother class
   registerWrapper( viewKeyStruct::numStepsString(), &m_numSteps ).
@@ -76,19 +86,43 @@ PVTDriver::PVTDriver( const string & name,
 
 void PVTDriver::postProcessInput()
 {
+  // Validate some inputs
+  GEOS_ERROR_IF( m_outputCompressibility != 0 && m_outputCompressibility != 1,
+                 getWrapperDataContext( viewKeyStruct::outputCompressibilityString() ) <<
+                 ": option can be either 0 (false) or 1 (true)" );
+
+  GEOS_ERROR_IF( m_outputPhaseComposition != 0 && m_outputPhaseComposition != 1,
+                 getWrapperDataContext( viewKeyStruct::outputPhaseCompositionString() ) <<
+                 ": option can be either 0 (false) or 1 (true)" );
+
   // get number of phases and components
 
-  ConstitutiveManager & constitutiveManager = this->getGroupByPath< ConstitutiveManager >( "/Problem/domain/Constitutive" );
-  MultiFluidBase & baseFluid = constitutiveManager.getGroup< MultiFluidBase >( m_fluidName );
+  constitutive::MultiFluidBase & baseFluid = getFluid();
 
   m_numPhases = baseFluid.numFluidPhases();
   m_numComponents = baseFluid.numFluidComponents();
 
-  // resize data table to fit number of timesteps and fluid phases:
-  // (numRows,numCols) = (numSteps+1,4+3*numPhases)
-  // column order = time, pressure, temp, totalDensity, phaseFraction_{1:NP}, phaseDensity_{1:NP}, phaseViscosity_{1:NP}
+  // Number of rows in numSteps+1
+  integer const numRows = m_numSteps+1;
 
-  m_table.resize( m_numSteps+1, 3*m_numPhases+4 );
+  // Number of columns depends on options
+  // Default column order = time, pressure, temp, totalDensity, phaseFraction_{1:NP}, phaseDensity_{1:NP}, phaseViscosity_{1:NP}
+  integer numCols = 3*m_numPhases+4;
+
+  // If the total compressibility is requested then add a column
+  if( m_outputCompressibility != 0 )
+  {
+    numCols++;
+  }
+
+  // If phase compositions are required we add {1:NP*NC} phase compositions
+  if( m_outputPhaseComposition != 0 )
+  {
+    numCols += m_numPhases * m_numComponents;
+  }
+
+  // resize data table to fit number of timesteps and fluid phases:
+  m_table.resize( numRows, numCols );
 
   // initialize functions
 
@@ -132,8 +166,7 @@ bool PVTDriver::execute( real64 const GEOS_UNUSED_PARAM( time_n ),
   // get the fluid out of the constitutive manager.
   // for the moment it is of type MultiFluidBase.
 
-  ConstitutiveManager & constitutiveManager = this->getGroupByPath< ConstitutiveManager >( "/Problem/domain/Constitutive" );
-  MultiFluidBase & baseFluid = constitutiveManager.getGroup< MultiFluidBase >( m_fluidName );
+  constitutive::MultiFluidBase & baseFluid = getFluid();
 
   // depending on logLevel, print some useful info
 
@@ -149,6 +182,8 @@ bool PVTDriver::execute( real64 const GEOS_UNUSED_PARAM( time_n ),
     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 Compressibility . " << m_outputCompressibility );
+    GEOS_LOG_RANK_0( "  Output Phase Comp. ..... " << m_outputPhaseComposition );
   }
 
   // create a dummy discretization with one quadrature point for
@@ -161,7 +196,7 @@ bool PVTDriver::execute( real64 const GEOS_UNUSED_PARAM( time_n ),
   discretization.resize( 1 );   // one element
   baseFluid.allocateConstitutiveData( discretization, 1 );   // one quadrature point
 
-  // pass the solid through the ConstitutivePassThru to downcast from the
+  // pass the fluid through the ConstitutivePassThru to downcast from the
   // base type to a known model type.  the lambda here then executes the
   // appropriate test driver. note that these calls will move data to device if available.
 
@@ -187,22 +222,42 @@ bool PVTDriver::execute( real64 const GEOS_UNUSED_PARAM( time_n ),
   return false;
 }
 
-
-
 void PVTDriver::outputResults()
 {
   // TODO: improve file path output to grab command line -o directory
   //       for the moment, we just use the specified m_outputFile directly
 
   FILE * fp = fopen( m_outputFile.c_str(), "w" );
 
-  fprintf( fp, "# column 1 = time\n" );
-  fprintf( fp, "# column 2 = pressure\n" );
-  fprintf( fp, "# column 3 = temperature\n" );
-  fprintf( fp, "# column 4 = density\n" );
-  fprintf( fp, "# columns %d-%d = phase fractions\n", 5, 4+m_numPhases );
-  fprintf( fp, "# columns %d-%d = phase densities\n", 5+m_numPhases, 4+2*m_numPhases );
-  fprintf( fp, "# columns %d-%d = phase viscosities\n", 5+2*m_numPhases, 4+3*m_numPhases );
+  integer columnIndex = 0;
+  fprintf( fp, "# column %d = time\n", ++columnIndex );
+  fprintf( fp, "# column %d = pressure\n", ++columnIndex );
+  fprintf( fp, "# column %d = temperature\n", ++columnIndex );
+  fprintf( fp, "# column %d = density\n", ++columnIndex );
+  if( m_outputCompressibility != 0 )
+  {
+    fprintf( fp, "# column %d = total compressibility\n", ++columnIndex );
+  }
+
+  auto const phaseNames = getFluid().phaseNames();
+
+  fprintf( fp, "# columns %d-%d = phase fractions\n", columnIndex+1, columnIndex + m_numPhases );
+  columnIndex += m_numPhases;
+  fprintf( fp, "# columns %d-%d = phase 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;
+
+  if( m_outputPhaseComposition != 0 )
+  {
+    string const componentNames = stringutilities::join( getFluid().componentNames(), ", " );
+    for( integer ip = 0; ip < m_numPhases; ++ip )
+    {
+      fprintf( fp, "# columns %d-%d = %s phase fractions [%s]\n", columnIndex+1, columnIndex + m_numComponents,
+               phaseNames[ip].c_str(), componentNames.c_str() );
+      columnIndex += m_numComponents;
+    }
+  }
 
   for( integer n=0; n<m_table.size( 0 ); ++n )
   {
@@ -225,8 +280,18 @@ void PVTDriver::compareWithBaseline()
 
   // discard file header
 
+  integer headerRows = 7;
+  if( m_outputCompressibility )
+  {
+    headerRows++;
+  }
+  if( m_outputPhaseComposition )
+  {
+    headerRows += getFluid().numFluidPhases();
+  }
+
   string line;
-  for( integer row=0; row < 7; ++row )
+  for( integer row=0; row < headerRows; ++row )
   {
     getline( file, line );
   }
@@ -245,9 +310,10 @@ void PVTDriver::compareWithBaseline()
       file >> value;
 
       real64 const error = fabs( m_table[row][col]-value ) / ( fabs( value )+1 );
-      GEOS_THROW_IF( error > MultiFluidConstants::baselineTolerance, "Results do not match baseline at data row " << row+1
-                                                                                                                  << " (row " << row+m_numColumns << " with header)"
-                                                                                                                  << " and column " << col+1, std::runtime_error );
+      GEOS_THROW_IF( error > MultiFluidConstants::baselineTolerance,
+                     "Results do not match baseline at data row " << row+1
+                                                                  << " (row " << row+headerRows << " with header)"
+                                                                  << " and column " << col+1, std::runtime_error );
     }
   }
 
@@ -266,6 +332,13 @@ void PVTDriver::compareWithBaseline()
   file.close();
 }
 
+constitutive::MultiFluidBase &
+PVTDriver::getFluid()
+{
+  ConstitutiveManager & constitutiveManager = this->getGroupByPath< ConstitutiveManager >( "/Problem/domain/Constitutive" );
+  MultiFluidBase & baseFluid = constitutiveManager.getGroup< MultiFluidBase >( m_fluidName );
+  return baseFluid;
+}
 
 REGISTER_CATALOG_ENTRY( TaskBase,
                         PVTDriver,

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

@@ -24,6 +24,11 @@
 namespace geos
 {
 
+namespace constitutive
+{
+class MultiFluidBase;
+}
+
 /**
  * @class PVTDriver
  *
@@ -68,6 +73,11 @@ class PVTDriver : public TaskBase
 
 private:
 
+  /**
+   * @brief Get the fluid model from the catalog
+   */
+  constitutive::MultiFluidBase & getFluid();
+
   /**
    * @struct viewKeyStruct holds char strings and viewKeys for fast lookup
    */
@@ -80,17 +90,20 @@ 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 * outputCompressibilityString() { return "outputCompressibility"; }
+    constexpr static char const * outputPhaseCompositionString() { return "outputPhaseComposition"; }
   };
 
   integer m_numSteps;      ///< Number of load steps
-  integer m_numColumns;    ///< Number of columns in data table (depends on number of fluid phases)
   integer m_numPhases;     ///< Number of fluid phases
   integer m_numComponents; ///< Number of fluid components
 
-  string m_fluidName;               ///< Fluid identifier
-  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)
+  string m_fluidName;                   ///< Fluid identifier
+  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_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
 
   array1d< real64 > m_feed;  ///< User specified feed composition
   array2d< real64 > m_table; ///< Table storing time-history of input/output