diff --git a/src/coreComponents/constitutive/fluid/multifluid/MultiFluidBase.hpp b/src/coreComponents/constitutive/fluid/multifluid/MultiFluidBase.hpp index 9396fe6ad3e..2a3e4c4f8e2 100644 --- a/src/coreComponents/constitutive/fluid/multifluid/MultiFluidBase.hpp +++ b/src/coreComponents/constitutive/fluid/multifluid/MultiFluidBase.hpp @@ -1020,27 +1020,16 @@ MultiFluidBase::KernelWrapper:: integer const numPhase = numPhases(); integer const numComp = numComponents(); - StackArray< real64, 4, maxNumDof *maxNumPhase, LAYOUT_PHASE_DC > dPhaseFrac( 1, 1, numPhase, numComp+2 ); - MultiFluidVarSlice< real64, 1, USD_PHASE - 2, USD_PHASE_DC - 2 > - phaseFracAndDeriv { phaseFrac, dPhaseFrac[0][0] }; - - StackArray< real64, 4, maxNumDof *maxNumPhase, LAYOUT_PHASE_DC > dPhaseMassDens( 1, 1, numPhase, numComp+2 ); - MultiFluidVarSlice< real64, 1, USD_PHASE - 2, USD_PHASE_DC - 2 > - phaseMassDensAndDeriv { phaseMassDens, dPhaseMassDens[0][0] }; + // Space for dummy derivatives + StackArray< real64, 4, 2*maxNumDof *maxNumPhase, LAYOUT_PHASE_DC > derivatives( 2, 1, numPhase, numComp+2 ); - StackArray< real64, 4, maxNumDof *maxNumPhase, LAYOUT_PHASE_DC > dPhaseEnthalpy( 1, 1, numPhase, numComp+2 ); - MultiFluidVarSlice< real64, 1, USD_PHASE - 2, USD_PHASE_DC - 2 > - phaseEnthalpyAndDeriv { phaseEnthalpy, dPhaseEnthalpy[0][0] }; - - StackArray< real64, 4, maxNumDof *maxNumPhase, LAYOUT_PHASE_DC > dPhaseInternalEnergy( 1, 1, numPhase, numComp+2 ); - MultiFluidVarSlice< real64, 1, USD_PHASE - 2, USD_PHASE_DC - 2 > - phaseInternalEnergyAndDeriv { phaseInternalEnergy, dPhaseInternalEnergy[0][0] }; + using SliceType = MultiFluidVarSlice< real64, 1, USD_PHASE - 2, USD_PHASE_DC - 2 >; computeInternalEnergy( pressure, - phaseFracAndDeriv, - phaseMassDensAndDeriv, - phaseEnthalpyAndDeriv, - phaseInternalEnergyAndDeriv ); + SliceType { phaseFrac, derivatives[0][0] }, + SliceType { phaseMassDens, derivatives[0][0] }, + SliceType { phaseEnthalpy, derivatives[0][0] }, + SliceType { phaseInternalEnergy, derivatives[1][0] } ); } GEOS_HOST_DEVICE diff --git a/src/coreComponents/constitutive/fluid/multifluid/compositional/CompositionalMultiphaseFluid.cpp b/src/coreComponents/constitutive/fluid/multifluid/compositional/CompositionalMultiphaseFluid.cpp index 56c24735d27..8e7ce268dbf 100644 --- a/src/coreComponents/constitutive/fluid/multifluid/compositional/CompositionalMultiphaseFluid.cpp +++ b/src/coreComponents/constitutive/fluid/multifluid/compositional/CompositionalMultiphaseFluid.cpp @@ -31,8 +31,8 @@ namespace geos namespace constitutive { -template< typename FLASH, typename PHASE1, typename PHASE2, typename PHASE3 > -CompositionalMultiphaseFluid< FLASH, PHASE1, PHASE2, PHASE3 >:: +template< typename FLASH, typename ... PHASES > +CompositionalMultiphaseFluid< FLASH, PHASES... >:: CompositionalMultiphaseFluid( string const & name, Group * const parent ) : MultiFluidBase( name, parent ), m_componentProperties( std::make_unique< compositional::ComponentProperties >( m_componentNames, m_componentMolarWeight ) ), @@ -76,8 +76,8 @@ CompositionalMultiphaseFluid( string const & name, Group * const parent ) .setRestartFlags( RestartFlags::NO_WRITE ); } -template< typename FLASH, typename PHASE1, typename PHASE2, typename PHASE3 > -integer CompositionalMultiphaseFluid< FLASH, PHASE1, PHASE2, PHASE3 >::getWaterPhaseIndex() const +template< typename FLASH, typename ... PHASES > +integer CompositionalMultiphaseFluid< FLASH, PHASES... >::getWaterPhaseIndex() const { auto const phaseTypes = getPhaseTypes(); integer const aqueous = static_cast< integer >(compositional::PhaseType::AQUEOUS); @@ -91,16 +91,16 @@ integer CompositionalMultiphaseFluid< FLASH, PHASE1, PHASE2, PHASE3 >::getWaterP return -1; } -template< typename FLASH, typename PHASE1, typename PHASE2, typename PHASE3 > -string CompositionalMultiphaseFluid< FLASH, PHASE1, PHASE2, PHASE3 >::catalogName() +template< typename FLASH, typename ... PHASES > +string CompositionalMultiphaseFluid< FLASH, PHASES... >::catalogName() { - return GEOS_FMT( "Compositional{}Fluid{}", - FLASH::catalogName(), - PHASE1::Viscosity::catalogName() ); + // Use the first phase viscosity + using ViscosityType = typename std::tuple_element_t< 0, std::tuple< PHASES... > >::Viscosity; + return GEOS_FMT( "Compositional{}Fluid{}", FLASH::catalogName(), ViscosityType::catalogName() ); } -template< typename FLASH, typename PHASE1, typename PHASE2, typename PHASE3 > -void CompositionalMultiphaseFluid< FLASH, PHASE1, PHASE2, PHASE3 >::initializeState() const +template< typename FLASH, typename ... PHASES > +void CompositionalMultiphaseFluid< FLASH, PHASES... >::initializeState() const { // Zero k-Values to force re-initialisation m_kValues.zero(); @@ -108,8 +108,8 @@ void CompositionalMultiphaseFluid< FLASH, PHASE1, PHASE2, PHASE3 >::initializeSt MultiFluidBase::initializeState(); } -template< typename FLASH, typename PHASE1, typename PHASE2, typename PHASE3 > -void CompositionalMultiphaseFluid< FLASH, PHASE1, PHASE2, PHASE3 >::allocateConstitutiveData( Group & parent, localIndex const numPts ) +template< typename FLASH, typename ... PHASES > +void CompositionalMultiphaseFluid< FLASH, PHASES... >::allocateConstitutiveData( Group & parent, localIndex const numPts ) { m_kValues.resize( 0, numPts, numFluidPhases()-1, numFluidComponents() ); @@ -119,8 +119,8 @@ void CompositionalMultiphaseFluid< FLASH, PHASE1, PHASE2, PHASE3 >::allocateCons m_kValues.zero(); } -template< typename FLASH, typename PHASE1, typename PHASE2, typename PHASE3 > -void CompositionalMultiphaseFluid< FLASH, PHASE1, PHASE2, PHASE3 >::postInputInitialization() +template< typename FLASH, typename ... PHASES > +void CompositionalMultiphaseFluid< FLASH, PHASES... >::postInputInitialization() { MultiFluidBase::postInputInitialization(); @@ -200,8 +200,8 @@ void CompositionalMultiphaseFluid< FLASH, PHASE1, PHASE2, PHASE3 >::postInputIni m_parameters->postInputInitialization( this, *m_componentProperties ); } -template< typename FLASH, typename PHASE1, typename PHASE2, typename PHASE3 > -void CompositionalMultiphaseFluid< FLASH, PHASE1, PHASE2, PHASE3 >::initializePostSubGroups() +template< typename FLASH, typename ... PHASES > +void CompositionalMultiphaseFluid< FLASH, PHASES... >::initializePostSubGroups() { MultiFluidBase::initializePostSubGroups(); @@ -209,10 +209,10 @@ void CompositionalMultiphaseFluid< FLASH, PHASE1, PHASE2, PHASE3 >::initializePo createModels(); } -template< typename FLASH, typename PHASE1, typename PHASE2, typename PHASE3 > +template< typename FLASH, typename ... PHASES > std::unique_ptr< ConstitutiveBase > -CompositionalMultiphaseFluid< FLASH, PHASE1, PHASE2, PHASE3 >::deliverClone( string const & name, - Group * const parent ) const +CompositionalMultiphaseFluid< FLASH, PHASES... >::deliverClone( string const & name, + Group * const parent ) const { std::unique_ptr< ConstitutiveBase > clone = MultiFluidBase::deliverClone( name, parent ); CompositionalMultiphaseFluid & newFluid = dynamicCast< CompositionalMultiphaseFluid & >( *clone ); @@ -220,15 +220,21 @@ CompositionalMultiphaseFluid< FLASH, PHASE1, PHASE2, PHASE3 >::deliverClone( str return clone; } -template< typename FLASH, typename PHASE1, typename PHASE2, typename PHASE3 > -typename CompositionalMultiphaseFluid< FLASH, PHASE1, PHASE2, PHASE3 >::KernelWrapper -CompositionalMultiphaseFluid< FLASH, PHASE1, PHASE2, PHASE3 >::createKernelWrapper() +template< typename FLASH, typename ... PHASES > +typename CompositionalMultiphaseFluid< FLASH, PHASES... >::KernelWrapper +CompositionalMultiphaseFluid< FLASH, PHASES... >::createKernelWrapper() +{ + return createKernelWrapper( std::index_sequence_for< PHASES... >{} ); +} + +template< typename FLASH, typename ... PHASES > +template< std::size_t... Is > +typename CompositionalMultiphaseFluid< FLASH, PHASES... >::KernelWrapper +CompositionalMultiphaseFluid< FLASH, PHASES... >::createKernelWrapper( std::index_sequence< Is... > ) { return KernelWrapper( *m_componentProperties, *m_flash, - *m_phase1, - *m_phase2, - *m_phase3, + *std::get< Is >( m_phases )..., m_phaseOrder.toViewConst(), m_componentMolarWeight, m_useMass, @@ -244,8 +250,8 @@ CompositionalMultiphaseFluid< FLASH, PHASE1, PHASE2, PHASE3 >::createKernelWrapp } // Create the fluid models -template< typename FLASH, typename PHASE1, typename PHASE2, typename PHASE3 > -void CompositionalMultiphaseFluid< FLASH, PHASE1, PHASE2, PHASE3 >::createModels() +template< typename FLASH, typename ... PHASES > +void CompositionalMultiphaseFluid< FLASH, PHASES... >::createModels() { m_phaseType = getPhaseTypes(); @@ -258,24 +264,23 @@ void CompositionalMultiphaseFluid< FLASH, PHASE1, PHASE2, PHASE3 >::createModels *m_parameters, m_phaseType.toViewConst() ); - m_phase1 = std::make_unique< PHASE1 >( GEOS_FMT( "{}_PhaseModel1", getName() ), - *m_componentProperties, - 0, - *m_parameters ); - - m_phase2 = std::make_unique< PHASE2 >( GEOS_FMT( "{}_PhaseModel2", getName() ), - *m_componentProperties, - 1, - *m_parameters ); + createPhaseModels( std::index_sequence_for< PHASES... >{} ); +} - m_phase3 = std::make_unique< PHASE3 >( GEOS_FMT( "{}_PhaseModel3", getName() ), - *m_componentProperties, - 2, - *m_parameters ); +template< typename FLASH, typename ... PHASES > +template< std::size_t... Is > +void CompositionalMultiphaseFluid< FLASH, PHASES... >::createPhaseModels( std::index_sequence< Is... > ) +{ + m_phases = std::make_tuple( + std::make_unique< PHASES >( GEOS_FMT( "{}_PhaseModel{}", getName(), Is + 1 ), + *m_componentProperties, + Is, + *m_parameters )... + ); } -template< typename FLASH, typename PHASE1, typename PHASE2, typename PHASE3 > -array1d< integer > CompositionalMultiphaseFluid< FLASH, PHASE1, PHASE2, PHASE3 >::getPhaseTypes() const +template< typename FLASH, typename ... PHASES > +array1d< integer > CompositionalMultiphaseFluid< FLASH, PHASES... >::getPhaseTypes() const { integer const numPhases = numFluidPhases(); array1d< integer > phaseTypes( numPhases ); @@ -287,44 +292,42 @@ array1d< integer > CompositionalMultiphaseFluid< FLASH, PHASE1, PHASE2, PHASE3 > } // Create the fluid models -template< typename FLASH, typename PHASE1, typename PHASE2, typename PHASE3 > +template< typename FLASH, typename ... PHASES > std::unique_ptr< compositional::ModelParameters > -CompositionalMultiphaseFluid< FLASH, PHASE1, PHASE2, PHASE3 >::createModelParameters() +CompositionalMultiphaseFluid< FLASH, PHASES... >::createModelParameters() { std::unique_ptr< compositional::ModelParameters > parameters; parameters = FLASH::createParameters( std::move( parameters )); - parameters = PHASE1::createParameters( std::move( parameters )); - parameters = PHASE2::createParameters( std::move( parameters )); - parameters = PHASE3::createParameters( std::move( parameters )); + ((parameters = PHASES::createParameters( std::move( parameters ) )), ...); return parameters; } // Explicit instantiation of the model template. template class CompositionalMultiphaseFluid< compositional::NegativeTwoPhaseFlashModel, - compositional::PhaseModel< compositional::CompositionalDensity, compositional::ConstantViscosity, compositional::NullModel >, - compositional::PhaseModel< compositional::CompositionalDensity, compositional::ConstantViscosity, compositional::NullModel > >; + compositional::PhaseModel< compositional::CompositionalDensity, compositional::ConstantViscosity >, + compositional::PhaseModel< compositional::CompositionalDensity, compositional::ConstantViscosity > >; template class CompositionalMultiphaseFluid< compositional::NegativeTwoPhaseFlashModel, - compositional::PhaseModel< compositional::CompositionalDensity, compositional::LohrenzBrayClarkViscosity, compositional::NullModel >, - compositional::PhaseModel< compositional::CompositionalDensity, compositional::LohrenzBrayClarkViscosity, compositional::NullModel > >; + compositional::PhaseModel< compositional::CompositionalDensity, compositional::LohrenzBrayClarkViscosity >, + compositional::PhaseModel< compositional::CompositionalDensity, compositional::LohrenzBrayClarkViscosity > >; template class CompositionalMultiphaseFluid< compositional::NegativeTwoPhaseFlashModel, - compositional::PhaseModel< compositional::PhillipsBrineDensity, compositional::PhillipsBrineViscosity, compositional::NullModel >, - compositional::PhaseModel< compositional::CompositionalDensity, compositional::LohrenzBrayClarkViscosity, compositional::NullModel > >; + compositional::PhaseModel< compositional::PhillipsBrineDensity, compositional::PhillipsBrineViscosity >, + compositional::PhaseModel< compositional::CompositionalDensity, compositional::LohrenzBrayClarkViscosity > >; template class CompositionalMultiphaseFluid< compositional::ImmiscibleWaterFlashModel, - compositional::PhaseModel< compositional::CompositionalDensity, compositional::LohrenzBrayClarkViscosity, compositional::NullModel >, - compositional::PhaseModel< compositional::CompositionalDensity, compositional::LohrenzBrayClarkViscosity, compositional::NullModel >, - compositional::PhaseModel< compositional::ImmiscibleWaterDensity, compositional::ImmiscibleWaterViscosity, compositional::NullModel > >; + compositional::PhaseModel< compositional::CompositionalDensity, compositional::LohrenzBrayClarkViscosity >, + compositional::PhaseModel< compositional::CompositionalDensity, compositional::LohrenzBrayClarkViscosity >, + compositional::PhaseModel< compositional::ImmiscibleWaterDensity, compositional::ImmiscibleWaterViscosity > >; template class CompositionalMultiphaseFluid< compositional::KValueFlashModel< 2 >, - compositional::PhaseModel< compositional::CompositionalDensity, compositional::LohrenzBrayClarkViscosity, compositional::NullModel >, - compositional::PhaseModel< compositional::CompositionalDensity, compositional::LohrenzBrayClarkViscosity, compositional::NullModel > >; + compositional::PhaseModel< compositional::CompositionalDensity, compositional::LohrenzBrayClarkViscosity >, + compositional::PhaseModel< compositional::CompositionalDensity, compositional::LohrenzBrayClarkViscosity > >; template class CompositionalMultiphaseFluid< compositional::KValueFlashModel< 2 >, - compositional::PhaseModel< compositional::PhillipsBrineDensity, compositional::PhillipsBrineViscosity, compositional::NullModel >, - compositional::PhaseModel< compositional::CompositionalDensity, compositional::LohrenzBrayClarkViscosity, compositional::NullModel > >; + compositional::PhaseModel< compositional::PhillipsBrineDensity, compositional::PhillipsBrineViscosity >, + compositional::PhaseModel< compositional::CompositionalDensity, compositional::LohrenzBrayClarkViscosity > >; REGISTER_CATALOG_ENTRY( ConstitutiveBase, CompositionalTwoPhaseConstantViscosity, diff --git a/src/coreComponents/constitutive/fluid/multifluid/compositional/CompositionalMultiphaseFluid.hpp b/src/coreComponents/constitutive/fluid/multifluid/compositional/CompositionalMultiphaseFluid.hpp index bd9faf0ce44..8d29973c628 100644 --- a/src/coreComponents/constitutive/fluid/multifluid/compositional/CompositionalMultiphaseFluid.hpp +++ b/src/coreComponents/constitutive/fluid/multifluid/compositional/CompositionalMultiphaseFluid.hpp @@ -43,23 +43,25 @@ namespace constitutive /** * @brief A general compositional fluid model. * @tparam FLASH Class describing the phase equilibrium model - * @tparam PHASE1 Class describing the phase property models for the first phase. - * @tparam PHASE2 Class describing the phase property models for the second phase. - * @tparam PHASE3 Class describing the phase property models for the possible third phase. + * @tparam PHASES Classes describing the phase property models for each phase. */ -template< typename FLASH, typename PHASE1, typename PHASE2, typename PHASE3 = compositional::NullPhaseModel > +template< typename FLASH, typename ... PHASES > class CompositionalMultiphaseFluid : public MultiFluidBase { public: using FlashModel = FLASH; - using Phase1Model = PHASE1; - using Phase2Model = PHASE2; - using Phase3Model = PHASE3; // Get the number of phases - static constexpr integer NUM_PHASES = FlashModel::KernelWrapper::getNumberOfPhases(); - // Currently restrict to 2 or 3 phases - static_assert( NUM_PHASES == 2 || NUM_PHASES == 3 ); + static constexpr integer NUM_PHASES = static_cast< integer >(sizeof...(PHASES)); + + // Ensure that the number of phases matches the flash object + static_assert( NUM_PHASES == FlashModel::KernelWrapper::getNumberOfPhases(), + "Number of phases should match the flash" ); + + // Either all phases are thermal or all are not + static_assert( (std::is_same_v< typename PHASES::Enthalpy, compositional::NullModel > && ...) || + (!std::is_same_v< typename PHASES::Enthalpy, compositional::NullModel > && ...), + "All phase models must either use NullModel for Enthalpy, or none should." ); using exec_policy = parallelDevicePolicy<>; @@ -77,7 +79,10 @@ class CompositionalMultiphaseFluid : public MultiFluidBase virtual void allocateConstitutiveData( dataRepository::Group & parent, localIndex const numPts ) override; - static constexpr bool isThermalType(){ return false; } + static constexpr bool isThermalType() + { + return (!std::is_same_v< typename PHASES::Enthalpy, compositional::NullModel > && ...); + } // TODO: This method should be implemented if an incorrect extrapolation of the pressure and temperature is encountered in the kernel /** @@ -102,7 +107,7 @@ class CompositionalMultiphaseFluid : public MultiFluidBase }; public: - using KernelWrapper = CompositionalMultiphaseFluidUpdates< FLASH, PHASE1, PHASE2, PHASE3 >; + using KernelWrapper = CompositionalMultiphaseFluidUpdates< FLASH, PHASES... >; /** * @brief Create an update kernel wrapper. @@ -120,6 +125,14 @@ class CompositionalMultiphaseFluid : public MultiFluidBase // Create the fluid models void createModels(); + // Create each of the phase models + template< std::size_t... Is > + void createPhaseModels( std::index_sequence< Is... > ); + + // Helper to create the kernel wrappers + template< std::size_t... Is > + KernelWrapper createKernelWrapper( std::index_sequence< Is... > ); + array1d< integer > getPhaseTypes() const; static std::unique_ptr< compositional::ModelParameters > createModelParameters(); @@ -132,9 +145,7 @@ class CompositionalMultiphaseFluid : public MultiFluidBase array1d< integer > m_phaseType{}; // Phase models - std::unique_ptr< PHASE1 > m_phase1{}; - std::unique_ptr< PHASE2 > m_phase2{}; - std::unique_ptr< PHASE3 > m_phase3{}; + std::tuple< std::unique_ptr< PHASES >... > m_phases{}; // Standard EOS component input std::unique_ptr< compositional::ComponentProperties > m_componentProperties{}; @@ -148,29 +159,29 @@ class CompositionalMultiphaseFluid : public MultiFluidBase using CompositionalTwoPhaseConstantViscosity = CompositionalMultiphaseFluid< compositional::NegativeTwoPhaseFlashModel, - compositional::PhaseModel< compositional::CompositionalDensity, compositional::ConstantViscosity, compositional::NullModel >, - compositional::PhaseModel< compositional::CompositionalDensity, compositional::ConstantViscosity, compositional::NullModel > >; + compositional::PhaseModel< compositional::CompositionalDensity, compositional::ConstantViscosity >, + compositional::PhaseModel< compositional::CompositionalDensity, compositional::ConstantViscosity > >; using CompositionalTwoPhaseLohrenzBrayClarkViscosity = CompositionalMultiphaseFluid< compositional::NegativeTwoPhaseFlashModel, - compositional::PhaseModel< compositional::CompositionalDensity, compositional::LohrenzBrayClarkViscosity, compositional::NullModel >, - compositional::PhaseModel< compositional::CompositionalDensity, compositional::LohrenzBrayClarkViscosity, compositional::NullModel > >; + compositional::PhaseModel< compositional::CompositionalDensity, compositional::LohrenzBrayClarkViscosity >, + compositional::PhaseModel< compositional::CompositionalDensity, compositional::LohrenzBrayClarkViscosity > >; using CompositionalTwoPhasePhillipsBrine = CompositionalMultiphaseFluid< compositional::NegativeTwoPhaseFlashModel, - compositional::PhaseModel< compositional::PhillipsBrineDensity, compositional::PhillipsBrineViscosity, compositional::NullModel >, - compositional::PhaseModel< compositional::CompositionalDensity, compositional::LohrenzBrayClarkViscosity, compositional::NullModel > >; + compositional::PhaseModel< compositional::PhillipsBrineDensity, compositional::PhillipsBrineViscosity >, + compositional::PhaseModel< compositional::CompositionalDensity, compositional::LohrenzBrayClarkViscosity > >; using CompositionalThreePhaseLohrenzBrayClarkViscosity = CompositionalMultiphaseFluid< compositional::ImmiscibleWaterFlashModel, - compositional::PhaseModel< compositional::CompositionalDensity, compositional::LohrenzBrayClarkViscosity, compositional::NullModel >, - compositional::PhaseModel< compositional::CompositionalDensity, compositional::LohrenzBrayClarkViscosity, compositional::NullModel >, - compositional::PhaseModel< compositional::ImmiscibleWaterDensity, compositional::ImmiscibleWaterViscosity, compositional::NullModel > >; + compositional::PhaseModel< compositional::CompositionalDensity, compositional::LohrenzBrayClarkViscosity >, + compositional::PhaseModel< compositional::CompositionalDensity, compositional::LohrenzBrayClarkViscosity >, + compositional::PhaseModel< compositional::ImmiscibleWaterDensity, compositional::ImmiscibleWaterViscosity > >; using CompositionalKValueLohrenzBrayClarkViscosity = CompositionalMultiphaseFluid< compositional::KValueFlashModel< 2 >, - compositional::PhaseModel< compositional::CompositionalDensity, compositional::LohrenzBrayClarkViscosity, compositional::NullModel >, - compositional::PhaseModel< compositional::CompositionalDensity, compositional::LohrenzBrayClarkViscosity, compositional::NullModel > >; + compositional::PhaseModel< compositional::CompositionalDensity, compositional::LohrenzBrayClarkViscosity >, + compositional::PhaseModel< compositional::CompositionalDensity, compositional::LohrenzBrayClarkViscosity > >; using CompositionalKValuePhillipsBrine = CompositionalMultiphaseFluid< compositional::KValueFlashModel< 2 >, - compositional::PhaseModel< compositional::PhillipsBrineDensity, compositional::PhillipsBrineViscosity, compositional::NullModel >, - compositional::PhaseModel< compositional::CompositionalDensity, compositional::LohrenzBrayClarkViscosity, compositional::NullModel > >; + compositional::PhaseModel< compositional::PhillipsBrineDensity, compositional::PhillipsBrineViscosity >, + compositional::PhaseModel< compositional::CompositionalDensity, compositional::LohrenzBrayClarkViscosity > >; } /* namespace constitutive */ diff --git a/src/coreComponents/constitutive/fluid/multifluid/compositional/CompositionalMultiphaseFluidUpdates.hpp b/src/coreComponents/constitutive/fluid/multifluid/compositional/CompositionalMultiphaseFluidUpdates.hpp index 6af1f3366dc..fb15d1a7c1c 100644 --- a/src/coreComponents/constitutive/fluid/multifluid/compositional/CompositionalMultiphaseFluidUpdates.hpp +++ b/src/coreComponents/constitutive/fluid/multifluid/compositional/CompositionalMultiphaseFluidUpdates.hpp @@ -21,6 +21,7 @@ #define GEOS_CONSTITUTIVE_FLUID_MULTIFLUID_COMPOSITIONAL_COMPOSITIONALMULTIPHASEFLUIDUPDATES_HPP_ #include "constitutive/fluid/multifluid/compositional/parameters/ComponentProperties.hpp" +#include "constitutive/fluid/multifluid/compositional/models/NullModel.hpp" #include "constitutive/fluid/multifluid/MultiFluidBase.hpp" #include "constitutive/fluid/multifluid/MultiFluidUtils.hpp" @@ -33,19 +34,21 @@ namespace constitutive /** * @brief Kernel wrapper class for CompositionalMultiphaseFluid. * @tparam FLASH Class describing the phase equilibrium model - * @tparam PHASE1 Class describing the phase property models for the first phase. - * @tparam PHASE2 Class describing the phase property models for the second phase. - * @tparam PHASE3 Class describing the phase property models for the possible third phase. + * @tparam PHASES Classes describing the phase property models for each phase. */ -template< typename FLASH, typename PHASE1, typename PHASE2, typename PHASE3 > +template< typename FLASH, typename ... PHASES > class CompositionalMultiphaseFluidUpdates final : public MultiFluidBase::KernelWrapper { +public: + using FlashModel = FLASH; + + // Get the number of phases + static constexpr integer NUM_PHASES = static_cast< integer >(sizeof...(PHASES)); + public: CompositionalMultiphaseFluidUpdates( compositional::ComponentProperties const & componentProperties, FLASH const & flash, - PHASE1 const & phase1, - PHASE2 const & phase2, - PHASE3 const & phase3, + PHASES const & ... phases, arrayView1d< integer const > const & phaseOrder, arrayView1d< real64 const > const & componentMolarWeight, bool const useMass, @@ -94,6 +97,20 @@ class CompositionalMultiphaseFluidUpdates final : public MultiFluidBase::KernelW MultiFluidBase::FluidProp::SliceType const totalDensity, MultiFluidBase::PhaseComp::SliceType::ValueType const & kValues ) const; + /** + * @brief Helper to unpack phase property models and compute phase densities and viscosities. + */ + template< std::size_t... Is > + GEOS_HOST_DEVICE + void computePhaseProperties( real64 const pressure, + real64 const temperature, + MultiFluidBase::PhaseComp::SliceType const phaseComponentFraction, + MultiFluidBase::PhaseProp::SliceType const phaseDensity, + MultiFluidBase::PhaseProp::SliceType const phaseMassDensity, + MultiFluidBase::PhaseProp::SliceType const phaseViscosity, + MultiFluidBase::PhaseProp::SliceType const phaseEnthalpy, + std::index_sequence< Is... > ) const; + /** * @brief Convert derivatives from phase mole fraction to total mole fraction * @details Given property derivatives @c dProperty where composition derivatives are with @@ -116,6 +133,12 @@ class CompositionalMultiphaseFluidUpdates final : public MultiFluidBase::KernelW GEOS_FORCE_INLINE static void setZero( real64 & val ){ val = 0.0; } + GEOS_HOST_DEVICE + static constexpr bool isThermalType() + { + return (!std::is_same_v< typename PHASES::Enthalpy, compositional::NullModel > && ...); + } + private: // The component properties compositional::ComponentProperties::KernelWrapper m_componentProperties; @@ -126,22 +149,18 @@ class CompositionalMultiphaseFluidUpdates final : public MultiFluidBase::KernelW // The ordering of phases arrayView1d< integer const > const m_phaseOrder; - // Phase model kernel wrappers - typename PHASE1::KernelWrapper m_phase1; - typename PHASE2::KernelWrapper m_phase2; - typename PHASE3::KernelWrapper m_phase3; + // Phase model kernel wrappers stored in a tuple + camp::tuple< typename PHASES::KernelWrapper... > m_phases; // Backup variables MultiFluidBase::PhaseComp::ViewValueType m_kValues; }; -template< typename FLASH, typename PHASE1, typename PHASE2, typename PHASE3 > -CompositionalMultiphaseFluidUpdates< FLASH, PHASE1, PHASE2, PHASE3 >:: +template< typename FLASH, typename ... PHASES > +CompositionalMultiphaseFluidUpdates< FLASH, PHASES... >:: CompositionalMultiphaseFluidUpdates( compositional::ComponentProperties const & componentProperties, FLASH const & flash, - PHASE1 const & phase1, - PHASE2 const & phase2, - PHASE3 const & phase3, + PHASES const &... phases, arrayView1d< integer const > const & phaseOrder, arrayView1d< real64 const > const & componentMolarWeight, bool const useMass, @@ -167,17 +186,15 @@ CompositionalMultiphaseFluidUpdates( compositional::ComponentProperties const & m_componentProperties( componentProperties.createKernelWrapper() ), m_flash( flash.createKernelWrapper() ), m_phaseOrder( phaseOrder ), - m_phase1( phase1.createKernelWrapper() ), - m_phase2( phase2.createKernelWrapper() ), - m_phase3( phase3.createKernelWrapper() ), + m_phases( phases.createKernelWrapper()... ), m_kValues( kValues ) {} -template< typename FLASH, typename PHASE1, typename PHASE2, typename PHASE3 > +template< typename FLASH, typename ... PHASES > GEOS_HOST_DEVICE GEOS_FORCE_INLINE void -CompositionalMultiphaseFluidUpdates< FLASH, PHASE1, PHASE2, PHASE3 >::compute( +CompositionalMultiphaseFluidUpdates< FLASH, PHASES... >::compute( real64 const pressure, real64 const temperature, arraySlice1d< real64 const, compflow::USD_COMP - 1 > const & composition, @@ -210,11 +227,62 @@ CompositionalMultiphaseFluidUpdates< FLASH, PHASE1, PHASE2, PHASE3 >::compute( kValues[0][0] ); } -template< typename FLASH, typename PHASE1, typename PHASE2, typename PHASE3 > +template< typename FLASH, typename ... PHASES > +template< std::size_t... Is > +GEOS_HOST_DEVICE +void CompositionalMultiphaseFluidUpdates< FLASH, PHASES... >::computePhaseProperties( + real64 const pressure, + real64 const temperature, + MultiFluidBase::PhaseComp::SliceType const phaseComponentFraction, + MultiFluidBase::PhaseProp::SliceType const phaseDensity, + MultiFluidBase::PhaseProp::SliceType const phaseMassDensity, + MultiFluidBase::PhaseProp::SliceType const phaseViscosity, + MultiFluidBase::PhaseProp::SliceType const phaseEnthalpy, + std::index_sequence< Is... > ) const +{ + // Density computations + ( camp::get< Is >( m_phases ).density.compute( + m_componentProperties, + pressure, + temperature, + phaseComponentFraction.value[m_phaseOrder[Is]].toSliceConst(), + phaseDensity.value[m_phaseOrder[Is]], + phaseDensity.derivs[m_phaseOrder[Is]], + phaseMassDensity.value[m_phaseOrder[Is]], + phaseMassDensity.derivs[m_phaseOrder[Is]], + m_useMass ), ... ); + + // Viscosity computations + ( camp::get< Is >( m_phases ).viscosity.compute( + m_componentProperties, + pressure, + temperature, + phaseComponentFraction.value[m_phaseOrder[Is]].toSliceConst(), + phaseMassDensity.value[m_phaseOrder[Is]], + phaseMassDensity.derivs[m_phaseOrder[Is]].toSliceConst(), + phaseViscosity.value[m_phaseOrder[Is]], + phaseViscosity.derivs[m_phaseOrder[Is]], + m_useMass ), ... ); + + if constexpr (isThermalType()) + { + // Enthalpy computations + ( camp::get< Is >( m_phases ).enthalpy.compute( + m_componentProperties, + pressure, + temperature, + phaseComponentFraction.value[m_phaseOrder[Is]].toSliceConst(), + phaseEnthalpy.value[m_phaseOrder[Is]], + phaseEnthalpy.derivs[m_phaseOrder[Is]], + m_useMass ), ... ); + } +} + +template< typename FLASH, typename ... PHASES > GEOS_HOST_DEVICE GEOS_FORCE_INLINE void -CompositionalMultiphaseFluidUpdates< FLASH, PHASE1, PHASE2, PHASE3 >::compute( +CompositionalMultiphaseFluidUpdates< FLASH, PHASES... >::compute( real64 const pressure, real64 const temperature, arraySlice1d< real64 const, compflow::USD_COMP - 1 > const & composition, @@ -230,9 +298,7 @@ CompositionalMultiphaseFluidUpdates< FLASH, PHASE1, PHASE2, PHASE3 >::compute( { integer constexpr maxNumComp = MultiFluidBase::MAX_NUM_COMPONENTS; integer constexpr maxNumDof = MultiFluidBase::MAX_NUM_COMPONENTS + 2; - integer constexpr maxNumPhase = MultiFluidBase::MAX_NUM_PHASES; integer const numComp = numComponents(); - integer const numPhase = numPhases(); integer const numDof = numComp + 2; // 1. Convert input mass fractions to mole fractions and keep derivatives @@ -264,73 +330,19 @@ CompositionalMultiphaseFluidUpdates< FLASH, PHASE1, PHASE2, PHASE3 >::compute( phaseFraction, phaseComponentFraction ); - // 3. Calculate the phase densities - m_phase1.density.compute( m_componentProperties, - pressure, - temperature, - phaseComponentFraction.value[m_phaseOrder[0]].toSliceConst(), - phaseDensity.value[m_phaseOrder[0]], - phaseDensity.derivs[m_phaseOrder[0]], - phaseMassDensity.value[m_phaseOrder[0]], - phaseMassDensity.derivs[m_phaseOrder[0]], - m_useMass ); - m_phase2.density.compute( m_componentProperties, - pressure, - temperature, - phaseComponentFraction.value[m_phaseOrder[1]].toSliceConst(), - phaseDensity.value[m_phaseOrder[1]], - phaseDensity.derivs[m_phaseOrder[1]], - phaseMassDensity.value[m_phaseOrder[1]], - phaseMassDensity.derivs[m_phaseOrder[1]], - m_useMass ); - if constexpr (2 < FLASH::KernelWrapper::getNumberOfPhases()) - { - m_phase3.density.compute( m_componentProperties, - pressure, - temperature, - phaseComponentFraction.value[m_phaseOrder[2]].toSliceConst(), - phaseDensity.value[m_phaseOrder[2]], - phaseDensity.derivs[m_phaseOrder[2]], - phaseMassDensity.value[m_phaseOrder[2]], - phaseMassDensity.derivs[m_phaseOrder[2]], - m_useMass ); - } - - // 4. Calculate the phase viscosities - m_phase1.viscosity.compute( m_componentProperties, - pressure, - temperature, - phaseComponentFraction.value[m_phaseOrder[0]].toSliceConst(), - phaseMassDensity.value[m_phaseOrder[0]], - phaseMassDensity.derivs[m_phaseOrder[0]].toSliceConst(), - phaseViscosity.value[m_phaseOrder[0]], - phaseViscosity.derivs[m_phaseOrder[0]], - m_useMass ); - m_phase2.viscosity.compute( m_componentProperties, - pressure, - temperature, - phaseComponentFraction.value[m_phaseOrder[1]].toSliceConst(), - phaseMassDensity.value[m_phaseOrder[1]], - phaseMassDensity.derivs[m_phaseOrder[1]].toSliceConst(), - phaseViscosity.value[m_phaseOrder[1]], - phaseViscosity.derivs[m_phaseOrder[1]], - m_useMass ); - if constexpr (2 < FLASH::KernelWrapper::getNumberOfPhases()) - { - m_phase3.viscosity.compute( m_componentProperties, - pressure, - temperature, - phaseComponentFraction.value[m_phaseOrder[2]].toSliceConst(), - phaseMassDensity.value[m_phaseOrder[2]], - phaseMassDensity.derivs[m_phaseOrder[2]].toSliceConst(), - phaseViscosity.value[m_phaseOrder[2]], - phaseViscosity.derivs[m_phaseOrder[2]], - m_useMass ); - } - - // 5. Convert derivatives from phase composition to total composition + // 3. Calculate phase properties: density, viscosity and enthalpy + computePhaseProperties( pressure, + temperature, + phaseComponentFraction, + phaseDensity, + phaseMassDensity, + phaseViscosity, + phaseEnthalpy, + std::make_index_sequence< NUM_PHASES >{} ); + + // 4. Convert derivatives from phase composition to total composition stackArray1d< real64, maxNumDof > workSpace( numDof ); - for( integer ip = 0; ip < FLASH::KernelWrapper::getNumberOfPhases(); ++ip ) + for( integer ip = 0; ip < NUM_PHASES; ++ip ) { convertDerivativesToTotalMoleFraction( numComp, phaseComponentFraction.derivs[ip].toSliceConst(), @@ -344,17 +356,34 @@ CompositionalMultiphaseFluidUpdates< FLASH, PHASE1, PHASE2, PHASE3 >::compute( phaseComponentFraction.derivs[ip].toSliceConst(), phaseViscosity.derivs[ip], workSpace ); + if constexpr (isThermalType()) + { + convertDerivativesToTotalMoleFraction( numComp, + phaseComponentFraction.derivs[ip].toSliceConst(), + phaseEnthalpy.derivs[ip], + workSpace ); + } + } + + // 5. Calculate the internal energy + if constexpr (isThermalType()) + { + computeInternalEnergy( pressure, + phaseFraction, + phaseMassDensity, + phaseEnthalpy, + phaseInternalEnergy ); } // 6. if mass variables used instead of molar, perform the conversion if( m_useMass ) { - real64 phaseMolecularWeight[maxNumPhase]{}; - real64 dPhaseMolecularWeight[maxNumPhase][maxNumDof]{}; + real64 phaseMolecularWeight[NUM_PHASES]{}; + real64 dPhaseMolecularWeight[NUM_PHASES][maxNumDof]{}; arrayView1d< real64 const > const & componentMolarWeight = m_componentProperties.m_componentMolarWeight; - for( integer ip = 0; ip < numPhase; ++ip ) + for( integer ip = 0; ip < NUM_PHASES; ++ip ) { auto const & phaseComposition = phaseComponentFraction.value[ip].toSliceConst(); auto const & dPhaseComposition = phaseComponentFraction.derivs[ip].toSliceConst(); @@ -380,7 +409,7 @@ CompositionalMultiphaseFluidUpdates< FLASH, PHASE1, PHASE2, PHASE3 >::compute( phaseInternalEnergy.derivs ); // Molar density equals mass density - for( integer ip = 0; ip < numPhase; ++ip ) + for( integer ip = 0; ip < NUM_PHASES; ++ip ) { phaseDensity.value[ip] = phaseMassDensity.value[ip]; for( integer idof = 0; idof < numDof; ++idof ) @@ -397,11 +426,11 @@ CompositionalMultiphaseFluidUpdates< FLASH, PHASE1, PHASE2, PHASE3 >::compute( totalDensity ); } -template< typename FLASH, typename PHASE1, typename PHASE2, typename PHASE3 > +template< typename FLASH, typename ... PHASES > GEOS_HOST_DEVICE GEOS_FORCE_INLINE void -CompositionalMultiphaseFluidUpdates< FLASH, PHASE1, PHASE2, PHASE3 >:: +CompositionalMultiphaseFluidUpdates< FLASH, PHASES... >:: update( localIndex const k, localIndex const q, real64 const pressure, @@ -422,11 +451,11 @@ update( localIndex const k, m_kValues[k][q] ); } -template< typename FLASH, typename PHASE1, typename PHASE2, typename PHASE3 > +template< typename FLASH, typename ... PHASES > template< int USD1, int USD2 > GEOS_HOST_DEVICE void -CompositionalMultiphaseFluidUpdates< FLASH, PHASE1, PHASE2, PHASE3 >:: +CompositionalMultiphaseFluidUpdates< FLASH, PHASES... >:: convertDerivativesToTotalMoleFraction( integer const numComps, arraySlice2d< real64 const, USD1 > const & dPhaseComposition, arraySlice1d< real64, USD2 > const & dProperty, diff --git a/src/coreComponents/constitutive/fluid/multifluid/compositional/models/PhaseModel.hpp b/src/coreComponents/constitutive/fluid/multifluid/compositional/models/PhaseModel.hpp index c99b723c0b0..06cc7779df0 100644 --- a/src/coreComponents/constitutive/fluid/multifluid/compositional/models/PhaseModel.hpp +++ b/src/coreComponents/constitutive/fluid/multifluid/compositional/models/PhaseModel.hpp @@ -39,7 +39,7 @@ class ModelParameters; * @tparam VISCOSITY Class describing the viscosity model * @tparam ENTHALPY Class describing the enthalpy model */ -template< typename DENSITY, typename VISCOSITY, typename ENTHALPY > +template< typename DENSITY, typename VISCOSITY, typename ENTHALPY = NullModel > struct PhaseModel { using Density = DENSITY; @@ -85,7 +85,6 @@ struct PhaseModel */ struct KernelWrapper { - /** * @brief Constructor for the kernel wrapper * @param[in] dens the density model @@ -120,7 +119,6 @@ struct PhaseModel /// Kernel wrapper for enthalpy updates typename Enthalpy::KernelWrapper enthalpy; - }; /** @@ -143,12 +141,8 @@ struct PhaseModel phaseParameters = Enthalpy::createParameters( std::move( phaseParameters ) ); return phaseParameters; } - }; -// A no-op phase model -using NullPhaseModel = PhaseModel< NullModel, NullModel, NullModel >; - } // namespace compositional } // namespace constitutive