feat: Immiscible water model

src/coreComponents/constitutive/CMakeLists.txt

@@ -97,8 +97,10 @@ set( constitutive_headers
      fluid/multifluid/compositional/models/CriticalVolume.hpp
      fluid/multifluid/compositional/models/EquationOfState.hpp
      fluid/multifluid/compositional/models/FunctionBase.hpp
+     fluid/multifluid/compositional/models/ImmiscibleWaterDensity.hpp
      fluid/multifluid/compositional/models/ImmiscibleWaterFlashModel.hpp
      fluid/multifluid/compositional/models/ImmiscibleWaterParameters.hpp
+     fluid/multifluid/compositional/models/ImmiscibleWaterViscosity.hpp
      fluid/multifluid/compositional/models/LohrenzBrayClarkViscosity.hpp
      fluid/multifluid/compositional/models/LohrenzBrayClarkViscosityImpl.hpp
      fluid/multifluid/compositional/models/NegativeTwoPhaseFlashModel.hpp
@@ -242,8 +244,10 @@ set( constitutive_sources
      fluid/multifluid/compositional/models/CompositionalDensity.cpp
      fluid/multifluid/compositional/models/ConstantViscosity.cpp
      fluid/multifluid/compositional/models/CriticalVolume.cpp
+     fluid/multifluid/compositional/models/ImmiscibleWaterDensity.cpp
      fluid/multifluid/compositional/models/ImmiscibleWaterFlashModel.cpp
      fluid/multifluid/compositional/models/ImmiscibleWaterParameters.cpp
+     fluid/multifluid/compositional/models/ImmiscibleWaterViscosity.cpp
      fluid/multifluid/compositional/models/LohrenzBrayClarkViscosity.cpp
      fluid/multifluid/compositional/models/NegativeTwoPhaseFlashModel.cpp
      fluid/multifluid/compositional/CompositionalMultiphaseFluid.cpp

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

@@ -53,6 +53,7 @@ void constitutiveUpdatePassThru( constitutive::MultiFluidBase const & fluid,
 #if !defined(GEOS_DEVICE_COMPILE)
                                CO2BrineEzrokhiThermalFluid,
                                CompositionalTwoPhaseLohrenzBrayClarkViscosity,
+                               CompositionalThreePhaseLohrenzBrayClarkViscosity,
 #endif
                                CompositionalTwoPhaseConstantViscosity
                                >::execute( fluid, std::forward< LAMBDA >( lambda ) );
@@ -75,6 +76,7 @@ void constitutiveUpdatePassThru( constitutive::MultiFluidBase & fluid,
 #if !defined(GEOS_DEVICE_COMPILE)
                                CO2BrineEzrokhiThermalFluid,
                                CompositionalTwoPhaseLohrenzBrayClarkViscosity,
+                               CompositionalThreePhaseLohrenzBrayClarkViscosity,
 #endif
                                CompositionalTwoPhaseConstantViscosity
                                >::execute( fluid, std::forward< LAMBDA >( lambda ) );

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

@@ -22,6 +22,7 @@
 #include "constitutive/fluid/multifluid/CO2Brine/functions/PVTFunctionHelpers.hpp"
 #include "constitutive/fluid/multifluid/MultiFluidFields.hpp"
 #include "codingUtilities/Utilities.hpp"
+#include "common/format/StringUtilities.hpp"
 
 namespace geos
 {
@@ -51,6 +52,7 @@ CompositionalMultiphaseFluid( string const & name, Group * const parent )
   m_parameters( createModelParameters() )
 {
   using InputFlags = dataRepository::InputFlags;
+  using RestartFlags = dataRepository::RestartFlags;
 
   getWrapperBase( viewKeyStruct::componentNamesString() ).setInputFlag( InputFlags::REQUIRED );
   getWrapperBase( viewKeyStruct::componentMolarWeightString() ).setInputFlag( InputFlags::REQUIRED );
@@ -80,13 +82,18 @@ CompositionalMultiphaseFluid( string const & name, Group * const parent )
 
   // Link parameters specific to each model
   m_parameters->registerParameters( this );
+
+  // Register extra wrappers to enable auto-cloning
+  registerWrapper( "phaseOrder", &m_phaseOrder )
+    .setSizedFromParent( 0 )
+    .setRestartFlags( RestartFlags::NO_WRITE );
 }
 
 template< typename FLASH, typename PHASE1, typename PHASE2, typename PHASE3 >
 integer CompositionalMultiphaseFluid< FLASH, PHASE1, PHASE2, PHASE3 >::getWaterPhaseIndex() const
 {
-  string const expectedWaterPhaseNames[] = { "water" };
-  return PVTProps::PVTFunctionHelpers::findName( m_phaseNames, expectedWaterPhaseNames, viewKeyStruct::phaseNamesString() );
+  integer const aqueous = static_cast< integer >(PhaseType::AQUEOUS);
+  return m_phaseOrder.size() > aqueous ? m_phaseOrder[aqueous] : -1;
 }
 
 template< typename FLASH, typename PHASE1, typename PHASE2, typename PHASE3 >
@@ -170,6 +177,12 @@ void CompositionalMultiphaseFluid< FLASH, PHASE1, PHASE2, PHASE3 >::postInputIni
     }
   }
 
+  // Determine the phase ordering
+  m_phaseOrder.resize( 3 );
+  m_phaseOrder[PhaseType::LIQUID] = findPhaseIndex( "oil,liq,liquid" );
+  m_phaseOrder[PhaseType::VAPOUR] = findPhaseIndex( "gas,vap,vapor,vapour" );
+  m_phaseOrder[PhaseType::AQUEOUS] = findPhaseIndex( "wat,water,aqueous" );
+
   m_parameters->postInputInitialization( this, *m_componentProperties );
 }
 
@@ -210,6 +223,7 @@ CompositionalMultiphaseFluid< FLASH, PHASE1, PHASE2, PHASE3 >::createKernelWrapp
                         *m_phase1,
                         *m_phase2,
                         *m_phase3,
+                        m_phaseOrder.toViewConst(),
                         m_componentMolarWeight,
                         m_useMass,
                         m_phaseFraction.toView(),
@@ -247,6 +261,22 @@ void CompositionalMultiphaseFluid< FLASH, PHASE1, PHASE2, PHASE3 >::createModels
                                          *m_parameters );
 }
 
+template< typename FLASH, typename PHASE1, typename PHASE2, typename PHASE3 >
+integer CompositionalMultiphaseFluid< FLASH, PHASE1, PHASE2, PHASE3 >::findPhaseIndex( string names ) const
+{
+  auto const nameContainer = stringutilities::tokenize( names, ",", true, false );
+
+  for( integer ip = 0; ip < numFluidPhases(); ++ip )
+  {
+    std::string const phaseName = stringutilities::toLower( m_phaseNames[ip] );
+    if( std::find( nameContainer.begin(), nameContainer.end(), phaseName ) != nameContainer.end())
+    {
+      return ip;
+    }
+  }
+  return -1;
+}
+
 // Create the fluid models
 template< typename FLASH, typename PHASE1, typename PHASE2, typename PHASE3 >
 std::unique_ptr< compositional::ModelParameters >
@@ -269,6 +299,11 @@ template class CompositionalMultiphaseFluid<
     compositional::NegativeTwoPhaseFlashModel,
     compositional::PhaseModel< compositional::CompositionalDensity, compositional::LohrenzBrayClarkViscosity, compositional::NullModel >,
     compositional::PhaseModel< compositional::CompositionalDensity, compositional::LohrenzBrayClarkViscosity, compositional::NullModel > >;
+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 > >;
 
 REGISTER_CATALOG_ENTRY( ConstitutiveBase,
                         CompositionalTwoPhaseConstantViscosity,
@@ -280,6 +315,11 @@ REGISTER_CATALOG_ENTRY( ConstitutiveBase,
                         string const &,
                         dataRepository::Group * const )
 
+REGISTER_CATALOG_ENTRY( ConstitutiveBase,
+                        CompositionalThreePhaseLohrenzBrayClarkViscosity,
+                        string const &,
+                        dataRepository::Group * const )
+
 } // namespace constitutive
 
 } // namespace geos

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

@@ -23,6 +23,9 @@
 #include "constitutive/fluid/multifluid/compositional/CompositionalMultiphaseFluidUpdates.hpp"
 #include "constitutive/fluid/multifluid/compositional/models/ConstantViscosity.hpp"
 #include "constitutive/fluid/multifluid/compositional/models/CompositionalDensity.hpp"
+#include "constitutive/fluid/multifluid/compositional/models/ImmiscibleWaterDensity.hpp"
+#include "constitutive/fluid/multifluid/compositional/models/ImmiscibleWaterFlashModel.hpp"
+#include "constitutive/fluid/multifluid/compositional/models/ImmiscibleWaterViscosity.hpp"
 #include "constitutive/fluid/multifluid/compositional/models/LohrenzBrayClarkViscosity.hpp"
 #include "constitutive/fluid/multifluid/compositional/models/NegativeTwoPhaseFlashModel.hpp"
 #include "constitutive/fluid/multifluid/compositional/models/ModelParameters.hpp"
@@ -110,15 +113,27 @@ class CompositionalMultiphaseFluid : public MultiFluidBase
 
   virtual void resizeFields( localIndex const size, localIndex const numPts ) override;
 
+  enum PhaseType : integer
+  {
+    LIQUID = 0,
+    VAPOUR = 1,
+    AQUEOUS = 2,
+  };
+
 private:
   // Create the fluid models
   void createModels();
 
+  integer findPhaseIndex( string names ) const;
+
   static std::unique_ptr< compositional::ModelParameters > createModelParameters();
 
   // Flash model
   std::unique_ptr< FLASH > m_flash{};
 
+  // Phase ordering
+  array1d< integer > m_phaseOrder;
+
   // Phase models
   std::unique_ptr< PHASE1 > m_phase1{};
   std::unique_ptr< PHASE2 > m_phase2{};
@@ -142,6 +157,11 @@ using CompositionalTwoPhaseLohrenzBrayClarkViscosity = CompositionalMultiphaseFl
   compositional::NegativeTwoPhaseFlashModel,
   compositional::PhaseModel< compositional::CompositionalDensity, compositional::LohrenzBrayClarkViscosity, compositional::NullModel >,
   compositional::PhaseModel< compositional::CompositionalDensity, compositional::LohrenzBrayClarkViscosity, compositional::NullModel > >;
+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 > >;
 
 } /* namespace constitutive */
 

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

@@ -46,6 +46,7 @@ class CompositionalMultiphaseFluidUpdates final : public MultiFluidBase::KernelW
                                        PHASE1 const & phase1,
                                        PHASE2 const & phase2,
                                        PHASE3 const & phase3,
+                                       arrayView1d< integer const > const & phaseOrder,
                                        arrayView1d< real64 const > const & componentMolarWeight,
                                        bool const useMass,
                                        MultiFluidBase::PhaseProp::ViewType phaseFrac,
@@ -122,6 +123,9 @@ class CompositionalMultiphaseFluidUpdates final : public MultiFluidBase::KernelW
   // Flash kernel wrapper
   typename FLASH::KernelWrapper m_flash;
 
+  // The ordering of phases
+  arrayView1d< integer const > const m_phaseOrder;
+
   // Phase model kernel wrappers
   typename PHASE1::KernelWrapper m_phase1;
   typename PHASE2::KernelWrapper m_phase2;
@@ -138,6 +142,7 @@ CompositionalMultiphaseFluidUpdates( compositional::ComponentProperties const &
                                      PHASE1 const & phase1,
                                      PHASE2 const & phase2,
                                      PHASE3 const & phase3,
+                                     arrayView1d< integer const > const & phaseOrder,
                                      arrayView1d< real64 const > const & componentMolarWeight,
                                      bool const useMass,
                                      MultiFluidBase::PhaseProp::ViewType phaseFrac,
@@ -161,6 +166,7 @@ CompositionalMultiphaseFluidUpdates( compositional::ComponentProperties const &
                                  std::move( totalDensity ) ),
   m_componentProperties( componentProperties.createKernelWrapper() ),
   m_flash( flash.createKernelWrapper() ),
+  m_phaseOrder( phaseOrder ),
   m_phase1( phase1.createKernelWrapper() ),
   m_phase2( phase2.createKernelWrapper() ),
   m_phase3( phase3.createKernelWrapper() ),
@@ -262,63 +268,63 @@ CompositionalMultiphaseFluidUpdates< FLASH, PHASE1, PHASE2, PHASE3 >::compute(
   m_phase1.density.compute( m_componentProperties,
                             pressure,
                             temperature,
-                            phaseCompFrac.value[0].toSliceConst(),
-                            phaseDens.value[0],
-                            phaseDens.derivs[0],
-                            phaseMassDensity.value[0],
-                            phaseMassDensity.derivs[0],
+                            phaseCompFrac.value[m_phaseOrder[0]].toSliceConst(),
+                            phaseDens.value[m_phaseOrder[0]],
+                            phaseDens.derivs[m_phaseOrder[0]],
+                            phaseMassDensity.value[m_phaseOrder[0]],
+                            phaseMassDensity.derivs[m_phaseOrder[0]],
                             m_useMass );
   m_phase2.density.compute( m_componentProperties,
                             pressure,
                             temperature,
-                            phaseCompFrac.value[1].toSliceConst(),
-                            phaseDens.value[1],
-                            phaseDens.derivs[1],
-                            phaseMassDensity.value[1],
-                            phaseMassDensity.derivs[1],
+                            phaseCompFrac.value[m_phaseOrder[1]].toSliceConst(),
+                            phaseDens.value[m_phaseOrder[1]],
+                            phaseDens.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,
-                              phaseCompFrac.value[2].toSliceConst(),
-                              phaseDens.value[2],
-                              phaseDens.derivs[2],
-                              phaseMassDensity.value[2],
-                              phaseMassDensity.derivs[2],
+                              phaseCompFrac.value[m_phaseOrder[2]].toSliceConst(),
+                              phaseDens.value[m_phaseOrder[2]],
+                              phaseDens.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,
-                              phaseCompFrac.value[0].toSliceConst(),
-                              phaseMassDensity.value[0],
-                              phaseMassDensity.derivs[0].toSliceConst(),
-                              phaseVisc.value[0],
-                              phaseVisc.derivs[0],
+                              phaseCompFrac.value[m_phaseOrder[0]].toSliceConst(),
+                              phaseMassDensity.value[m_phaseOrder[0]],
+                              phaseMassDensity.derivs[m_phaseOrder[0]].toSliceConst(),
+                              phaseVisc.value[m_phaseOrder[0]],
+                              phaseVisc.derivs[m_phaseOrder[0]],
                               m_useMass );
   m_phase2.viscosity.compute( m_componentProperties,
                               pressure,
                               temperature,
-                              phaseCompFrac.value[1].toSliceConst(),
-                              phaseMassDensity.value[1],
-                              phaseMassDensity.derivs[1].toSliceConst(),
-                              phaseVisc.value[1],
-                              phaseVisc.derivs[1],
+                              phaseCompFrac.value[m_phaseOrder[1]].toSliceConst(),
+                              phaseMassDensity.value[m_phaseOrder[1]],
+                              phaseMassDensity.derivs[m_phaseOrder[1]].toSliceConst(),
+                              phaseVisc.value[m_phaseOrder[1]],
+                              phaseVisc.derivs[m_phaseOrder[1]],
                               m_useMass );
   if constexpr (2 < FLASH::KernelWrapper::getNumberOfPhases())
   {
     m_phase3.viscosity.compute( m_componentProperties,
                                 pressure,
                                 temperature,
-                                phaseCompFrac.value[2].toSliceConst(),
-                                phaseMassDensity.value[2],
-                                phaseMassDensity.derivs[2].toSliceConst(),
-                                phaseVisc.value[2],
-                                phaseVisc.derivs[2],
+                                phaseCompFrac.value[m_phaseOrder[2]].toSliceConst(),
+                                phaseMassDensity.value[m_phaseOrder[2]],
+                                phaseMassDensity.derivs[m_phaseOrder[2]].toSliceConst(),
+                                phaseVisc.value[m_phaseOrder[2]],
+                                phaseVisc.derivs[m_phaseOrder[2]],
                                 m_useMass );
   }