feat: Add Newton solves for stability and flash

src/coreComponents/constitutive/CMakeLists.txt

@@ -101,6 +101,7 @@ set( constitutive_headers
      fluid/multifluid/compositional/functions/SoreideWhitsonEOSPhaseModel_impl.hpp
      fluid/multifluid/compositional/functions/SoreideWhitsonEOSPhaseModel.hpp
      fluid/multifluid/compositional/functions/StabilityTest.hpp
+     fluid/multifluid/compositional/functions/StabilityTest_impl.hpp
      fluid/multifluid/compositional/functions/Utilities.hpp
      fluid/multifluid/compositional/models/CompositionalDensity.hpp
      fluid/multifluid/compositional/models/CompositionalEnthalpy.hpp

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

@@ -38,6 +38,25 @@ GEOS_HOST_DEVICE
 GEOS_FORCE_INLINE
 static void setZero( real64 & value ){ value = 0.0; }
 
+/**
+ * @brief Compies a composition from one array to another
+ * @param[in] numComps number of components
+ * @param[in] inputComposition input composition to be copied
+ * @param[out] outputComposition output composition to be copied into
+ */
+template< integer USD1, integer USD2 >
+GEOS_HOST_DEVICE
+GEOS_FORCE_INLINE
+static void copyComposition( integer const numComps,
+                             arraySlice1d< real64 const, USD1 > const & inputComposition,
+                             arraySlice1d< real64, USD2 > const & outputComposition )
+{
+  for( integer ic = 0; ic < numComps; ++ic )
+  {
+    outputComposition[ic] = inputComposition[ic];
+  }
+}
+
 /**
  * @brief Normalise a composition in place to ensure that the components add up to unity
  * @param[in] numComps number of components
@@ -63,6 +82,31 @@ static real64 normalizeComposition( integer const numComps,
   return totalMoles;
 }
 
+/**
+ * @brief Find the reference component - the index of the component with the largest composition
+ * @param[in] numComps number of components
+ * @param[in] composition composition to be assessed
+ * @return the index of the component with the largest composition
+ */
+template< integer USD >
+GEOS_HOST_DEVICE
+GEOS_FORCE_INLINE
+static integer getReferenceComponent( integer const numComps,
+                                      arraySlice1d< real64 const, USD > const & composition )
+{
+  real64 maxComposition = -LvArray::NumericLimits< real64 >::max;
+  integer refComponent = -1;
+  for( integer ic = 0; ic < numComps; ++ic )
+  {
+    if( maxComposition < composition[ic] )
+    {
+      refComponent = ic;
+      maxComposition = composition[ic];
+    }
+  }
+  return refComponent;
+}
+
 /**
  * @brief Checks if a composition array has numerically small values
  * @param[in] numComps number of values in the composition array
@@ -113,6 +157,34 @@ static integer calculatePresentComponents( integer const numComps,
   return presentCount;
 }
 
+/**
+ * @brief Calculate which components are absent.
+ * @details Creates a list of indices whose components have zero mole fraction.
+ * @param[in] numComps number of components
+ * @param[in] composition the composition of the fluid
+ * @param[out] absentComponents the list of absent components
+ * @return the number of absent components
+ */
+template< integer USD, typename ARRAY >
+GEOS_HOST_DEVICE
+GEOS_FORCE_INLINE
+static integer calculateAbsentComponents( integer const numComps,
+                                          arraySlice1d< real64 const, USD > const & composition,
+                                          ARRAY & absentComponents )
+{
+  // Check for machine-zero feed values
+  integer absentCount = 0;
+  for( integer ic = 0; ic < numComps; ++ic )
+  {
+    if( composition[ic] < MultiFluidConstants::minForSpeciesPresence )
+    {
+      absentComponents[absentCount++] = ic;
+    }
+  }
+  absentComponents.resize( absentCount );
+  return absentCount;
+}
+
 /**
  * @brief Determines whether a given composition represents a pure or nearly pure mixture.
  *
@@ -157,15 +229,15 @@ static bool checkPureMixture( integer const numComps,
 }
 
 /**
- * @brief Solve the lineat system for the derivatives of the flash
+ * @brief Solve the linear system for the derivatives of the flash
  * @param[in/out] A the coefficient matrix. Destroyed after call
  * @param[in/out] X the rhs and solution
  * @return @c true if the problem is well solved @c false otherwise
  */
-template< int USD >
+template< int USD, typename RHS_TYPE >
 GEOS_HOST_DEVICE
 static bool solveLinearSystem( arraySlice2d< real64, USD > const & A,
-                               arraySlice2d< real64, USD > const & X )
+                               RHS_TYPE const & X )
 {
 #if defined(GEOS_DEVICE_COMPILE)
   return denseLinearAlgebra::solveGaussianElimination( A, X );

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

@@ -87,9 +87,15 @@ class NegativeTwoPhaseFlashModelUpdate final : public FunctionBaseUpdate
     integer almostPureComponent = -1;
     checkPureMixture( m_numComponents, compFraction, pureComponent, almostPureComponent );
 
+    for( integer ic = 0; ic < m_numComponents; ++ic )
+    {
+      incipientComposition[ic] = compFraction[ic];
+    }
+
     // Check if k-Values need to be initialised
     auto kVapourLiquid = kValues[0];
     bool const needInitialisation = (pureComponent < 0) && hasZero( m_numComponents, kVapourLiquid.toSliceConst() );
+
     if( needInitialisation )
     {
       if( m_flashData.liquidEos == EquationOfStateType::SoreideWhitson )

src/coreComponents/constitutive/fluid/multifluid/compositional/parameters/FlashParameters.cpp

@@ -21,6 +21,7 @@
 #include "ComponentProperties.hpp"
 
 #include "constitutive/fluid/multifluid/MultiFluidBase.hpp"
+#include "common/MpiWrapper.hpp"
 
 namespace geos
 {
@@ -39,7 +40,8 @@ FlashParameters::FlashParameters( std::unique_ptr< ModelParameters > parameters
   m_continuousParameters[STABILITY_THRESHOLD] = -1.0e-8;
   m_continuousParameters[STABILITY_TOLERANCE] = 1.0e-8;
   m_continuousParameters[FLASH_TOLERANCE] = 1.0e-8;
-  m_continuousParameters[NEGATIVE_FLASH_TOLERANCE] = -1.0e-1;
+  m_continuousParameters[NEGATIVE_FLASH_TOLERANCE] = 1.0e-5;
+  m_continuousParameters[SSI_TOLERANCE] = 1.0e-3;
 
   m_discreteParameters[STABILITY_MAX_ITERATIONS] = 300;
   m_discreteParameters[FLASH_MAX_ITERATIONS] = 300;
@@ -69,7 +71,7 @@ void FlashParameters::registerParametersImpl( MultiFluidBase * fluid )
   fluid->registerWrapper( viewKeyStruct::stabilityMaxIterationsString(), &m_discreteParameters[STABILITY_MAX_ITERATIONS] ).
     setInputFlag( dataRepository::InputFlags::OPTIONAL ).
     setDefaultValue( m_discreteParameters[STABILITY_MAX_ITERATIONS] ).
-    setDescription( "The maximum number of successive substitution iterations used during the stability test. "
+    setDescription( "The maximum number of iterations used during the stability test. "
                     "A zero or negative value skips the stability test" );
 
   fluid->registerWrapper( viewKeyStruct::flashMaxIterationsString(), &m_discreteParameters[FLASH_MAX_ITERATIONS] ).
@@ -81,6 +83,17 @@ void FlashParameters::registerParametersImpl( MultiFluidBase * fluid )
     setInputFlag( dataRepository::InputFlags::OPTIONAL ).
     setDefaultValue( m_continuousParameters[FLASH_TOLERANCE] ).
     setDescription( "The tolerance to use to determine convergece of the flash" );
+
+  fluid->registerWrapper( viewKeyStruct::ssiToleranceString(), &m_continuousParameters[SSI_TOLERANCE] ).
+    setInputFlag( dataRepository::InputFlags::OPTIONAL ).
+    setDefaultValue( m_continuousParameters[SSI_TOLERANCE] ).
+    setDescription( "The tolerance to use to switch from SSI iterations to Newton iterations" );
+
+  fluid->registerWrapper( viewKeyStruct::negativeFlashToleranceString(), &m_continuousParameters[NEGATIVE_FLASH_TOLERANCE] ).
+    setInputFlag( dataRepository::InputFlags::OPTIONAL ).
+    setDefaultValue( m_continuousParameters[NEGATIVE_FLASH_TOLERANCE] ).
+    setDescription( "The tolerance within the physical region [0,1] which the phase fraction must be in "
+                    "order to allow a switch from SSI iterations to Newton iterations" );
 }
 
 void FlashParameters::postInputInitializationImpl( MultiFluidBase const * fluid,
@@ -120,6 +133,28 @@ void FlashParameters::postInputInitializationImpl( MultiFluidBase const * fluid,
   real64 const flashTolerance = m_continuousParameters[FLASH_TOLERANCE];
   checkLowerBound( flashTolerance, epsilon, viewKeyStruct::flashToleranceString() );
 
+  // SSI tolerance should be positive
+  real64 const ssiTolerance = m_continuousParameters[SSI_TOLERANCE];
+  checkLowerBound( ssiTolerance, epsilon, viewKeyStruct::ssiToleranceString() );
+
+  // If the SSI tolerance is lower than the stability tolerance then it will be ignored for stability
+  if( ssiTolerance < stabilityTolerance )
+  {
+    GEOS_WARNING_IF( MpiWrapper::commRank() == 0,
+                     GEOS_FMT( "The SSI tolerance {} is less than the stability tolerance {}. "
+                               "The SSI tolerance will be ignored for stability calculations.",
+                               ssiTolerance, stabilityTolerance ));
+  }
+
+  // If the SSI tolerance is lower than the flash tolerance then it will be ignored for flash
+  if( ssiTolerance < flashTolerance )
+  {
+    GEOS_WARNING_IF( MpiWrapper::commRank() == 0,
+                     GEOS_FMT( "The SSI tolerance {} is less than the flash tolerance {}. "
+                               "The SSI tolerance will be ignored for flash calculations.",
+                               ssiTolerance, flashTolerance ));
+  }
+
   // Max number of flash iterations should be at least 1
   integer const flashMaxIterations = m_discreteParameters[FLASH_MAX_ITERATIONS];
   checkLowerBound( flashMaxIterations, 1, viewKeyStruct::flashMaxIterationsString() );

src/coreComponents/constitutive/fluid/multifluid/compositional/parameters/FlashParameters.hpp

@@ -63,8 +63,11 @@ class FlashParameters : public ModelParameters
   /** The tolerance to use to determine convergece to a stationary point in the stability test. */
   static integer constexpr STABILITY_TOLERANCE = 3;
 
+  /** The tolerance to use to switch from SSI iterations to Newton iterations. */
+  static integer constexpr SSI_TOLERANCE = 4;
+
   /* Number of continuous variables */
-  static integer constexpr NUM_FLOAT = 4;
+  static integer constexpr NUM_FLOAT = 5;
 
 public:
   FlashParameters( std::unique_ptr< ModelParameters > parameters );
@@ -78,7 +81,9 @@ class FlashParameters : public ModelParameters
     static constexpr char const * flashMaxIterationsString() { return "flashMaxIterations"; }
     static constexpr char const * stabilityThresholdString() { return "stabilityThreshold"; }
     static constexpr char const * stabilityToleranceString() { return "stabilityTolerance"; }
+    static constexpr char const * ssiToleranceString() { return "ssiTolerance"; }
     static constexpr char const * stabilityMaxIterationsString() { return "stabilityMaxIterations"; }
+    static constexpr char const * negativeFlashToleranceString() { return "negativeFlashTolerance"; }
   };
 
   array1d< real64 > m_continuousParameters;

src/coreComponents/constitutive/unitTests/CMakeLists.txt

@@ -30,6 +30,7 @@ set( gtest_geosx_tests
      testMultiFluidSelector.cpp
      testNegativeTwoPhaseFlash.cpp
      testNegativeTwoPhaseFlash9Comp.cpp
+     testNegativeTwoPhaseFlashJacobian.cpp
      testNegativeTwoPhaseFlashModel.cpp
      testParticleFluidEnums.cpp
      testPressureTemperatureCoordinates.cpp
@@ -39,6 +40,8 @@ set( gtest_geosx_tests
      testSoreideWhitsonFlash.cpp
      testStabilityTest2Comp.cpp
      testStabilityTest9Comp.cpp
+     testStabilityTestDerivatives.cpp
+     testStabilityTestJacobian.cpp
      testRachfordRice.cpp
      testTwoPhaseImmiscibleFluid.cpp )
 

src/coreComponents/constitutive/unitTests/testNegativeTwoPhaseFlash9Comp.cpp

@@ -90,6 +90,7 @@ class NegativeTwoPhaseFlashTest9CompFixture :  public ::testing::TestWithParam<
 
     auto const parameters = FlashParameters::create( std::make_unique< ModelParameters >() );
     auto const * flashParameters = parameters->get< FlashParameters >();
+    flashParameters->m_discreteParameters[FlashParameters::FLASH_MAX_ITERATIONS] = 100;
 
     bool status = NegativeTwoPhaseFlash::compute(
       numComps,