feat: Add temperature correction to density volume shift

src/coreComponents/constitutive/fluid/multifluid/compositional/CompositionalMultiphaseFluid.cpp

@@ -70,10 +70,6 @@ CompositionalMultiphaseFluid( string const & name, Group * const parent )
     setInputFlag( InputFlags::REQUIRED ).
     setDescription( "Component acentric factors" );
 
-  registerWrapper( viewKeyStruct::componentVolumeShiftString(), &m_componentProperties->m_componentVolumeShift ).
-    setInputFlag( InputFlags::OPTIONAL ).
-    setDescription( "Component volume shifts" );
-
   registerWrapper( viewKeyStruct::componentBinaryCoeffString(), &m_componentProperties->m_componentBinaryCoeff ).
     setInputFlag( InputFlags::OPTIONAL ).
     setDescription( "Table of binary interaction coefficients" );
@@ -132,19 +128,11 @@ void CompositionalMultiphaseFluid< FLASH, PHASE1, PHASE2, PHASE3 >::postInputIni
     GEOS_THROW_IF_NE_MSG( array.size(), expected,
                           GEOS_FMT( "{}: invalid number of values in attribute '{}'", getFullName(), attribute ),
                           InputError );
-
   };
   checkInputSize( m_componentProperties->m_componentCriticalPressure, NC, viewKeyStruct::componentCriticalPressureString() );
   checkInputSize( m_componentProperties->m_componentCriticalTemperature, NC, viewKeyStruct::componentCriticalTemperatureString() );
   checkInputSize( m_componentProperties->m_componentAcentricFactor, NC, viewKeyStruct::componentAcentricFactorString() );
 
-  if( m_componentProperties->m_componentVolumeShift.empty() )
-  {
-    m_componentProperties->m_componentVolumeShift.resize( NC );
-    m_componentProperties->m_componentVolumeShift.zero();
-  }
-  checkInputSize( m_componentProperties->m_componentVolumeShift, NC, viewKeyStruct::componentVolumeShiftString() );
-
   array2d< real64 > & componentBinaryCoeff = m_componentProperties->m_componentBinaryCoeff;
   if( componentBinaryCoeff.empty() )
   {

src/coreComponents/constitutive/fluid/multifluid/compositional/CompositionalMultiphaseFluid.hpp

@@ -92,7 +92,6 @@ class CompositionalMultiphaseFluid : public MultiFluidBase
     static constexpr char const * componentCriticalPressureString() { return "componentCriticalPressure"; }
     static constexpr char const * componentCriticalTemperatureString() { return "componentCriticalTemperature"; }
     static constexpr char const * componentAcentricFactorString() { return "componentAcentricFactor"; }
-    static constexpr char const * componentVolumeShiftString() { return "componentVolumeShift"; }
     static constexpr char const * componentBinaryCoeffString() { return "componentBinaryCoeff"; }
   };
 

src/coreComponents/constitutive/fluid/multifluid/compositional/functions/CubicEOSPhaseModel.hpp

@@ -155,10 +155,12 @@ struct CubicEOSPhaseModel
    * @details Computes the dimensional form of the volume shifts given the user defined non-dimensional form.
    * @param[in] numComps The number of components
    * @param[in] componentProperties The compositional model properties
+   * @param[in] volumeShift The input non-dimensional volume shifts
    * @param[out] dimensionalVolumeShift The calculated dimensional volume shifts
    */
   GEOS_FORCE_INLINE
   static void calculateDimensionalVolumeShift( ComponentProperties const & componentProperties,
+                                               arraySlice1d< real64 const > const & volumeShift,
                                                arraySlice1d< real64 > const & dimensionalVolumeShift );
 
   /**
@@ -645,12 +647,13 @@ template< typename EOS_TYPE >
 void
 CubicEOSPhaseModel< EOS_TYPE >::
 calculateDimensionalVolumeShift( ComponentProperties const & componentProperties,
+                                 arraySlice1d< real64 const > const & volumeShift,
                                  arraySlice1d< real64 > const & dimensionalVolumeShift )
 {
   integer const numComps = componentProperties.getNumberOfComponents();
   for( integer ic = 0; ic < numComps; ++ic )
   {
-    real64 const Vs = componentProperties.getComponentVolumeShift()[ic];
+    real64 const Vs = volumeShift[ic];
     real64 const Pc = componentProperties.getComponentCriticalPressure()[ic];
     real64 const Tc = componentProperties.getComponentCriticalTemperature()[ic];
     real64 constexpr omegaB = EOS_TYPE::omegaB;

src/coreComponents/constitutive/fluid/multifluid/compositional/models/ComponentProperties.hpp

@@ -61,21 +61,18 @@ class ComponentProperties final
   arrayView1d< real64 > const & getComponentCriticalPressure() const { return m_componentCriticalPressure; }
   arrayView1d< real64 > const & getComponentCriticalTemperature() const { return m_componentCriticalTemperature; }
   arrayView1d< real64 > const & getComponentAcentricFactor() const { return m_componentAcentricFactor; }
-  arrayView1d< real64 > const & getComponentVolumeShift() const { return m_componentVolumeShift; }
 
   struct KernelWrapper
   {
     KernelWrapper( arrayView1d< real64 const > const & componentMolarWeight,
                    arrayView1d< real64 const > const & componentCriticalPressure,
                    arrayView1d< real64 const > const & componentCriticalTemperature,
                    arrayView1d< real64 const > const & componentAcentricFactor,
-                   arrayView1d< real64 const > const & componentVolumeShift,
                    arrayView2d< real64 const > const & componentBinaryCoeff ):
       m_componentMolarWeight ( componentMolarWeight ),
       m_componentCriticalPressure ( componentCriticalPressure ),
       m_componentCriticalTemperature( componentCriticalTemperature ),
       m_componentAcentricFactor( componentAcentricFactor ),
-      m_componentVolumeShift( componentVolumeShift ),
       m_componentBinaryCoeff( componentBinaryCoeff )
     {}
 
@@ -92,7 +89,6 @@ class ComponentProperties final
       m_componentCriticalPressure.move( space, touch );
       m_componentCriticalTemperature.move( space, touch );
       m_componentAcentricFactor.move( space, touch );
-      m_componentVolumeShift.move( space, touch );
       m_componentBinaryCoeff.move( space, touch );
     }
 
@@ -101,7 +97,6 @@ class ComponentProperties final
     arrayView1d< real64 const > m_componentCriticalPressure;
     arrayView1d< real64 const > m_componentCriticalTemperature;
     arrayView1d< real64 const > m_componentAcentricFactor;
-    arrayView1d< real64 const > m_componentVolumeShift;
     arrayView2d< real64 const > m_componentBinaryCoeff;
   };
 
@@ -115,7 +110,6 @@ class ComponentProperties final
                           m_componentCriticalPressure,
                           m_componentCriticalTemperature,
                           m_componentAcentricFactor,
-                          m_componentVolumeShift,
                           m_componentBinaryCoeff );
   }
 
@@ -126,7 +120,6 @@ class ComponentProperties final
   array1d< real64 > m_componentCriticalPressure;
   array1d< real64 > m_componentCriticalTemperature;
   array1d< real64 > m_componentAcentricFactor;
-  array1d< real64 > m_componentVolumeShift;
   array2d< real64 > m_componentBinaryCoeff;
 };
 

src/coreComponents/constitutive/fluid/multifluid/compositional/models/CompositionalDensity.cpp

@@ -27,6 +27,7 @@ namespace constitutive
 
 namespace compositional
 {
+
 CompositionalDensity::CompositionalDensity( string const & name,
                                             ComponentProperties const & componentProperties,
                                             integer const phaseIndex,
@@ -37,10 +38,13 @@ CompositionalDensity::CompositionalDensity( string const & name,
   string const eosName = equationOfState->m_equationsOfStateNames[phaseIndex];
   m_equationOfState = EnumStrings< EquationOfStateType >::fromString( eosName );
 
+  Parameters const * densityParameters = modelParameters.get< Parameters >();
+
   // Calculate the dimensional volume shift
   m_componentDimensionalVolumeShift.resize( componentProperties.getNumberOfComponents());
   calculateDimensionalVolumeShift( componentProperties,
                                    m_equationOfState,
+                                   densityParameters->m_componentVolumeShift.toSliceConst(),
                                    m_componentDimensionalVolumeShift );
 }
 
@@ -53,25 +57,70 @@ CompositionalDensity::createKernelWrapper() const
 std::unique_ptr< ModelParameters >
 CompositionalDensity::createParameters( std::unique_ptr< ModelParameters > parameters )
 {
-  return EquationOfState::create( std::move( parameters ) );
+  auto params = EquationOfState::create( std::move( parameters ) );
+  return Parameters::create( std::move( params ) );
 }
 
 void CompositionalDensity::calculateDimensionalVolumeShift( ComponentProperties const & componentProperties,
                                                             EquationOfStateType const & equationOfState,
+                                                            arraySlice1d< real64 const > componentVolumeShift,
                                                             arraySlice1d< real64 > componentDimensionalVolumeShift )
 {
   if( equationOfState == EquationOfStateType::PengRobinson )
   {
     CubicEOSPhaseModel< PengRobinsonEOS >::calculateDimensionalVolumeShift( componentProperties,
+                                                                            componentVolumeShift,
                                                                             componentDimensionalVolumeShift );
   }
   else if( equationOfState == EquationOfStateType::SoaveRedlichKwong )
   {
     CubicEOSPhaseModel< SoaveRedlichKwongEOS >::calculateDimensionalVolumeShift( componentProperties,
+                                                                                 componentVolumeShift,
                                                                                  componentDimensionalVolumeShift );
   }
 }
 
+// Compositional density parameters
+
+CompositionalDensity::Parameters::Parameters( std::unique_ptr< ModelParameters > parameters ):
+  ModelParameters( std::move( parameters ) )
+{}
+
+std::unique_ptr< ModelParameters >
+CompositionalDensity::Parameters::create( std::unique_ptr< ModelParameters > parameters )
+{
+  if( parameters && parameters->get< Parameters >() != nullptr )
+  {
+    return parameters;
+  }
+  return std::make_unique< Parameters >( std::move( parameters ) );
+}
+
+void CompositionalDensity::Parameters::registerParametersImpl( MultiFluidBase * fluid )
+{
+  fluid->registerWrapper( viewKeyStruct::componentVolumeShiftString(), &m_componentVolumeShift ).
+    setInputFlag( dataRepository::InputFlags::OPTIONAL ).
+    setDescription( "Component volume shifts" );
+}
+
+void CompositionalDensity::Parameters::postInputInitializationImpl( MultiFluidBase const * fluid,
+                                                                    ComponentProperties const & componentProperties )
+{
+  GEOS_UNUSED_VAR( componentProperties );
+
+  integer const numComponents = fluid->numFluidComponents();
+
+  if( m_componentVolumeShift.empty() )
+  {
+    m_componentVolumeShift.resize( numComponents );
+    m_componentVolumeShift.zero();
+  }
+  GEOS_THROW_IF_NE_MSG( m_componentVolumeShift.size(), numComponents,
+                        GEOS_FMT( "{}: invalid number of values in attribute '{}'", fluid->getFullName(),
+                                  viewKeyStruct::componentVolumeShiftString() ),
+                        InputError );
+}
+
 } // namespace compositional
 
 } // namespace constitutive

src/coreComponents/constitutive/fluid/multifluid/compositional/models/CompositionalDensity.hpp

@@ -102,9 +102,30 @@ class CompositionalDensity : public FunctionBase
   // Create parameters unique to this model
   static std::unique_ptr< ModelParameters > createParameters( std::unique_ptr< ModelParameters > parameters );
 
+  class Parameters : public ModelParameters
+  {
+public:
+    Parameters( std::unique_ptr< ModelParameters > parameters );
+    ~Parameters() override = default;
+
+    static std::unique_ptr< ModelParameters > create( std::unique_ptr< ModelParameters > parameters );
+
+    struct viewKeyStruct
+    {
+      static constexpr char const * componentVolumeShiftString() { return "componentVolumeShift"; }
+    };
+
+    array1d< real64 > m_componentVolumeShift;
+
+protected:
+    void registerParametersImpl( MultiFluidBase * fluid ) override;
+    void postInputInitializationImpl( MultiFluidBase const * fluid, ComponentProperties const & componentProperties ) override;
+  };
+
 private:
   static void calculateDimensionalVolumeShift( ComponentProperties const & componentProperties,
                                                EquationOfStateType const & equationOfState,
+                                               arraySlice1d< real64 const > componentVolumeShift,
                                                arraySlice1d< real64 > componentDimensionalVolumeShift );
 
 private:

src/coreComponents/constitutive/unitTests/TestFluid.hpp

@@ -50,9 +50,8 @@ struct Fluid
   static constexpr integer Vc = 2;    // Critical colume
   static constexpr integer Ac = 3;    // Accentric factor
   static constexpr integer Mw = 4;    // Molecular weight
-  static constexpr integer Vs = 5;    // Volume shift
 
-  static std::array< real64, 66 > data;
+  static std::array< real64, 55 > data;
 };
 
 template< int NC >
@@ -81,7 +80,6 @@ class TestFluid
     createArray( testFluid->criticalVolume, components, Fluid::Vc, Fluid::data );
     createArray( testFluid->acentricFactor, components, Fluid::Ac, Fluid::data );
     createArray( testFluid->molecularWeight, components, Fluid::Mw, Fluid::data );
-    createArray( testFluid->volumeShift, components, Fluid::Vs, Fluid::data );
     testFluid->binaryCoeff.resize( NC, NC );
     return testFluid;
   }
@@ -111,7 +109,6 @@ class TestFluid
       createArray( m_component_properties->m_componentCriticalPressure, criticalPressure );
       createArray( m_component_properties->m_componentCriticalTemperature, criticalTemperature );
       createArray( m_component_properties->m_componentAcentricFactor, acentricFactor );
-      createArray( m_component_properties->m_componentVolumeShift, volumeShift );
       m_component_properties->m_componentBinaryCoeff.resize( NC, NC );
       for( integer ic = 0; ic < NC; ++ic )
       {
@@ -139,7 +136,6 @@ class TestFluid
   array1d< real64 > criticalVolume;
   array1d< real64 > acentricFactor;
   array1d< real64 > molecularWeight;
-  array1d< real64 > volumeShift;
   array2d< real64 > binaryCoeff;
 
 private:
@@ -174,7 +170,7 @@ class TestFluid
   }
 };
 
-std::array< real64, 66 > Fluid::data = {
+std::array< real64, 55 > Fluid::data = {
   // -- Pc
   2.2050e+07, 7.3750e+06, 3.4000e+06, 8.9630e+06, 1.2960e+06, 4.8721e+06,
   4.2481e+06, 3.6400e+06, 4.5990e+06, 2.5300e+06, 1.4600e+06,
@@ -190,9 +186,6 @@ std::array< real64, 66 > Fluid::data = {
   // -- Mw
   1.8015e-02, 4.4010e-02, 2.8013e-02, 3.4100e-02, 1.6043e-02, 3.0070e-02,
   4.4097e-02, 5.8124e-02, 7.2151e-02, 1.1423e-01, 1.4228e-01,
-  // -- Vs
-  0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00,
-  0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00,
 };
 
 }// testing

src/coreComponents/constitutive/unitTests/testCompositionalDensity.cpp

@@ -72,6 +72,10 @@ class CompositionalDensityTestFixture :  public ::testing::TestWithParam< Densit
     string const eosName = EnumStrings< EquationOfStateType >::toString( EOS_TYPE );
     equationOfState->m_equationsOfStateNames.emplace_back( eosName );
 
+    auto densityParameters = const_cast< CompositionalDensity::Parameters * >(m_parameters->get< CompositionalDensity::Parameters >());
+    densityParameters->m_componentVolumeShift.resize( NC );
+    densityParameters->m_componentVolumeShift.zero();
+
     m_density = std::make_unique< CompositionalDensity >( "PhaseDensity", componentProperties, 0, *m_parameters );
   }
 

src/coreComponents/constitutive/unitTests/testCompositionalProperties.cpp

@@ -206,7 +206,8 @@ class CompositionalPropertiesTestDataTestFixture : public ::testing::TestWithPar
   {
     auto const componentProperties = this->m_fluid->createKernelWrapper();
     auto const binaryInteractionCoefficients = componentProperties.m_componentBinaryCoeff;
-    auto const volumeShift = componentProperties.m_componentVolumeShift;
+    stackArray1d< real64, NC > volumeShift( NC );
+    volumeShift.zero();
 
     real64 compressibilityFactor = 0.0;
     stackArray1d< real64, numComps > aPureCoefficient( numComps );

src/coreComponents/constitutive/unitTests/testKValueInitialization.cpp

@@ -85,7 +85,6 @@ class WilsonKValueInitializationTestFixture :
       criticalPressure,
       criticalTemperature,
       acentricFactor,
-      discarded,
       discarded2d );
   }
 

src/coreComponents/constitutive/unitTests/testLohrenzBrayClarkViscosity.cpp

@@ -80,6 +80,10 @@ class LohrenzBrayClarkViscosityTestFixture :  public ::testing::TestWithParam< V
     string const eosName = EnumStrings< EquationOfStateType >::toString( EquationOfStateType::PengRobinson );
     equationOfState->m_equationsOfStateNames.emplace_back( eosName );
 
+    auto densityParameters = const_cast< CompositionalDensity::Parameters * >(m_parameters->get< CompositionalDensity::Parameters >());
+    densityParameters->m_componentVolumeShift.resize( NC );
+    densityParameters->m_componentVolumeShift.zero();
+
     m_density = std::make_unique< CompositionalDensity >( "PhaseDensity", componentProperties, 0, *m_parameters );
     m_viscosity = std::make_unique< LohrenzBrayClarkViscosity >( "PhaseViscosity", componentProperties, 0, *m_parameters );
   }

src/coreComponents/constitutive/unitTests/testNegativeTwoPhaseFlash9Comp.cpp

@@ -258,7 +258,6 @@ std::unique_ptr< TestFluid< numComps > > NegativeTwoPhaseFlashTest9CompFixture<
   TestFluid< numComps >::populateArray( fluid->criticalVolume, Feed< 9 >{9.3999e-05, 9.0001e-05, 9.7999e-05, 1.4800e-04, 2.0000e-04, 2.5800e-04, 3.1000e-04, 3.5100e-04, 6.8243e-04} );
   TestFluid< numComps >::populateArray( fluid->acentricFactor, Feed< 9 >{0.225, 0.04, 0.013, 0.0986, 0.1524, 0.1956, 0.2413, 0.299, 0.5618} );
   TestFluid< numComps >::populateArray( fluid->molecularWeight, Feed< 9 >{44.01e-3, 28.01e-3, 16.04e-3, 30.07e-3, 44.1e-3, 58.12e-3, 72.15e-3, 84e-3, 173e-3} );
-  TestFluid< numComps >::populateArray( fluid->volumeShift, Feed< 9 >{ -0.04958, -0.136012, -0.1486264, -0.10863408, -0.08349872, -0.06331568, -0.04196464, -0.0150072, 0.0000 } );
   fluid->setBinaryCoefficients( Feed< 36 >{
         1.0000e-02,
         0.0000e+00, 3.7320e-03,

src/coreComponents/constitutive/unitTests/testStabilityTest9Comp.cpp

@@ -157,7 +157,6 @@ std::unique_ptr< TestFluid< numComps > > StabilityTestTest9CompFixture< EOS_TYPE
   TestFluid< numComps >::populateArray( fluid->criticalVolume, Feed< 9 >{9.3999e-05, 9.0001e-05, 9.7999e-05, 1.4800e-04, 2.0000e-04, 2.5800e-04, 3.1000e-04, 3.5100e-04, 6.8243e-04} );
   TestFluid< numComps >::populateArray( fluid->acentricFactor, Feed< 9 >{0.225, 0.04, 0.013, 0.0986, 0.1524, 0.1956, 0.2413, 0.299, 0.5618} );
   TestFluid< numComps >::populateArray( fluid->molecularWeight, Feed< 9 >{44.01e-3, 28.01e-3, 16.04e-3, 30.07e-3, 44.1e-3, 58.12e-3, 72.15e-3, 84e-3, 173e-3} );
-  TestFluid< numComps >::populateArray( fluid->volumeShift, Feed< 9 >{ -0.04958, -0.136012, -0.1486264, -0.10863408, -0.08349872, -0.06331568, -0.04196464, -0.0150072, 0.0000 } );
   fluid->setBinaryCoefficients( Feed< 36 >{
         1.0000e-02,
         0.0000e+00, 3.7320e-03,

src/coreComponents/unitTests/constitutiveTests/testMultiFluidCompositionalMultiphase.cpp

@@ -171,7 +171,7 @@ struct Fluid< FluidModel, 5 >
     fill< 5 >( criticalTemperature, {304.1280, 126.1920, 190.5640, 305.3300, 504.2160} );
     array1d< real64 > & acentricFactor = fluid.getReference< array1d< real64 > >( FluidModel::viewKeyStruct::componentAcentricFactorString() );
     fill< 5 >( acentricFactor, {0.223000, 0.037200, 0.010400, 0.099100, 0.250274} );
-    array1d< real64 > & volumeShift = fluid.getReference< array1d< real64 > >( FluidModel::viewKeyStruct::componentVolumeShiftString() );
+    array1d< real64 > & volumeShift = fluid.getReference< array1d< real64 > >( CompositionalDensity::Parameters::viewKeyStruct::componentVolumeShiftString() );
     fill< 5 >( volumeShift, {1.845465e-01, -1.283880e-01, 9.225800e-02, 6.458060e-02, 0.000000e+00} );
     array2d< real64 > & binaryCoeff = fluid.getReference< array2d< real64 > >( FluidModel::viewKeyStruct::componentBinaryCoeffString() );
     fillBinaryCoeffs< 5 >( binaryCoeff, {0.0, 0.1, 0.03, 0.139, 0.032, 0.0, 0.12, 0.03, 0.0, 0.0} );