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

@@ -262,6 +262,9 @@ class MultiFluidBase : public ConstitutiveBase
     GEOS_HOST_DEVICE arrayView2d< real64 const, multifluid::USD_FLUID > totalDensity() const
     { return m_totalDensity.value; }
 
+    GEOS_HOST_DEVICE arrayView3d< real64 const, multifluid::USD_FLUID_DC > dTotalDensity() const
+    { return m_totalDensity.derivs; }
+
     GEOS_HOST_DEVICE arrayView3d< real64 const, multifluid::USD_PHASE > phaseEnthalpy() const
     { return m_phaseEnthalpy.value; }
 

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

@@ -56,6 +56,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 ).
@@ -75,19 +84,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;
+
+  // 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;
 
-  m_table.resize( m_numSteps+1, 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
 
@@ -131,8 +164,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
 
@@ -148,6 +180,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
@@ -160,7 +194,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.
 
@@ -186,22 +220,39 @@ 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 );
+  }
+  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 )
+  {
+    auto const phaseNames = getFluid().phaseNames();
+    for( integer ip = 0; ip < m_numPhases; ++ip )
+    {
+      fprintf( fp, "# columns %d-%d = %s phase fractions\n", columnIndex+1, columnIndex + m_numComponents,
+               phaseNames[ip].c_str() );
+      columnIndex += m_numComponents;
+    }
+  }
 
   for( integer n=0; n<m_table.size( 0 ); ++n )
   {
@@ -224,8 +275,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,7 +306,7 @@ void PVTDriver::compareWithBaseline()
 
       real64 const error = fabs( m_table[row][col]-value ) / ( fabs( value )+1 );
       GEOS_THROW_IF( error > m_baselineTol, "Results do not match baseline at data row " << row+1
-                                                                                         << " (row " << row+m_numColumns << " with header)"
+                                                                                         << " (row " << row+headerRows << " with header)"
                                                                                          << " and column " << col+1, std::runtime_error );
     }
   }
@@ -265,6 +326,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

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

@@ -21,6 +21,8 @@
 
 #include "PVTDriver.hpp"
 
+#include "constitutive/fluid/multifluid/Layouts.hpp"
+
 namespace geos
 {
 
@@ -46,10 +48,17 @@ void PVTDriver::runTest( FLUID_TYPE & fluid, arrayView2d< real64 > const & table
   GEOS_ASSERT_EQ( numComponents, m_feed.size() );
   array2d< real64, compflow::LAYOUT_COMP > const compositionValues( 1, numComponents );
 
+  // extract the molecular weights
+  array1d< real64 > const componentMolarWeights( numComponents );
+  for( integer i = 0; i < numComponents; ++i )
+  {
+    componentMolarWeights[i] = fluid.componentMolarWeights()[i];
+  }
+
   real64 sum = 0.0;
   for( integer i = 0; i < numComponents; ++i )
   {
-    compositionValues[0][i] = m_feed[i] * fluid.componentMolarWeights()[i];
+    compositionValues[0][i] = m_feed[i] * componentMolarWeights[i];
     sum += compositionValues[0][i];
   }
   for( integer i = 0; i < numComponents; ++i )
@@ -65,18 +74,60 @@ void PVTDriver::runTest( FLUID_TYPE & fluid, arrayView2d< real64 > const & table
   integer const numSteps = m_numSteps;
   using ExecPolicy = typename FLUID_TYPE::exec_policy;
   forAll< ExecPolicy >( composition.size( 0 ),
-                        [numPhases, numSteps, kernelWrapper, table, composition] GEOS_HOST_DEVICE ( localIndex const i )
+                        [numPhases, numComponents, numSteps, kernelWrapper, table, composition,
+                         molarWeights=componentMolarWeights.toViewConst(),
+                         outputCompressibility=m_outputCompressibility,
+                         outputPhaseComposition=m_outputPhaseComposition]
+                        GEOS_HOST_DEVICE ( localIndex const i )
   {
+    constexpr real64 epsilon = LvArray::NumericLimits< real64 >::epsilon;
+
+    // Pressure derivative index
+    constexpr integer dPIndex = constitutive::multifluid::DerivativeOffset::dP;
+
+    // Index for start of phase properties
+    integer const PHASE = outputCompressibility != 0 ? TEMP + 3 : TEMP + 2;
+
     for( integer n = 0; n <= numSteps; ++n )
     {
       kernelWrapper.update( i, 0, table( n, PRES ), table( n, TEMP ), composition[i] );
-      table( n, TEMP + 1 ) = kernelWrapper.totalDensity()( i, 0 );
+      real64 const totalDensity = kernelWrapper.totalDensity()( i, 0 );
+      table( n, TEMP + 1 ) = totalDensity;
+
+      if( outputCompressibility != 0 )
+      {
+        real64 compressibility = kernelWrapper.dTotalDensity()( i, 0, dPIndex ) / totalDensity;
+        table( n, TEMP + 2 ) =  LvArray::math::max( 0.0, compressibility );
+      }
 
       for( integer p = 0; p < numPhases; ++p )
       {
-        table( n, TEMP + 2 + p ) = kernelWrapper.phaseFraction()( i, 0, p );
-        table( n, TEMP + 2 + p + numPhases ) = kernelWrapper.phaseDensity()( i, 0, p );
-        table( n, TEMP + 2 + p + 2 * numPhases ) = kernelWrapper.phaseViscosity()( i, 0, 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 );
+      }
+      if( outputPhaseComposition != 0 )
+      {
+        for( integer p = 0; p < numPhases; ++p )
+        {
+          integer const compStartIndex = PHASE + 3 * numPhases + p * numComponents;
+
+          // Convert mass fractions to molar fractions
+          real64 molarSum = 0.0;
+          for( integer ic = 0; ic < numComponents; ++ic )
+          {
+            real64 const massFraction = kernelWrapper.phaseCompFraction()( i, 0, p, ic );
+            table( n, compStartIndex + ic ) = massFraction / molarWeights[ic];
+            molarSum += table( n, compStartIndex + ic );
+          }
+
+          real64 const invMolarSum = epsilon < molarSum ? 1.0/molarSum : 0.0;
+
+          for( integer ic = 0; ic < numComponents; ++ic )
+          {
+            table( n, compStartIndex + ic ) *= invMolarSum;
+          }
+        }
       }
     }
   } );
@@ -85,5 +136,4 @@ void PVTDriver::runTest( FLUID_TYPE & fluid, arrayView2d< real64 > const & table
 
 }
 
-
 #endif /* GEOS_CONSTITUTIVE_PVTDRIVERRUNTEST_HPP_ */

src/coreComponents/schema/docs/PVTDriver.rst

@@ -1,17 +1,19 @@
 
 
-================== ============ ======== ============================================= 
-Name               Type         Default  Description                                   
-================== ============ ======== ============================================= 
-baseline           path         none     Baseline file                                 
-feedComposition    real64_array required Feed composition array [mol fraction]         
-fluid              string       required Fluid to test                                 
-logLevel           integer      0        Log level                                     
-name               string       required A name is required for any non-unique nodes   
-output             string       none     Output file                                   
-pressureControl    string       required Function controlling pressure time history    
-steps              integer      required Number of load steps to take                  
-temperatureControl string       required Function controlling temperature time history 
-================== ============ ======== ============================================= 
+====================== ============ ======== ================================================================ 
+Name                   Type         Default  Description                                                      
+====================== ============ ======== ================================================================ 
+baseline               path         none     Baseline file                                                    
+feedComposition        real64_array required Feed composition array [mol fraction]                            
+fluid                  string       required Fluid to test                                                    
+logLevel               integer      0        Log level                                                        
+name                   string       required A name is required for any non-unique nodes                      
+output                 string       none     Output file                                                      
+outputCompressibility  integer      0        Flag to indicate that the total compressibility should be output 
+outputPhaseComposition integer      0        Flag to indicate that phase compositions should be output        
+pressureControl        string       required Function controlling pressure time history                       
+steps                  integer      required Number of load steps to take                                     
+temperatureControl     string       required Function controlling temperature time history                    
+====================== ============ ======== ================================================================