diff --git a/src/coreComponents/constitutive/CMakeLists.txt b/src/coreComponents/constitutive/CMakeLists.txt index 39ddb39771b..86e6c4e7a70 100644 --- a/src/coreComponents/constitutive/CMakeLists.txt +++ b/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 diff --git a/src/coreComponents/constitutive/fluid/multifluid/compositional/functions/NegativeTwoPhaseFlash.hpp b/src/coreComponents/constitutive/fluid/multifluid/compositional/functions/NegativeTwoPhaseFlash.hpp index 1509701600a..0940f158ebc 100644 --- a/src/coreComponents/constitutive/fluid/multifluid/compositional/functions/NegativeTwoPhaseFlash.hpp +++ b/src/coreComponents/constitutive/fluid/multifluid/compositional/functions/NegativeTwoPhaseFlash.hpp @@ -32,38 +32,62 @@ namespace constitutive namespace compositional { +/** + * @brief Struct implementing a negative two-phase flash calculation. + * + * @details This struct provides methods to perform a two-phase flash calculation + * for reservoir fluid systems where the mixture is thermodynamically unstable. + * It uses the algorithm described in Michelsen (1982) + * + * The flash calculation is performed in two stages: + * - Successive Substitution Iteration (SSI): Used for initial convergence. + * - Newton-Raphson Iteration: Used to refine the solution and accelerate convergence. + * + * The algorithm solves the Rachford-Rice equation to determine the vapor phase mole fraction, + * and iteratively updates K-values and phase compositions to satisfy fugacity equilibrium. + * + * Missing components (i.e., those with zero mole fraction in the feed) are excluded from + * fugacity and Jacobian calculations but handled explicitly in the residual system. + * + * Reference: Michelsen, M.L. (1982). The Isothermal Flash Problem. Part II. Phase Split Calculation. + * Fluid Phase Equilibria, 9(1), 21–40. https://doi.org/10.1016/0378-3812(82)85002-4 + */ struct NegativeTwoPhaseFlash { using Deriv = constitutive::multifluid::DerivativeOffset; public: /** - * @brief Perform a two-phase EOS flash calculation. + * @brief Perform a negative two-phase flash calculation for unstable mixtures. * - * @details - * Executes a negative flash calculation using an equation of state (EOS) model - * to determine the equilibrium phase split of a multicomponent mixture at a given - * pressure and temperature. The function computes the vapor phase mole fraction, - * as well as the compositions of the vapor and liquid phases, based on the provided - * initial estimates and component properties. Uses successive substitution. + * @details This method solves the isothermal two-phase flash problem for mixtures + * that are thermodynamically unstable. + * + * The algorithm proceeds in two stages: + * - Successive Substitution Iteration (SSI): Used for initial convergence of K-values. + * - Newton-Raphson Iteration: Refines the solution using a Jacobian matrix and residual vector. * - * @tparam USD1 Unique stride descriptor for the K-values array. - * @tparam USD2 Unique stride descriptor for the phase composition arrays. + * The Jacobian matrix includes partial derivatives of fugacity with respect to composition + * and vapor fraction. Missing components (zero mole fraction) are handled by setting + * diagonal Jacobian entries to unity to ensure numerical stability. + * + * @tparam USD1 Unit Stride Dimension for the K-values array. + * @tparam USD2 Unit Stride Dimension for the composition arrays. * * @param[in] numComps Number of components in the mixture. * @param[in] pressure System pressure [Pa]. * @param[in] temperature System temperature [K]. - * @param[in] composition Overall composition of the mixture (z_i). + * @param[in] composition Overall feed composition (mole fractions). * @param[in] componentProperties Thermodynamic properties of each component. - * @param[in] flashData Precomputed data and parameters required for the flash calculation. - * @param[in] continuousFlashParameters List of continuous (floating-point) parameters for the flash. - * @param[in] discreteFlashParameters List of discrete (integer) parameters for the flash. - * @param[in,out] kValues Phase equilibrium ratios (K-values), updated during the flash. - * @param[out] vapourPhaseMoleFraction Calculated vapor phase mole fraction (V). - * @param[out] liquidComposition Calculated composition of the liquid phase (x_i). - * @param[out] vapourComposition Calculated composition of the vapor phase (y_i). - * - * @return True if the flash calculation was successful; false otherwise. + * @param[in] flashData Flash calculation context including EOS types. + * @param[in] continuousFlashParameters Continuous parameters for flash calculation. + * @param[in] discreteFlashParameters Discrete parameters for flash calculation. + * @param[in,out] kValues Initial and updated K-values (vapor-liquid equilibrium ratios). + * @param[out] vapourPhaseMoleFraction Computed vapor phase mole fraction. + * @param[out] liquidComposition Computed liquid phase composition. + * @param[out] vapourComposition Computed vapor phase composition. + * + * @return True if the flash calculation converged; false otherwise. */ template< int USD1, int USD2 > GEOS_HOST_DEVICE @@ -81,19 +105,29 @@ struct NegativeTwoPhaseFlash arraySlice1d< real64, USD2 > const & vapourComposition ); /** - * @brief Calculate derivatives from the two-phase negative flash - * @param[in] numComps number of components - * @param[in] pressure pressure - * @param[in] temperature temperature - * @param[in] composition composition of the mixture - * @param[in] componentProperties The compositional component properties - * @param[in] flashData The parameters required for the flash - * @param[in] vapourFraction the calculated vapour (gas) mole fraction - * @param[in] liquidComposition the calculated liquid phase composition - * @param[in] vapourComposition the calculated vapour phase composition - * @param[out] vapourFractionDerivs derivatives of the calculated vapour (gas) mole fraction - * @param[out] liquidCompositionDerivs derivatives of the calculated liquid phase composition - * @param[out] vapourCompositionDerivs derivatives of the calculated vapour phase composition + * @brief Compute derivatives of the vapor fraction and phase compositions for a negative two-phase flash. + * + * @details This method computes the sensitivities of the vapor fraction and the liquid/vapor phase compositions + * with respect to pressure, temperature, and overall composition. It uses implicit differentiation of the + * equilibrium conditions and mass balance constraints to solve for the derivatives. The Jacobian matrix is + * constructed from the linearized system of equations derived from fugacity equilibrium and Rachford-Rice relations. + * + * @tparam USD1 Unit Stride Dimension for input phase composition arrays. + * @tparam USD2 Unit Stride Dimension for vapor fraction derivative array. + * @tparam USD3 Unit Stride Dimension for phase composition derivative arrays. + * + * @param[in] numComps Number of components in the mixture. + * @param[in] pressure System pressure [Pa]. + * @param[in] temperature System temperature [K]. + * @param[in] composition Overall feed composition. + * @param[in] componentProperties Thermodynamic properties of each component. + * @param[in] flashData Flash calculation context including EOS types. + * @param[in] vapourFraction Vapor phase mole fraction. + * @param[in] liquidComposition Composition of the liquid phase. + * @param[in] vapourComposition Composition of the vapor phase. + * @param[out] vapourFractionDerivs Derivatives of vapor fraction with respect to pressure, temperature, and composition. + * @param[out] liquidCompositionDerivs Derivatives of liquid phase composition. + * @param[out] vapourCompositionDerivs Derivatives of vapor phase composition. */ template< integer USD1, integer USD2, integer USD3 > GEOS_HOST_DEVICE @@ -110,6 +144,112 @@ struct NegativeTwoPhaseFlash arraySlice2d< real64, USD3 > const & liquidCompositionDerivs, arraySlice2d< real64, USD3 > const & vapourCompositionDerivs ); + /** + * @brief Compute the residual vector for a two-phase flash calculation. + * + * @details This method evaluates the phase equilibrium condition by computing + * the residuals of fugacity balance between liquid and vapor phases for each component. + * It also includes a mass balance residual for the vapor phase mole fraction. + * The compositions are calculated using the Rachford-Rice formulation and normalized. + * + * The fugacity coefficients are computed for both phases using the specified equation of state. + * The residuals are defined as: + * R_i = x_i * \phi_i^L - y_i * \phi_i^V + * R_N = \sum_{i=1}^N (y_i - x_i) + * where \phi_i^L and \phi_i^V are fugacity coefficients in liquid and vapor phases. + * + * @tparam USD1 Unit Stride Dimension for k-values array. + * @tparam USD2 Unit Stride Dimension for composition arrays. + * @tparam USD3 Unit Stride Dimension for residual array. + * + * @param[in] numComps Number of components in the mixture. + * @param[in] pressure System pressure [Pa]. + * @param[in] temperature System temperature [K]. + * @param[in] composition Overall feed composition (mole fractions). + * @param[in] componentProperties Thermodynamic properties of components. + * @param[in] flashData Flash calculation context including EOS types. + * @param[in] kValues Equilibrium ratios (K-values) for each component. + * @param[in] vapourPhaseMoleFraction Mole fraction of the vapor phase. + * @param[out] liquidComposition Computed composition of the liquid phase. + * @param[out] vapourComposition Computed composition of the vapor phase. + * @param[out] logLiquidFugacity Log fugacity coefficients for the liquid phase. + * @param[out] logVapourFugacity Log fugacity coefficients for the vapor phase. + * @param[out] residual Residual vector containing fugacity and mass balance errors. + * + * @return L2 norm of the residual vector. + */ + template< int USD1, int USD2, int USD3 > + GEOS_HOST_DEVICE + static real64 calculateResidual( integer const numComps, + real64 const pressure, + real64 const temperature, + arraySlice1d< real64 const > const & composition, + ComponentProperties::KernelWrapper const & componentProperties, + FlashData const & flashData, + arraySlice1d< real64 const, USD1 > const & kValues, + real64 const & vapourPhaseMoleFraction, + arraySlice1d< real64, USD2 > const & liquidComposition, + arraySlice1d< real64, USD2 > const & vapourComposition, + arraySlice1d< real64 > const & logLiquidFugacity, + arraySlice1d< real64 > const & logVapourFugacity, + arraySlice1d< real64, USD3 > const & residual ); + + /** + * @brief Calculates the residual vector and Jacobian matrix for a two-phase flash system. + * + * @details This method computes the nonlinear residuals and Jacobian matrix for a negative two-phase flash + * system using fugacity equilibrium and mass balance constraints. It is used in Newton-Raphson iterations + * to refine the solution of the flash problem. The primary variables in the Jacobian are the logarithm of + * the K-values (log(K)) and the vapor phase mole fraction. The Jacobian is constructed with respect to these + * variables. Missing components (i.e., those with zero mole fraction in the feed) are excluded from the + * fugacity and Jacobian calculations but are handled explicitly by setting diagonal Jacobian entries to unity. + * + * The method normalizes phase compositions, computes fugacity coefficients and their derivatives, and + * assembles the Jacobian matrix using chain rule and normalization sensitivities. + * + * @tparam USD1 Unit Stride Dimension for the K-values array. + * @tparam USD2 Unit Stride Dimension for the composition arrays. + * @tparam USD3 Unit Stride Dimension for the residual vector. + * @tparam USD4 Unit Stride Dimension for the Jacobian matrix. + * + * @param[in] numComps Number of components in the mixture. + * @param[in] pressure System pressure [Pa]. + * @param[in] temperature System temperature [K]. + * @param[in] composition Overall feed composition (mole fractions). + * @param[in] componentProperties Thermodynamic properties of each component. + * @param[in] flashData Flash calculation context including EOS types. + * @param[in] kValues Equilibrium ratios (K-values) for each component. + * @param[in] vapourPhaseMoleFraction Mole fraction of the vapor phase. + * @param[in] liquidComposition Composition of the liquid phase. + * @param[in] vapourComposition Composition of the vapor phase. + * @param[in] logLiquidFugacity Logarithm of fugacity coefficients in the liquid phase. + * @param[in] logVapourFugacity Logarithm of fugacity coefficients in the vapor phase. + * @param[in] logLiquidFugacityDerivs Derivatives of log fugacity in the liquid phase. + * @param[in] logVapourFugacityDerivs Derivatives of log fugacity in the vapor phase. + * @param[out] residual Residual vector containing fugacity and mass balance errors. + * @param[out] jacobian Jacobian matrix of the residuals with respect to log(K) and vapor fraction. + * + * @return Error norm of the residual vector. + */ + template< int USD1, int USD2, int USD3, int USD4 > + GEOS_HOST_DEVICE + static real64 calculateResidualAndJacobian( integer const numComps, + real64 const pressure, + real64 const temperature, + arraySlice1d< real64 const > const & composition, + ComponentProperties::KernelWrapper const & componentProperties, + FlashData const & flashData, + arraySlice1d< real64 const, USD1 > const & kValues, + real64 const & vapourPhaseMoleFraction, + arraySlice1d< real64, USD2 > const & liquidComposition, + arraySlice1d< real64, USD2 > const & vapourComposition, + arraySlice1d< real64 > const & logLiquidFugacity, + arraySlice1d< real64 > const & logVapourFugacity, + arraySlice2d< real64 > const & logLiquidFugacityDerivs, + arraySlice2d< real64 > const & logVapourFugacityDerivs, + arraySlice1d< real64, USD3 > const & residual, + arraySlice2d< real64, USD4 > const & jacobian ); + private: /** * @brief Calculate the derivatives of the flash. @@ -138,19 +278,18 @@ struct NegativeTwoPhaseFlash */ template< integer USD1, integer USD2, integer USD3 > GEOS_HOST_DEVICE - static void computeDerivatives( - integer const numComps, - arraySlice1d< real64 const > const & totalComposition, - real64 const & phase1Fraction, - arraySlice1d< real64 const, USD1 > const & phase1Composition, - arraySlice1d< real64 const, USD1 > const & phase2Composition, - arraySlice1d< real64 const > const & phase1Fugacity, - arraySlice1d< real64 const > const & phase2Fugacity, - arraySlice2d< real64 const > const & phase1LogFugacityDerivs, - arraySlice2d< real64 const > const & phase2LogFugacityDerivs, - arraySlice1d< real64, USD2 > const & phase1FractionDerivs, - arraySlice2d< real64, USD3 > const & phase1CompositionDerivs, - arraySlice2d< real64, USD3 > const & phase2CompositionDerivs ); + static void computeDerivatives( integer const numComps, + arraySlice1d< real64 const > const & totalComposition, + real64 const & phase1Fraction, + arraySlice1d< real64 const, USD1 > const & phase1Composition, + arraySlice1d< real64 const, USD1 > const & phase2Composition, + arraySlice1d< real64 const > const & phase1Fugacity, + arraySlice1d< real64 const > const & phase2Fugacity, + arraySlice2d< real64 const > const & phase1LogFugacityDerivs, + arraySlice2d< real64 const > const & phase2LogFugacityDerivs, + arraySlice1d< real64, USD2 > const & phase1FractionDerivs, + arraySlice2d< real64, USD3 > const & phase1CompositionDerivs, + arraySlice2d< real64, USD3 > const & phase2CompositionDerivs ); /** * @brief A test for when the result of the negative flash should be truncated immediately diff --git a/src/coreComponents/constitutive/fluid/multifluid/compositional/functions/NegativeTwoPhaseFlash_impl.hpp b/src/coreComponents/constitutive/fluid/multifluid/compositional/functions/NegativeTwoPhaseFlash_impl.hpp index f251d632117..8b8faf70ae6 100644 --- a/src/coreComponents/constitutive/fluid/multifluid/compositional/functions/NegativeTwoPhaseFlash_impl.hpp +++ b/src/coreComponents/constitutive/fluid/multifluid/compositional/functions/NegativeTwoPhaseFlash_impl.hpp @@ -50,20 +50,29 @@ bool NegativeTwoPhaseFlash::compute( integer const numComps, arraySlice1d< real64, USD2 > const & vapourComposition ) { constexpr integer maxNumComps = MultiFluidConstants::MAX_NUM_COMPONENTS; - StackArray< real64, 2, 3*maxNumComps > workSpace( 3, maxNumComps ); + constexpr integer maxNumDofs = MultiFluidConstants::MAX_NUM_COMPONENTS + 2; + + integer const numDofs = numComps + 2; + + StackArray< real64, 2, 3*maxNumComps > workSpace( 3, numComps ); + StackArray< real64, 1, maxNumComps+1 > residualSpace( numComps+1 ); arraySlice1d< real64 > logLiquidFugacity = workSpace[0]; arraySlice1d< real64 > logVapourFugacity = workSpace[1]; - arraySlice1d< real64 > fugacityRatios = workSpace[2]; - stackArray1d< integer, maxNumComps > componentIndices( numComps ); + arraySlice1d< real64 > residual = residualSpace.toSlice(); auto const & kVapourLiquid = kValues[0]; + stackArray1d< integer, maxNumComps > availableComponents( numComps ); + stackArray1d< integer, maxNumComps > unavailableComponents( numComps ); - calculatePresentComponents( numComps, composition, componentIndices ); - - auto const presentComponents = componentIndices.toSliceConst(); + calculatePresentComponents( numComps, composition, availableComponents ); + calculateAbsentComponents( numComps, composition, unavailableComponents ); + auto const presentComponents = availableComponents.toSliceConst(); + auto const absentComponents = unavailableComponents.toSliceConst(); // Extract flash parameters integer const maxIterations = discreteFlashParameters[FlashParameters::FLASH_MAX_ITERATIONS]; real64 const flashTolerance = continuousFlashParameters[FlashParameters::FLASH_TOLERANCE]; + real64 const flashSSITolerance = continuousFlashParameters[FlashParameters::SSI_TOLERANCE]; + real64 const negativeFlashTolerance = continuousFlashParameters[FlashParameters::NEGATIVE_FLASH_TOLERANCE]; // Initialise compositions to feed composition for( integer ic = 0; ic < numComps; ++ic ) @@ -72,61 +81,91 @@ bool NegativeTwoPhaseFlash::compute( integer const numComps, vapourComposition[ic] = composition[ic]; } + real64 const & V = vapourPhaseMoleFraction; + bool converged = false; - for( localIndex iterationCount = 0; iterationCount < maxIterations; ++iterationCount ) + localIndex iterationCount = 0; + + // Start SSI iterations + for(; ( !converged ) && (iterationCount < maxIterations); ++iterationCount ) { // Solve Rachford-Rice Equation vapourPhaseMoleFraction = RachfordRice::solve( kVapourLiquid.toSliceConst(), composition, presentComponents ); - // Assign phase compositions - for( integer ic = 0; ic < numComps; ++ic ) - { - liquidComposition[ic] = composition[ic] / ( 1.0 + vapourPhaseMoleFraction * ( kVapourLiquid[ic] - 1.0 ) ); - vapourComposition[ic] = kVapourLiquid[ic] * liquidComposition[ic]; - } - - normalizeComposition( numComps, liquidComposition ); - normalizeComposition( numComps, vapourComposition ); - - FugacityCalculator::computeLogFugacity( numComps, + real64 const error = calculateResidual( numComps, pressure, temperature, - liquidComposition.toSliceConst(), + composition, componentProperties, - flashData.liquidEos, flashData, - logLiquidFugacity ); - FugacityCalculator::computeLogFugacity( numComps, - pressure, - temperature, - vapourComposition.toSliceConst(), - componentProperties, - flashData.vapourEos, - flashData, - logVapourFugacity ); + kVapourLiquid.toSliceConst(), + vapourPhaseMoleFraction, + liquidComposition, + vapourComposition, + logLiquidFugacity, + logVapourFugacity, + residual ); - // Compute fugacity ratios and calculate the error - real64 error = 0.0; + // Update K-values for( integer const ic : presentComponents ) { - fugacityRatios[ic] = ( logLiquidFugacity[ic] - logVapourFugacity[ic] ) + log( liquidComposition[ic] ) - log( vapourComposition[ic] ); - error += (fugacityRatios[ic]*fugacityRatios[ic]); + kVapourLiquid[ic] = exp( logLiquidFugacity[ic] - logVapourFugacity[ic] ); } - error = LvArray::math::sqrt( error ); - // Compute fugacity ratios and check convergence - converged = (error < flashTolerance); + // Check convergence + converged = ( error < flashTolerance ); - if( converged ) + // Don't switch to Newton if we are in the "negative flash" range + // If we have reached the SSI tolerance and are not non-pysical switch to Newton + if( error < flashSSITolerance && -negativeFlashTolerance < V && -negativeFlashTolerance < 1.0-V ) { break; } + } + + // Start Newton iterations + StackArray< real64, 3, 2*maxNumComps * maxNumDofs > logFugacityDerivs( 2, numComps, numDofs ); + StackArray< real64, 2, (maxNumComps+1) * (maxNumComps+1) > jacobian( numComps+1, numComps+1 ); + arraySlice2d< real64 > const logLiquidFugacityDerivs = logFugacityDerivs[0]; + arraySlice2d< real64 > const logVapourFugacityDerivs = logFugacityDerivs[1]; + + for( iterationCount++; ( !converged ) && (iterationCount < maxIterations); ++iterationCount ) + { + real64 const error = calculateResidualAndJacobian( numComps, + pressure, + temperature, + composition, + componentProperties, + flashData, + kVapourLiquid.toSliceConst(), + vapourPhaseMoleFraction, + liquidComposition, + vapourComposition, + logLiquidFugacity, + logVapourFugacity, + logLiquidFugacityDerivs, + logVapourFugacityDerivs, + residual, + jacobian.toSlice() ); + + // Account for missing components in Jacobian + for( integer const ic : absentComponents ) + { + jacobian( ic, ic ) = 1.0; + } + + // Solve for next step + solveLinearSystem( jacobian.toSlice(), residual ); // Update K-values - for( integer ic = 0; ic < numComps; ++ic ) + for( integer const ic : presentComponents ) { - kVapourLiquid[ic] *= exp( fugacityRatios[ic] ); + kVapourLiquid[ic] *= exp( -residual[ic] ); } + vapourPhaseMoleFraction -= residual[numComps]; + + // Check convergence + converged = ( error < flashTolerance ); } // Test if we have converged to a null or trivial solution @@ -150,19 +189,18 @@ bool NegativeTwoPhaseFlash::compute( integer const numComps, template< integer USD1, integer USD2, integer USD3 > GEOS_HOST_DEVICE -void NegativeTwoPhaseFlash::computeDerivatives( - integer const numComps, - real64 const pressure, - real64 const temperature, - arraySlice1d< real64 const > const & composition, - ComponentProperties::KernelWrapper const & componentProperties, - FlashData const & flashData, - real64 const & vapourFraction, - arraySlice1d< real64 const, USD1 > const & liquidComposition, - arraySlice1d< real64 const, USD1 > const & vapourComposition, - arraySlice1d< real64, USD2 > const & vapourFractionDerivs, - arraySlice2d< real64, USD3 > const & liquidCompositionDerivs, - arraySlice2d< real64, USD3 > const & vapourCompositionDerivs ) +void NegativeTwoPhaseFlash::computeDerivatives( integer const numComps, + real64 const pressure, + real64 const temperature, + arraySlice1d< real64 const > const & composition, + ComponentProperties::KernelWrapper const & componentProperties, + FlashData const & flashData, + real64 const & vapourFraction, + arraySlice1d< real64 const, USD1 > const & liquidComposition, + arraySlice1d< real64 const, USD1 > const & vapourComposition, + arraySlice1d< real64, USD2 > const & vapourFractionDerivs, + arraySlice2d< real64, USD3 > const & liquidCompositionDerivs, + arraySlice2d< real64, USD3 > const & vapourCompositionDerivs ) { constexpr integer maxNumComps = MultiFluidConstants::MAX_NUM_COMPONENTS; constexpr integer maxNumDofs = MultiFluidConstants::MAX_NUM_COMPONENTS + 2; @@ -261,21 +299,248 @@ void NegativeTwoPhaseFlash::computeDerivatives( } } +template< int USD1, int USD2, int USD3 > +GEOS_HOST_DEVICE +real64 NegativeTwoPhaseFlash::calculateResidual( integer const numComps, + real64 const pressure, + real64 const temperature, + arraySlice1d< real64 const > const & composition, + ComponentProperties::KernelWrapper const & componentProperties, + FlashData const & flashData, + arraySlice1d< real64 const, USD1 > const & kValues, + real64 const & vapourPhaseMoleFraction, + arraySlice1d< real64, USD2 > const & liquidComposition, + arraySlice1d< real64, USD2 > const & vapourComposition, + arraySlice1d< real64 > const & logLiquidFugacity, + arraySlice1d< real64 > const & logVapourFugacity, + arraySlice1d< real64, USD3 > const & residual ) +{ + // Calculate phase and normalize compositions + real64 massBalance = 0.0; + for( integer ic = 0; ic < numComps; ++ic ) + { + liquidComposition[ic] = composition[ic] / ( 1.0 + vapourPhaseMoleFraction * ( kValues[ic] - 1.0 ) ); + vapourComposition[ic] = kValues[ic] * liquidComposition[ic]; + massBalance += (vapourComposition[ic] - liquidComposition[ic]); + } + + normalizeComposition( numComps, liquidComposition ); + normalizeComposition( numComps, vapourComposition ); + + FugacityCalculator::computeLogFugacity( numComps, + pressure, + temperature, + liquidComposition.toSliceConst(), + componentProperties, + flashData.liquidEos, + flashData, + logLiquidFugacity ); + FugacityCalculator::computeLogFugacity( numComps, + pressure, + temperature, + vapourComposition.toSliceConst(), + componentProperties, + flashData.vapourEos, + flashData, + logVapourFugacity ); + + // Compute fugacity ratios and calculate the error + real64 error = 0.0; + for( integer ic = 0; ic < numComps; ++ic ) + { + residual[ic] = liquidComposition[ic]*LvArray::math::exp( logLiquidFugacity[ic] ) - vapourComposition[ic]*LvArray::math::exp( logVapourFugacity[ic] ); + error += (residual[ic]*residual[ic]); + } + residual[numComps] = massBalance; + error += (residual[numComps]*residual[numComps]); + error = LvArray::math::sqrt( error ); + + return error; +} + +template< int USD1, int USD2, int USD3, int USD4 > +GEOS_HOST_DEVICE +real64 NegativeTwoPhaseFlash::calculateResidualAndJacobian( integer const numComps, + real64 const pressure, + real64 const temperature, + arraySlice1d< real64 const > const & composition, + ComponentProperties::KernelWrapper const & componentProperties, + FlashData const & flashData, + arraySlice1d< real64 const, USD1 > const & kValues, + real64 const & vapourPhaseMoleFraction, + arraySlice1d< real64, USD2 > const & liquidComposition, + arraySlice1d< real64, USD2 > const & vapourComposition, + arraySlice1d< real64 > const & logLiquidFugacity, + arraySlice1d< real64 > const & logVapourFugacity, + arraySlice2d< real64 > const & logLiquidFugacityDerivs, + arraySlice2d< real64 > const & logVapourFugacityDerivs, + arraySlice1d< real64, USD3 > const & residual, + arraySlice2d< real64, USD4 > const & jacobian ) +{ + real64 const & V = vapourPhaseMoleFraction; + + // Calculate phase compositions --- + for( integer ic = 0; ic < numComps; ++ic ) + { + real64 const denominator = 1.0 + V * (kValues[ic] - 1.0); + liquidComposition[ic] = composition[ic] / denominator; + vapourComposition[ic] = kValues[ic] * liquidComposition[ic]; + } + + real64 const sumX = normalizeComposition( numComps, liquidComposition ); + real64 const sumY = normalizeComposition( numComps, vapourComposition ); + + // Calculate the fugacity and derivatives + FugacityCalculator::computeLogFugacityDerivatives( numComps, + pressure, + temperature, + liquidComposition.toSliceConst(), + componentProperties, + flashData.liquidEos, + flashData, + logLiquidFugacity, + logLiquidFugacityDerivs ); + FugacityCalculator::computeLogFugacityDerivatives( numComps, + pressure, + temperature, + vapourComposition.toSliceConst(), + componentProperties, + flashData.vapourEos, + flashData, + logVapourFugacity, + logVapourFugacityDerivs ); + + // Convert from log (use same memory space: log fugacity no longer required) + arraySlice1d< real64 > const & liquidFugacity = logLiquidFugacity; + arraySlice1d< real64 > const & vapourFugacity = logVapourFugacity; + for( integer ic = 0; ic < numComps; ++ic ) + { + liquidFugacity[ic] = LvArray::math::exp( logLiquidFugacity[ic] ); + vapourFugacity[ic] = LvArray::math::exp( logVapourFugacity[ic] ); + } + + // Calculate the residual vector + real64 error = 0.0; + residual[numComps] = 0.0; + for( integer ic = 0; ic < numComps; ++ic ) + { + real64 const zi = composition[ic]; + real64 const xi = liquidComposition[ic]; + real64 const yi = vapourComposition[ic]; + real64 const Ki = kValues[ic]; + residual[ic] = xi * liquidFugacity[ic] - yi * vapourFugacity[ic]; + residual[numComps] += zi * (Ki - 1.0) / (1.0 + V * (Ki - 1.0)); + error += (residual[ic]*residual[ic]); + } + error += (residual[numComps]*residual[numComps]); + error = LvArray::math::sqrt( error ); + + real64 const invSumX = 1.0/sumX; + real64 const invSumY = 1.0/sumY; + + real64 sum_d_xs_d_V = 0.0; + real64 sum_d_ys_d_V = 0.0; + for( int ic = 0; ic < numComps; ++ic ) + { + real64 const zi = composition[ic]; + real64 const Ki = kValues[ic]; + + real64 const invKi = 1.0 / ( 1.0 + V * (Ki - 1.0) ); + real64 const d_xsi_d_V = -zi * (Ki - 1.0) * invKi * invKi; + real64 const d_ysi_d_V = Ki * d_xsi_d_V; + + sum_d_xs_d_V += d_xsi_d_V; + sum_d_ys_d_V += d_ysi_d_V; + } + + // Start calculating the Jacobian + real64 sum_j_N_V = 0.0; + for( integer ic = 0; ic < numComps; ++ic ) + { + real64 const zi = composition[ic]; + real64 const xi = liquidComposition[ic]; + real64 const yi = vapourComposition[ic]; + real64 const Ki = kValues[ic]; + + real64 const invKi = 1.0 / ( 1.0 + V * (Ki - 1.0) ); + real64 const d_xsi_d_V = -zi * (Ki - 1.0) * invKi * invKi; + real64 const d_ysi_d_V = Ki * d_xsi_d_V; + + real64 const d_xi_d_V = (d_xsi_d_V - xi * sum_d_xs_d_V) * invSumX; + real64 const d_yi_d_V = (d_ysi_d_V - yi * sum_d_ys_d_V) * invSumY; + + real64 sum_dXid_xj_dxj_d_V = 0.0; + real64 sum_dYid_yj_dyj_d_V = 0.0; + + for( integer jc = 0; jc < numComps; ++jc ) + { + real64 const zj = composition[jc]; + real64 const xj = liquidComposition[jc]; + real64 const yj = vapourComposition[jc]; + real64 const Kj = kValues[jc]; + + real64 const dphiLi_dxj = logLiquidFugacityDerivs( ic, Deriv::dC + jc ); + real64 const dphiVi_dyj = logVapourFugacityDerivs( ic, Deriv::dC + jc ); + + real64 const invKj = 1.0 / ( 1.0 + V * (Kj - 1.0) ); + real64 const d_xsj_d_V = -zj * (Kj - 1.0) * invKj * invKj; + real64 const d_ysj_d_V = Kj * d_xsj_d_V; + + real64 const d_xsj_d_kj = -zj * V * Kj * invKj * invKj; + real64 const d_ysj_d_kj = Kj * (xj * sumX + d_xsj_d_kj); + + real64 const d_xj_d_V = (d_xsj_d_V - xj * sum_d_xs_d_V) * invSumX; + real64 const d_yj_d_V = (d_ysj_d_V - yj * sum_d_ys_d_V) * invSumY; + + sum_dXid_xj_dxj_d_V += dphiLi_dxj * d_xj_d_V; + sum_dYid_yj_dyj_d_V += dphiVi_dyj * d_yj_d_V; + + real64 const dij = (ic == jc) ? 1.0 : 0.0; + real64 const d_xi_d_kj = d_xsj_d_kj * (dij - xi) * invSumX; + real64 const d_yi_d_kj = d_ysj_d_kj * (dij - yi) * invSumY; + + real64 sum_dXid_xk_dxk_d_kj = 0.0; + real64 sum_dYid_yk_dyk_d_kj = 0.0; + for( integer kc = 0; kc < numComps; ++kc ) + { + real64 const djk = (jc == kc) ? 1.0 : 0.0; + real64 const d_xk_d_kj = d_xsj_d_kj * (djk - liquidComposition[kc]) * invSumX; + real64 const d_yk_d_kj = d_ysj_d_kj * (djk - vapourComposition[kc]) * invSumY; + + real64 const dphiLi_dxk = logLiquidFugacityDerivs( ic, Deriv::dC + kc ); + real64 const dphiVi_dyk = logVapourFugacityDerivs( ic, Deriv::dC + kc ); + + sum_dXid_xk_dxk_d_kj += dphiLi_dxk * d_xk_d_kj; + sum_dYid_yk_dyk_d_kj += dphiVi_dyk * d_yk_d_kj; + } + jacobian( ic, jc ) = liquidFugacity[ic] * (d_xi_d_kj + xi * sum_dXid_xk_dxk_d_kj) - + vapourFugacity[ic] * (d_yi_d_kj + yi * sum_dYid_yk_dyk_d_kj); + } + jacobian( ic, numComps ) = liquidFugacity[ic] * (d_xi_d_V + xi * sum_dXid_xj_dxj_d_V) - + vapourFugacity[ic] * (d_yi_d_V + yi * sum_dYid_yj_dyj_d_V); + jacobian( numComps, ic ) = zi * Ki * invKi* invKi; + + sum_j_N_V -= zi * (Ki - 1.0) * (Ki - 1.0) * invKi * invKi; + } + jacobian( numComps, numComps ) = sum_j_N_V; + + return error; +} + template< integer USD1, integer USD2, integer USD3 > GEOS_HOST_DEVICE -void NegativeTwoPhaseFlash::computeDerivatives( - integer const numComps, - arraySlice1d< real64 const > const & totalComposition, - real64 const & phase1Fraction, - arraySlice1d< real64 const, USD1 > const & phase1Composition, - arraySlice1d< real64 const, USD1 > const & phase2Composition, - arraySlice1d< real64 const > const & phase1Fugacity, - arraySlice1d< real64 const > const & phase2Fugacity, - arraySlice2d< real64 const > const & phase1LogFugacityDerivs, - arraySlice2d< real64 const > const & phase2LogFugacityDerivs, - arraySlice1d< real64, USD2 > const & phase1FractionDerivs, - arraySlice2d< real64, USD3 > const & phase1CompositionDerivs, - arraySlice2d< real64, USD3 > const & phase2CompositionDerivs ) +void NegativeTwoPhaseFlash::computeDerivatives( integer const numComps, + arraySlice1d< real64 const > const & totalComposition, + real64 const & phase1Fraction, + arraySlice1d< real64 const, USD1 > const & phase1Composition, + arraySlice1d< real64 const, USD1 > const & phase2Composition, + arraySlice1d< real64 const > const & phase1Fugacity, + arraySlice1d< real64 const > const & phase2Fugacity, + arraySlice2d< real64 const > const & phase1LogFugacityDerivs, + arraySlice2d< real64 const > const & phase2LogFugacityDerivs, + arraySlice1d< real64, USD2 > const & phase1FractionDerivs, + arraySlice2d< real64, USD3 > const & phase1CompositionDerivs, + arraySlice2d< real64, USD3 > const & phase2CompositionDerivs ) { constexpr integer maxNumDofs = MultiFluidConstants::MAX_NUM_COMPONENTS + 2; constexpr integer maxNumVals = MultiFluidConstants::MAX_NUM_COMPONENTS+1; @@ -302,7 +567,7 @@ void NegativeTwoPhaseFlash::computeDerivatives( real64 const xj = phase2Composition[j]; real64 const yj = phase1Composition[j]; - real64 const c1_ij = xi * phi_2_i * phase2LogFugacityDerivs( i, Deriv::dC+j ); + real64 const c1_ij = xi * phi_2_i * phase2LogFugacityDerivs( i, Deriv::dC+j ); real64 const c2_ij = -yi * phi_1_i * phase1LogFugacityDerivs( i, Deriv::dC+j ); A( i, j ) = c2_ij - factor1 * c1_ij; diff --git a/src/coreComponents/constitutive/fluid/multifluid/compositional/functions/StabilityTest.hpp b/src/coreComponents/constitutive/fluid/multifluid/compositional/functions/StabilityTest.hpp index 3e017c005ea..8fef4857509 100644 --- a/src/coreComponents/constitutive/fluid/multifluid/compositional/functions/StabilityTest.hpp +++ b/src/coreComponents/constitutive/fluid/multifluid/compositional/functions/StabilityTest.hpp @@ -38,29 +38,70 @@ namespace constitutive namespace compositional { +/** + * @struct StabilityTest + * @brief Performs stability analysis for a reservoir fluid. + * + * This class implements a stability test to determine if a single-phase fluid mixture is thermodynamically + * stable or if it will split into multiple phases. The algorithm follows the tangent-plane-distance (TPD) + * method, as described in the reference by Michelsen (1982). + * + * The method involves finding a trial composition ($y_i$) that minimizes the tangent-plane-distance ($TPD$) + * at a given pressure and temperature. The TPD for a trial phase with composition $y$ relative to a sample + * phase with composition $z$ is given by: + * $$ TPD(y) = \sum_{i=1}^{nc} y_i \left[ \ln(y_i) - \ln(z_i) + \ln(\phi}_i(y)) - \ln(\phi_i(z)) \right] $$ + * where $nc$ is the number of components, $z_i$ and $y_i$ are the mole fractions of component $i$ in + * the sample and trial phases, respectively, and $\phi_i$ and $\hat{\phi}_i$ are the fugacity coefficients + * of component $i$ in the sample and trial phases, respectively. + * + * The stability analysis is performed by solving for the stationary points of the TPD function. The solution + * is found by a combination of successive substitution iterations (SSI) and Newton-Raphson iterations. + * SSI is used initially to provide a good starting point for the faster, but more sensitive, Newton-Raphson + * method. The mixture is considered unstable if the minimum TPD is less than zero, indicating that a more + * stable phase exists. + * + * Reference: Michelsen, M.L. (1982). The Isothermal Flash Problem. Part I. Stability. + * Fluid Phase Equilibria, 9(1), 1–19. https://doi.org/10.1016/0378-3812(82)85001-2 + */ struct StabilityTest { private: static constexpr integer maxNumComps = MultiFluidConstants::MAX_NUM_COMPONENTS; + static constexpr integer maxNumDofs = MultiFluidConstants::MAX_NUM_COMPONENTS + 2; using Deriv = constitutive::multifluid::DerivativeOffset; + public: /** - * @brief Perform a two-phase stability test - * @param[in] numComps number of components - * @param[in] pressure pressure - * @param[in] temperature temperature - * @param[in] composition composition of the mixture - * @param[in] componentProperties The compositional component properties - * @param[in] flashData The parameters required for the flash - * @param[in] kValues the k-values to use to create test samples - * @param[in] continuousFlashParameters List of continuous (float) parameters for flash - * @param[in] discreteFlashParameters List of discrete (integer) parameters for flash - * @param[out] unstableMixture The stability test result indicating if the mixture is unstable - * @param[out] incipientEquationOfState The equation of state of the incipient phase - * @param[out] incipientComposition The composition of the incipient phase. This will be - * set to the incipient phase composition at the stationary point that is furthest - * from the trivial solution. - * @return a flag indicating that the stability test is successful i.e. mixture is unstable or all test points converged + * @brief Perform a two-phase stability test using the Tangent Plane Distance (TPD) method. + * + * The stability analysis is performed by solving for the stationary points of the TPD function. The solution + * is found by a combination of successive substitution iterations (SSI) and Newton-Raphson iterations. + * SSI is used initially to provide a good starting point for the faster, but more sensitive, Newton-Raphson + * method. The mixture is considered unstable if the minimum TPD is less than zero, indicating that a more + * stable phase exists. + * + * This implementation supports asymmetric equations of state, where the EOS used for the + * incipient phase may differ from that of the original mixture. This is handled by testing + * both configurations independently. + * + * @tparam USD1 Unit Stride Dimension for the composition array. + * @tparam USD2 Unit Stride Dimension for the K-values array. + * + * @param[in] numComps Number of components in the mixture. + * @param[in] pressure Pressure of the system [Pa]. + * @param[in] temperature Temperature of the system [K]. + * @param[in] composition Mole fractions of each component in the mixture (must be normalized). + * @param[in] componentProperties Thermodynamic properties of each component. + * @param[in] flashData Flash calculation context, including EOS types for liquid and vapor phases. + * @param[in] kValues Initial equilibrium ratios (K-values) used to generate trial compositions. + * @param[in] continuousFlashParameters Continuous parameters (e.g., tolerances, thresholds) for the flash calculation. + * @param[in] discreteFlashParameters Discrete parameters (e.g., max iterations) for the flash calculation. + * @param[out] unstableMixture Flag indicating whether the mixture is unstable. + * @param[out] incipientEquationOfState EOS type of the incipient phase (liquid or vapor). + * @param[out] incipientComposition Composition of the incipient phase at the stationary point furthest from the trivial solution. + * + * @return True if the mixture is unstable or all trial compositions converge to stationary points; false otherwise. + * */ template< integer USD1, integer USD2 > GEOS_HOST_DEVICE @@ -75,173 +116,34 @@ struct StabilityTest arraySlice1d< integer const > const & discreteFlashParameters, bool & unstableMixture, EquationOfStateType & incipientEquationOfState, - arraySlice1d< real64 > const & incipientComposition ) - { - stackArray2d< real64, 4*maxNumComps > workSpace( 4, numComps ); - arraySlice1d< real64 > logFugacity = workSpace[0]; - arraySlice1d< real64 > normalizedComposition = workSpace[1]; - arraySlice1d< real64 > logTrialComposition = workSpace[2]; - arraySlice1d< real64 > hyperplane = workSpace[3]; // h-parameter - stackArray1d< integer, maxNumComps > availableComponents( numComps ); - - //integer const presentCount = - calculatePresentComponents( numComps, composition, availableComponents ); - auto const presentComponents = availableComponents.toSliceConst(); - - LvArray::forValuesInSlice( workSpace.toSlice(), setZero ); - - // Extract flash parameters - integer const maxIterations = discreteFlashParameters[FlashParameters::STABILITY_MAX_ITERATIONS]; - real64 const stabilityThreshold = continuousFlashParameters[FlashParameters::STABILITY_THRESHOLD]; - real64 const stabilityTolerance = continuousFlashParameters[FlashParameters::STABILITY_TOLERANCE]; - - integer const numConfigSteps = flashData.liquidEos == flashData.vapourEos ? 1 : 2; - - // Initialise by assuming a stable mixture - unstableMixture = false; - - // Flag to indicate all trial compositions converged - bool allConverged = true; - - // Measure of distance to trivial solution - real64 maxDistanceToTrivialSolution = 0.0; - - for( integer configStep = 0; configStep < numConfigSteps; ++configStep ) - { - EquationOfStateType const sampleEos = (configStep == 0) ? flashData.liquidEos : flashData.vapourEos; - EquationOfStateType const incipientEos = (configStep == 0) ? flashData.vapourEos : flashData.liquidEos; - - // Calculate the hyperplane parameter - // h_i = log( z_i ) + log( phi_i ) - FugacityCalculator::computeLogFugacity( numComps, - pressure, - temperature, - composition, - componentProperties, - sampleEos, - flashData, - logFugacity ); - for( integer const ic : presentComponents ) - { - hyperplane[ic] = LvArray::math::log( composition[ic] ) + logFugacity[ic]; - } - - for( real64 const alpha : { 1.0, -1.0, 0.0 } ) - { - // Initialise next sample - if( LvArray::math::abs( alpha ) < MultiFluidConstants::epsilon ) - { - for( integer const ic : presentComponents ) - { - logTrialComposition[ic] = 0.0; - normalizedComposition[ic] = 1.0; - } - } - else - { - for( integer const ic : presentComponents ) - { - logTrialComposition[ic] = LvArray::math::log( composition[ic] ) + alpha*LvArray::math::log( kValues[ic] ); - normalizedComposition[ic] = LvArray::math::exp( logTrialComposition[ic] ); - } - } - - // Start iterations for this sample - bool converged = false; - - for( localIndex iterationCount = 0; iterationCount < maxIterations; ++iterationCount ) - { - // Normalise the composition and calculate the fugacity - normalizeComposition( numComps, normalizedComposition ); - FugacityCalculator::computeLogFugacity( numComps, - pressure, - temperature, - normalizedComposition.toSliceConst(), - componentProperties, - incipientEos, - flashData, - logFugacity ); - - // Calculate the stationarity condition - real64 error = 0.0; - for( integer const ic : presentComponents ) - { - real64 const dG = logTrialComposition[ic] + logFugacity[ic] - hyperplane[ic]; - error += (dG*dG); - } - error = LvArray::math::sqrt( error ); - - // Check stationarity - if( error < stabilityTolerance ) - { - converged = true; - break; - } - - // Update to next step - for( integer const ic : presentComponents ) - { - logTrialComposition[ic] = hyperplane[ic] - logFugacity[ic]; - normalizedComposition[ic] = LvArray::math::exp( logTrialComposition[ic] ); - } - } - allConverged = allConverged && converged; - - // Calculate the tangent-plane-distance (TPD) and distance to the trivial solution - real64 tpd = 1.0; - real64 distance = 0.0; - for( integer const ic : presentComponents ) - { - real64 const dZ = normalizedComposition[ic] - composition[ic]; - real64 const dG = logTrialComposition[ic] + logFugacity[ic] - hyperplane[ic]; - tpd += normalizedComposition[ic] * (dG - 1.0); - distance += (dZ*dZ); - } - - if( maxDistanceToTrivialSolution < distance ) - { - maxDistanceToTrivialSolution = distance; - incipientEquationOfState = incipientEos; - for( integer ic = 0; ic < numComps; ++ic ) - { - incipientComposition[ic] = normalizedComposition[ic]; - } - } - // The mixture is unstable if either the TPD is negative or if the TPD is zero but the incipient composition - // is different from the sample composition - // stabilityThreshold is negative (default -1e-8) - if( tpd < stabilityThreshold ) - { - unstableMixture = true; - } - else if((tpd < -stabilityThreshold) && (-stabilityThreshold < distance)) - { - unstableMixture = true; - } - if( unstableMixture ) - { - return true; - } - } - } - - // The test is successful if either we have an unstable mixture or all test compositions converged to stationarity - return unstableMixture || allConverged; - } + arraySlice1d< real64 > const & incipientComposition ); /** - * @brief Compute derivatives of the incipient phase composition - * @param[in] numComps number of components - * @param[in] pressure pressure - * @param[in] temperature temperature - * @param[in] composition composition of the mixture - * @param[in] componentProperties The compositional component properties - * @param[in] sampleEquationOfState The equation of state for the sample - * @param[in] sampleEquationOfState The equation of state for the incipient phase - * @param[in] flashData The parameters required for the flash - * @param[in] incipientComposition The composition of the incipient phase - * @param[out] incipientCompositionDerivs Derivatives of the composition of the incipient phase - * @param[out] compositionDerivs Workspace + * @brief Compute derivatives of the incipient phase composition with respect to pressure, temperature, and overall composition. + * + * @details This method calculates the sensitivity of the incipient phase composition resulting from a phase stability test. + * It computes the partial derivatives of the incipient composition with respect to pressure, temperature, and each component + * of the overall mixture composition. + * + * The implementation solves a linear system derived from the stationarity condition of the Tangent Plane Distance (TPD) + * function. It uses fugacity derivatives from both the sample and incipient phases, and enforces the constraint that the + * incipient composition must sum to unity. + * + * @tparam USD1 Unit Stride Dimension for the composition array. + * @tparam USD2 Unit Stride Dimension for the derivative arrays. + * + * @param[in] numComps Number of components in the mixture. + * @param[in] pressure Pressure of the system [Pa]. + * @param[in] temperature Temperature of the system [K]. + * @param[in] composition Mole fractions of each component in the original mixture. + * @param[in] componentProperties Thermodynamic properties of each component. + * @param[in] sampleEquationOfState EOS used for the original mixture. + * @param[in] incipientEquationOfState EOS used for the incipient phase. + * @param[in] flashData Flash calculation context and EOS configuration. + * @param[in] incipientComposition Composition of the incipient phase. + * @param[out] incipientCompositionDerivs Derivatives of the incipient composition with respect to pressure, temperature, and composition. + * @param[out] compositionDerivs Workspace for storing derivatives of the sample composition. + * */ template< integer USD1, integer USD2 > GEOS_HOST_DEVICE @@ -255,86 +157,101 @@ struct StabilityTest FlashData const & flashData, arraySlice1d< real64 const > const & incipientComposition, arraySlice2d< real64, USD2 > const & incipientCompositionDerivs, - arraySlice2d< real64, USD2 > const & compositionDerivs ) - { - integer constexpr maxNumRows = MultiFluidConstants::MAX_NUM_COMPONENTS + 1; - integer constexpr maxDofs = MultiFluidConstants::MAX_NUM_COMPONENTS + 2; - - integer const numDofs = 2 + numComps; - - StackArray< real64, 1, maxNumComps > logFugacity( numComps ); - auto const & logFugacityIncipientDerivs = incipientCompositionDerivs; - auto const & logFugacitySampleDerivs = compositionDerivs; - - FugacityCalculator::computeLogFugacityDerivatives( numComps, - pressure, - temperature, - composition, - componentProperties, - sampleEquationOfState, - flashData, - logFugacity.toSlice(), - logFugacitySampleDerivs ); - FugacityCalculator::computeLogFugacityDerivatives( numComps, - pressure, - temperature, - incipientComposition, - componentProperties, - incipientEquationOfState, - flashData, - logFugacity.toSlice(), - logFugacityIncipientDerivs ); - - StackArray< real64, 2, maxNumRows * maxNumRows, MatrixLayout::COL_MAJOR_PERM > A( numComps + 1, numComps + 1 ); - StackArray< real64, 2, maxNumRows * maxDofs, MatrixLayout::COL_MAJOR_PERM > X( numComps + 1, numDofs ); + arraySlice2d< real64, USD2 > const & compositionDerivs ); - for( integer ic = 0; ic < numComps; ++ic ) - { - real64 const yi = incipientComposition[ic]; - for( integer jc = 0; jc < numComps; ++jc ) - { - real64 const df_dyj = logFugacityIncipientDerivs( ic, Deriv::dC+jc ); - A( ic, jc ) = yi*df_dyj; - } - A( ic, ic ) += 1.0; - A( ic, numComps ) = yi; - A( numComps, ic ) = 1.0; - } - A( numComps, numComps ) = 0.0; - - for( integer ic = 0; ic < numComps; ++ic ) - { - real64 const yi = incipientComposition[ic]; - for( integer const idof : {Deriv::dP, Deriv::dT} ) - { - X( ic, idof ) = yi*(-logFugacityIncipientDerivs( ic, idof ) + logFugacitySampleDerivs( ic, idof )); - } - for( integer jc = 0; jc < numComps; ++jc ) - { - integer const idof = Deriv::dC + jc; - X( ic, idof ) = yi*logFugacitySampleDerivs( ic, idof ); - } - if( MultiFluidConstants::epsilon < composition[ic] ) - { - X( ic, Deriv::dC+ic ) += yi/composition[ic]; - } - } - for( integer idof = 0; idof < numDofs; ++idof ) - { - X( numComps, idof ) = 0.0; - } - - // Solve linear system - solveLinearSystem( A.toSlice(), X.toSlice() ); + /** + * @brief Calculates the residual vector and Tangent Plane Distance (TPD) for a trial composition. + * + * @details This method evaluates the deviation of a trial composition from stationarity + * in the context of phase stability analysis. It computes the residual vector used in + * iterative solvers (SSI/Newton) and the tangent plane distance (TPD) which indicates + * whether the trial composition has lower Gibbs energy than the original mixture. + * + * The trial composition is provided in logarithmic form and exponentiated internally. + * Missing components (not present in presentComponents) are ignored in the residual + * and TPD calculation. The composition is normalized after exponentiation. + * + * @tparam USD1 Unit Stride Dimension for logTestComposition array. + * @tparam USD2 Unit Stride Dimension for testComposition array. + * @tparam USD3 Unit Stride Dimension for residual array. + * + * @param[in] numComps Number of components in the mixture. + * @param[in] presentComponents Indices of components present in the mixture. + * @param[in] pressure System pressure [Pa]. + * @param[in] temperature System temperature [K]. + * @param[in] componentProperties Thermodynamic properties of each component. + * @param[in] equationOfState Equation of state used for fugacity calculations. + * @param[in] flashData Flash calculation context. + * @param[in] logTestComposition Logarithm of trial composition. + * @param[out] tangentPlaneDistance Computed tangent plane distance (TPD). + * @param[out] testComposition Trial composition (exponentiated and normalized). + * @param[out] testLogFugacity Log fugacity coefficients of the trial composition. + * @param[out] residual Residual vector for iterative solver. + * + * @return Error norm of the residual vector. + * + */ + template< int USD1, int USD2, int USD3 > + GEOS_HOST_DEVICE + static real64 calculateResidual( integer const numComps, + arraySlice1d< integer const > const & presentComponents, + real64 const pressure, + real64 const temperature, + ComponentProperties::KernelWrapper const & componentProperties, + EquationOfStateType const equationOfState, + FlashData const & flashData, + arraySlice1d< real64 const, USD1 > const & logTestComposition, + real64 & tangentPlaneDistance, + arraySlice1d< real64, USD2 > const & testComposition, + arraySlice1d< real64 > const & testLogFugacity, + arraySlice1d< real64, USD3 > const & residual ); - for( integer idof = 0; idof < numDofs; ++idof ) - { - for( integer ic = 0; ic < numComps; ++ic ) - { - incipientCompositionDerivs( ic, idof ) = X( ic, idof ); - } - } - } + /** + * @brief Calculates the residual vector and Jacobian matrix for a trial composition in phase stability analysis. + * + * @details This method evaluates the stationarity of a trial composition using the Tangent Plane Distance (TPD) criterion. + * It computes the residual vector and the Jacobian matrix required for Newton-Raphson iterations. + * The fugacity and its derivatives are calculated using the specified equation of state. + * Missing components (not present in the trial composition) are ignored in the residual and Jacobian calculations. + * + * @tparam USD1 Unit Stride Dimension for logTestComposition array. + * @tparam USD2 Unit Stride Dimension for testComposition array. + * @tparam USD3 Unit Stride Dimension for residual array. + * @tparam USD4 Unit Stride Dimension for jacobian array. + * + * @param[in] numComps Number of components in the mixture. + * @param[in] presentComponents Indices of components present in the trial composition. + * @param[in] pressure Pressure of the system [Pa]. + * @param[in] temperature Temperature of the system [K]. + * @param[in] componentProperties Thermodynamic properties of each component. + * @param[in] equationOfState Equation of state used for fugacity calculations. + * @param[in] flashData Flash calculation context and parameters. + * @param[in] logTestComposition Logarithm of the trial composition. + * @param[out] tangentPlaneDistance Computed TPD value for the trial composition. + * @param[out] testComposition Normalized trial composition. + * @param[out] testFugacity Log fugacity coefficients of the trial composition. + * @param[out] logTestFugacityDerivs Derivatives of log fugacity with respect to composition. + * @param[out] residual Residual vector used in Newton-Raphson iterations. + * @param[out] jacobian Jacobian matrix of the residual vector. + * + * @return Error norm of the residual vector. + */ + template< int USD1, int USD2, int USD3, int USD4 > + GEOS_HOST_DEVICE + static real64 calculateResidualAndJacobian( integer const numComps, + arraySlice1d< integer const > const & presentComponents, + real64 const pressure, + real64 const temperature, + ComponentProperties::KernelWrapper const & componentProperties, + EquationOfStateType const equationOfState, + FlashData const & flashData, + arraySlice1d< real64 const, USD1 > const & logTestComposition, + real64 & tangentPlaneDistance, + arraySlice1d< real64, USD2 > const & testComposition, + arraySlice1d< real64 > const & testFugacity, + arraySlice2d< real64 > const & logTestFugacityDerivs, + arraySlice1d< real64, USD3 > const & residual, + arraySlice2d< real64, USD4 > const & jacobian ); }; } // namespace compositional @@ -343,4 +260,7 @@ struct StabilityTest } // namespace geos +// Include the implementation +#include "StabilityTest_impl.hpp" + #endif //GEOS_CONSTITUTIVE_FLUID_MULTIFLUID_COMPOSITIONAL_FUNCTIONS_STABILITYTEST_HPP_ diff --git a/src/coreComponents/constitutive/fluid/multifluid/compositional/functions/StabilityTest_impl.hpp b/src/coreComponents/constitutive/fluid/multifluid/compositional/functions/StabilityTest_impl.hpp new file mode 100644 index 00000000000..f74bf07d254 --- /dev/null +++ b/src/coreComponents/constitutive/fluid/multifluid/compositional/functions/StabilityTest_impl.hpp @@ -0,0 +1,494 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file StabilityTest_impl.hpp + */ + +#ifndef GEOS_CONSTITUTIVE_FLUID_MULTIFLUID_COMPOSITIONAL_FUNCTIONS_STABILITYTEST_IMPL_HPP_ +#define GEOS_CONSTITUTIVE_FLUID_MULTIFLUID_COMPOSITIONAL_FUNCTIONS_STABILITYTEST_IMPL_HPP_ + +#include "StabilityTest.hpp" + +namespace geos +{ + +namespace constitutive +{ + +namespace compositional +{ + +template< integer USD1, integer USD2 > +GEOS_HOST_DEVICE +bool StabilityTest::compute( integer const numComps, + real64 const pressure, + real64 const temperature, + arraySlice1d< real64 const, USD1 > const & composition, + ComponentProperties::KernelWrapper const & componentProperties, + FlashData const & flashData, + arraySlice1d< real64 const, USD2 > const & kValues, + arraySlice1d< real64 const > const & continuousFlashParameters, + arraySlice1d< integer const > const & discreteFlashParameters, + bool & unstableMixture, + EquationOfStateType & incipientEquationOfState, + arraySlice1d< real64 > const & incipientComposition ) +{ + integer const numDofs = 2 + numComps; + + stackArray2d< real64, 5*maxNumComps > workSpace( 5, numComps ); + stackArray2d< real64, maxNumComps *maxNumComps > jacobian( numComps, numComps ); + stackArray2d< real64, maxNumComps *maxNumDofs > derivatives( numComps, numDofs ); + arraySlice1d< real64 > logFugacity = workSpace[0]; + arraySlice1d< real64 > logTestComposition = workSpace[1]; + arraySlice1d< real64 > testComposition = workSpace[2]; + arraySlice1d< real64 > residual = workSpace[3]; + arraySlice1d< real64 > hyperplane = workSpace[4]; // h-parameter + stackArray1d< integer, maxNumComps > availableComponents( numComps ); + stackArray1d< integer, maxNumComps > unavailableComponents( numComps ); + + calculatePresentComponents( numComps, composition, availableComponents ); + calculateAbsentComponents( numComps, composition, unavailableComponents ); + auto const presentComponents = availableComponents.toSliceConst(); + auto const absentComponents = unavailableComponents.toSliceConst(); + + LvArray::forValuesInSlice( workSpace.toSlice(), setZero ); + + // Extract flash parameters + integer const maxIterations = discreteFlashParameters[FlashParameters::STABILITY_MAX_ITERATIONS]; + real64 const stabilityThreshold = continuousFlashParameters[FlashParameters::STABILITY_THRESHOLD]; + real64 const stabilityTolerance = continuousFlashParameters[FlashParameters::STABILITY_TOLERANCE]; + real64 const stabilitySSITolerance = continuousFlashParameters[FlashParameters::SSI_TOLERANCE]; + + // Initialise by assuming a stable mixture + unstableMixture = false; + + // Flag to indicate all trial compositions converged + bool allConverged = true; + + // Measure of distance to trivial solution + real64 maxDistanceToTrivialSolution = -LvArray::NumericLimits< real64 >::max; + auto const storeIncipientComposition = [&]( EquationOfStateType const eos ) + { + incipientEquationOfState = eos; + for( integer ic = 0; ic < numComps; ++ic ) + { + incipientComposition[ic] = testComposition[ic]; + } + }; + + // If the 2 EOS are not the same, we will do 2 passes swapping the 2 EOS if required + integer const numConfigSteps = flashData.liquidEos == flashData.vapourEos ? 1 : 2; + + for( integer configStep = 0; configStep < numConfigSteps; ++configStep ) + { + EquationOfStateType const sampleEos = (configStep == 0) ? flashData.liquidEos : flashData.vapourEos; + EquationOfStateType const incipientEos = (configStep == 0) ? flashData.vapourEos : flashData.liquidEos; + + real64 tpd = 0.0; + + // Calculate the hyperplane parameter + // h_i = log( z_i ) + log( phi_i ) + LvArray::forValuesInSlice( hyperplane, setZero ); + for( integer const ic : presentComponents ) + { + logTestComposition[ic] = LvArray::math::log( composition[ic] ); + } + calculateResidual( numComps, + presentComponents, + pressure, + temperature, + componentProperties, + sampleEos, + flashData, + logTestComposition.toSliceConst(), + tpd, + testComposition, + logFugacity, + hyperplane ); + + for( int const alpha : { 1, -1, 0 } ) + { + // If alpha is zero then it means we want to test a constant composition + if( alpha == 0 ) + { + for( integer const ic : presentComponents ) + { + logTestComposition[ic] = 0.0; + } + } + else + { + for( integer const ic : presentComponents ) + { + real64 const logK = LvArray::math::log( kValues[ic] ); + logTestComposition[ic] = LvArray::math::log( composition[ic] ) + (1.0/alpha) * logK; + } + } + + // Start iterations for this sample + bool converged = false; + + localIndex iterationCount = 0; + + // Start with SSI iterations + for(; (!converged) && (iterationCount < maxIterations); ++iterationCount ) + { + for( integer const ic : presentComponents ) + { + residual[ic] = -hyperplane[ic]; + } + real64 error = calculateResidual( numComps, + presentComponents, + pressure, + temperature, + componentProperties, + incipientEos, + flashData, + logTestComposition.toSliceConst(), + tpd, + testComposition, + logFugacity, + residual ); + // Update to next step + for( integer const ic : presentComponents ) + { + logTestComposition[ic] -= residual[ic]; + } + + // Check stationarity + converged = ( error < stabilityTolerance ); + if( error < stabilitySSITolerance ) + { + break; + } + } + + // Start with Newton iterations + for( iterationCount++; (!converged) && (iterationCount < maxIterations); ++iterationCount ) + { + for( integer const ic : presentComponents ) + { + residual[ic] = -hyperplane[ic]; + } + real64 error = calculateResidualAndJacobian( numComps, + presentComponents, + pressure, + temperature, + componentProperties, + incipientEos, + flashData, + logTestComposition.toSliceConst(), + tpd, + testComposition, + logFugacity, + derivatives.toSlice(), + residual, + jacobian.toSlice() ); + + // Ensure that the Jacobian is well formed for the absent components + for( integer const ic : absentComponents ) + { + jacobian( ic, ic ) = 1.0; + } + solveLinearSystem( jacobian.toSlice(), residual ); + // Update to next step + for( integer const ic : presentComponents ) + { + logTestComposition[ic] -= residual[ic]; + } + + // Check stationarity + converged = ( error < stabilityTolerance ); + } + + allConverged = allConverged && converged; + // Calculate the tangent-plane-distance (TPD) and distance to the trivial solution + real64 distance = 0.0; + for( integer const ic : presentComponents ) + { + real64 const dZ = testComposition[ic] - composition[ic]; + distance += (dZ*dZ); + } + + // The mixture is unstable if either the TPD is negative or if the TPD is zero but the incipient composition + // is different from the sample composition + // stabilityThreshold is negative (default -1e-8) + if( tpd < stabilityThreshold ) + { + unstableMixture = true; + } + else if((tpd < -stabilityThreshold) && (-stabilityThreshold < distance)) + { + unstableMixture = true; + } + + if( unstableMixture ) + { + storeIncipientComposition( incipientEos ); + return true; + } + else if( maxDistanceToTrivialSolution < distance ) + { + maxDistanceToTrivialSolution = distance; + storeIncipientComposition( incipientEos ); + } + } + } + + // The test is successful if either we have an unstable mixture or all test compositions converged to stationarity + return unstableMixture || allConverged; +} + +template< integer USD1, integer USD2 > +GEOS_HOST_DEVICE +void StabilityTest::computeDerivatives( integer const numComps, + real64 const pressure, + real64 const temperature, + arraySlice1d< real64 const, USD1 > const & composition, + ComponentProperties::KernelWrapper const & componentProperties, + EquationOfStateType const & sampleEquationOfState, + EquationOfStateType const & incipientEquationOfState, + FlashData const & flashData, + arraySlice1d< real64 const > const & incipientComposition, + arraySlice2d< real64, USD2 > const & incipientCompositionDerivs, + arraySlice2d< real64, USD2 > const & compositionDerivs ) +{ + integer constexpr maxNumRows = MultiFluidConstants::MAX_NUM_COMPONENTS + 1; + + integer const numDofs = 2 + numComps; + + StackArray< real64, 2, 3*maxNumComps > workSpace( 3, numComps ); + arraySlice1d< real64 > logFugacitySample = workSpace[0]; + arraySlice1d< real64 > logFugacityIncipient = workSpace[1]; + arraySlice1d< real64 > normalisedComposition = workSpace[2]; + StackArray< real64, 2, maxNumRows * maxNumRows, MatrixLayout::COL_MAJOR_PERM > A( numComps + 1, numComps + 1 ); + StackArray< real64, 2, maxNumRows * maxNumDofs, MatrixLayout::COL_MAJOR_PERM > X( numComps + 1, numDofs ); + + auto const & logFugacityIncipientDerivs = incipientCompositionDerivs; + auto const & logFugacitySampleDerivs = compositionDerivs; + + // Normalise the composition vector + copyComposition( numComps, composition, normalisedComposition ); + real64 const totalAmount = normalizeComposition( numComps, normalisedComposition ); + real64 const invTotalAmount = totalAmount < MultiFluidConstants::epsilon ? 0.0 : 1.0 / totalAmount; + FugacityCalculator::computeLogFugacityDerivatives( numComps, + pressure, + temperature, + normalisedComposition.toSliceConst(), + componentProperties, + sampleEquationOfState, + flashData, + logFugacitySample, + logFugacitySampleDerivs ); + + // Include the normalisation in the composition derivatives (chain rule): + // dx_j/dz_k = (delta_jk - x_j) / sum(z) + // dlnphi_i/dz_k = ( dlnphi_i/dx_k - sum_j( dlnphi_i/dx_j x_j ) ) / sum(z) + for( integer ic = 0; ic < numComps; ++ic ) + { + real64 weightedSum = 0.0; + for( integer jc = 0; jc < numComps; ++jc ) + { + weightedSum += logFugacitySampleDerivs( ic, Deriv::dC+jc ) * normalisedComposition[jc]; + } + for( integer kc = 0; kc < numComps; ++kc ) + { + logFugacitySampleDerivs( ic, Deriv::dC+kc ) = (logFugacitySampleDerivs( ic, Deriv::dC+kc ) - weightedSum) * invTotalAmount; + } + } + + FugacityCalculator::computeLogFugacityDerivatives( numComps, + pressure, + temperature, + incipientComposition, + componentProperties, + incipientEquationOfState, + flashData, + logFugacityIncipient, + logFugacityIncipientDerivs ); + + integer const refComp = getReferenceComponent( numComps, incipientComposition ); + real64 const logRatioRef = refComp < 0 ? 0.0 : logFugacitySample[refComp] - logFugacityIncipient[refComp]; + real64 const yOverZRef = refComp < 0 ? 0.0 : incipientComposition[refComp] / composition[refComp]; + + for( integer ic = 0; ic < numComps; ++ic ) + { + real64 const yi = incipientComposition[ic]; + for( integer jc = 0; jc < numComps; ++jc ) + { + real64 const df_dyj = logFugacityIncipientDerivs( ic, Deriv::dC+jc ); + A( ic, jc ) = yi*df_dyj; + } + A( ic, ic ) += 1.0; + A( ic, numComps ) = yi; + A( numComps, ic ) = 1.0; + } + A( numComps, numComps ) = 0.0; + + for( integer ic = 0; ic < numComps; ++ic ) + { + real64 const yi = incipientComposition[ic]; + for( integer const idof : {Deriv::dP, Deriv::dT} ) + { + X( ic, idof ) = yi*(-logFugacityIncipientDerivs( ic, idof ) + logFugacitySampleDerivs( ic, idof )); + } + for( integer jc = 0; jc < numComps; ++jc ) + { + integer const idof = Deriv::dC + jc; + X( ic, idof ) = yi*logFugacitySampleDerivs( ic, idof ); + } + if( MultiFluidConstants::minForSpeciesPresence < composition[ic] ) + { + X( ic, Deriv::dC+ic ) += yi/composition[ic]; + } + else + { + // Absent component: y_k = r_k z_k => dy_k/dz_k = r_k (one-sided derivative at z_k = 0) + real64 const logRatio = logFugacitySample[ic] - logFugacityIncipient[ic]; + X( ic, Deriv::dC+ic ) = yOverZRef * LvArray::math::exp( logRatio - logRatioRef ); + } + } + for( integer idof = 0; idof < numDofs; ++idof ) + { + X( numComps, idof ) = 0.0; + } + + // Solve linear system + solveLinearSystem( A.toSlice(), X.toSlice() ); + + for( integer idof = 0; idof < numDofs; ++idof ) + { + for( integer ic = 0; ic < numComps; ++ic ) + { + incipientCompositionDerivs( ic, idof ) = X( ic, idof ); + } + } +} + +template< int USD1, int USD2, int USD3 > +GEOS_HOST_DEVICE +real64 StabilityTest::calculateResidual( integer const numComps, + arraySlice1d< integer const > const & presentComponents, + real64 const pressure, + real64 const temperature, + ComponentProperties::KernelWrapper const & componentProperties, + EquationOfStateType const equationOfState, + FlashData const & flashData, + arraySlice1d< real64 const, USD1 > const & logTestComposition, + real64 & tangentPlaneDistance, + arraySlice1d< real64, USD2 > const & testComposition, + arraySlice1d< real64 > const & testLogFugacity, + arraySlice1d< real64, USD3 > const & residual ) +{ + for( integer const ic : presentComponents ) + { + testComposition[ic] = LvArray::math::exp( logTestComposition[ic] ); + } + real64 const sumY = normalizeComposition( numComps, testComposition ); + + FugacityCalculator::computeLogFugacity( numComps, + pressure, + temperature, + testComposition.toSliceConst(), + componentProperties, + equationOfState, + flashData, + testLogFugacity ); + + real64 error = 0.0; + tangentPlaneDistance = 1.0; + for( integer const ic : presentComponents ) + { + residual[ic] += logTestComposition[ic] + testLogFugacity[ic]; + tangentPlaneDistance += sumY * testComposition[ic] * (residual[ic] - 1.0); + error += (residual[ic]*residual[ic]); + } + error = LvArray::math::sqrt( error ); + return error; +} + +template< int USD1, int USD2, int USD3, int USD4 > +GEOS_HOST_DEVICE +real64 StabilityTest::calculateResidualAndJacobian( integer const numComps, + arraySlice1d< integer const > const & presentComponents, + real64 const pressure, + real64 const temperature, + ComponentProperties::KernelWrapper const & componentProperties, + EquationOfStateType const equationOfState, + FlashData const & flashData, + arraySlice1d< real64 const, USD1 > const & logTestComposition, + real64 & tangentPlaneDistance, + arraySlice1d< real64, USD2 > const & testComposition, + arraySlice1d< real64 > const & testLogFugacity, + arraySlice2d< real64 > const & logTestFugacityDerivs, + arraySlice1d< real64, USD3 > const & residual, + arraySlice2d< real64, USD4 > const & jacobian ) +{ + for( integer const ic : presentComponents ) + { + testComposition[ic] = LvArray::math::exp( logTestComposition[ic] ); + } + real64 const sumY = normalizeComposition( numComps, testComposition ); + + FugacityCalculator::computeLogFugacityDerivatives( numComps, + pressure, + temperature, + testComposition.toSliceConst(), + componentProperties, + equationOfState, + flashData, + testLogFugacity, + logTestFugacityDerivs ); + + LvArray::forValuesInSlice( jacobian, setZero ); + + tangentPlaneDistance = 1.0; + real64 error = 0.0; + for( integer const ic : presentComponents ) + { + for( integer const jc : presentComponents ) + { + real64 const yj = testComposition[jc]; + + real64 const dphi_ic_dyjc = logTestFugacityDerivs( ic, Deriv::dC + jc ); + + real64 sum_dphi_ic_dykc = dphi_ic_dyjc; + + for( integer const kc : presentComponents ) + { + real64 const yk = testComposition[kc]; + real64 const dphi_ic_dykc = logTestFugacityDerivs( ic, Deriv::dC + kc ); + sum_dphi_ic_dykc -= (dphi_ic_dykc * yk); + } + jacobian( ic, jc ) = yj * sum_dphi_ic_dykc; + } + jacobian( ic, ic ) += 1.0; + + residual[ic] += logTestComposition[ic] + testLogFugacity[ic]; + tangentPlaneDistance += sumY * testComposition[ic] * (residual[ic] - 1.0); + error += (residual[ic]*residual[ic]); + } + error = LvArray::math::sqrt( error ); + return error; +} + +} // namespace compositional + +} // namespace constitutive + +} // namespace geos + +#endif //GEOS_CONSTITUTIVE_FLUID_MULTIFLUID_COMPOSITIONAL_FUNCTIONS_STABILITYTEST_IMPL_HPP_ diff --git a/src/coreComponents/constitutive/fluid/multifluid/compositional/functions/Utilities.hpp b/src/coreComponents/constitutive/fluid/multifluid/compositional/functions/Utilities.hpp index d0430020bc7..89e483b4730 100644 --- a/src/coreComponents/constitutive/fluid/multifluid/compositional/functions/Utilities.hpp +++ b/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 ); diff --git a/src/coreComponents/constitutive/fluid/multifluid/compositional/models/NegativeTwoPhaseFlashModel.hpp b/src/coreComponents/constitutive/fluid/multifluid/compositional/models/NegativeTwoPhaseFlashModel.hpp index 9370e92a70c..955951c201d 100644 --- a/src/coreComponents/constitutive/fluid/multifluid/compositional/models/NegativeTwoPhaseFlashModel.hpp +++ b/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 ) diff --git a/src/coreComponents/constitutive/fluid/multifluid/compositional/parameters/FlashParameters.cpp b/src/coreComponents/constitutive/fluid/multifluid/compositional/parameters/FlashParameters.cpp index 7ecf7542d18..ab40bbbd80a 100644 --- a/src/coreComponents/constitutive/fluid/multifluid/compositional/parameters/FlashParameters.cpp +++ b/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() ); diff --git a/src/coreComponents/constitutive/fluid/multifluid/compositional/parameters/FlashParameters.hpp b/src/coreComponents/constitutive/fluid/multifluid/compositional/parameters/FlashParameters.hpp index fe430d83bbc..3e225e6508b 100644 --- a/src/coreComponents/constitutive/fluid/multifluid/compositional/parameters/FlashParameters.hpp +++ b/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; diff --git a/src/coreComponents/constitutive/unitTests/CMakeLists.txt b/src/coreComponents/constitutive/unitTests/CMakeLists.txt index 94276e254da..916034e85a9 100644 --- a/src/coreComponents/constitutive/unitTests/CMakeLists.txt +++ b/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 ) diff --git a/src/coreComponents/constitutive/unitTests/testNegativeTwoPhaseFlash9Comp.cpp b/src/coreComponents/constitutive/unitTests/testNegativeTwoPhaseFlash9Comp.cpp index 00a60a52b7a..48b350819bf 100644 --- a/src/coreComponents/constitutive/unitTests/testNegativeTwoPhaseFlash9Comp.cpp +++ b/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, diff --git a/src/coreComponents/constitutive/unitTests/testNegativeTwoPhaseFlashJacobian.cpp b/src/coreComponents/constitutive/unitTests/testNegativeTwoPhaseFlashJacobian.cpp new file mode 100644 index 00000000000..b4d52b40319 --- /dev/null +++ b/src/coreComponents/constitutive/unitTests/testNegativeTwoPhaseFlashJacobian.cpp @@ -0,0 +1,288 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +// Source includes +#include "codingUtilities/UnitTestUtilities.hpp" +#include "constitutive/fluid/multifluid/compositional/functions/NegativeTwoPhaseFlash.hpp" +#include "constitutive/fluid/multifluid/compositional/functions/KValueInitialization.hpp" +#include "TestFluid.hpp" +#include "TestFluidUtilities.hpp" + +using namespace geos::constitutive::compositional; + +namespace geos +{ +namespace testing +{ + +template< int NC > +struct FluidData +{}; + +template<> +struct FluidData< 2 > +{ + static std::unique_ptr< TestFluid< 2 > > createFluid() + { + auto fluid = TestFluid< 2 >::create( {Fluid::CO2, Fluid::C5H12} ); + return fluid; + } +}; + +template<> +struct FluidData< 4 > +{ + static std::unique_ptr< TestFluid< 4 > > createFluid() + { + return TestFluid< 4 >::create( {Fluid::N2, Fluid::CH4, Fluid::CO2, Fluid::H2O} ); + } +}; + +template<> +struct FluidData< 9 > +{ + static std::unique_ptr< TestFluid< 9 > > createFluid() + { + std::unique_ptr< TestFluid< 9 > > fluid = TestFluid< 9 >::create( {0, 0, 0, 0, 0, 0, 0, 0, 0} ); + // Manually populate + TestFluid< 9 >::populateArray( fluid->criticalPressure, Feed< 9 >{73.8659e5, 33.9439e5, 46.0421e5, 48.8387e5, 42.4552e5, 37.47e5, 33.5892e5, 30.1037e5, 20.549e5} ); + TestFluid< 9 >::populateArray( fluid->criticalTemperature, Feed< 9 >{304.7, 126.2, 190.6, 305.43, 369.8, 419.5, 465.9, 507.5, 678.8} ); + TestFluid< 9 >::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< 9 >::populateArray( fluid->acentricFactor, Feed< 9 >{0.225, 0.04, 0.013, 0.0986, 0.1524, 0.1956, 0.2413, 0.299, 0.5618} ); + TestFluid< 9 >::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< 9 >::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, + 0.0000e+00, 1.0000e-02, 0.0000e+00, + 0.0000e+00, 1.0000e-02, 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, 0.0000e+00, 0.0000e+00, + 1.0000e-02, 0.0000e+00, 2.8000e-02, 1.0000e-02, 1.0000e-02, 0.0000e+00, 0.0000e+00, + 1.0000e-02, 0.0000e+00, 4.5320e-02, 1.0000e-02, 1.0000e-02, 0.0000e+00, 0.0000e+00, 0.0000e+00 + } ); + return fluid; + } +}; + +template< int NC > +using FlashData = std::tuple< + real64 const, // pressure + real64 const, // temperature + Feed< NC > const, // total composition + real64 const // vapour fraction + >; + +template< int NC, EquationOfStateType LIQUID_EOS_TYPE, EquationOfStateType VAPOUR_EOS_TYPE > +class NegativeTwoPhaseFlashJacobianTestFixture : public ::testing::TestWithParam< FlashData< NC > > +{ + static constexpr real64 relTol = 1.0e-5; + static constexpr real64 absTol = 1.0e-7; + static constexpr int numComps = NC; + static constexpr int numDofs = NC + 2; + using Deriv = geos::constitutive::multifluid::DerivativeOffset; +public: + NegativeTwoPhaseFlashJacobianTestFixture() + : m_fluid( FluidData< NC >::createFluid() ) + {} + + ~NegativeTwoPhaseFlashJacobianTestFixture() = default; + + void testJacobian( FlashData< NC > const & data ) + { + // Number of values in each column + constexpr integer numValues = 1 + numComps; + + auto componentProperties = this->m_fluid->createKernelWrapper(); + + constitutive::compositional::FlashData flashData; + flashData.liquidEos = LIQUID_EOS_TYPE; + flashData.vapourEos = VAPOUR_EOS_TYPE; + + real64 const pressure = std::get< 0 >( data ); + real64 const temperature = std::get< 1 >( data ); + stackArray1d< real64, numComps > composition; + TestFluid< NC >::createArray( composition, std::get< 2 >( data )); + real64 const vapourFraction = std::get< 3 >( data ); + + StackArray< real64, 2, 5*numComps > workSpace( 5, numComps ); + StackArray< real64, 3, 2*numComps *numDofs > derivativeWorkSpace( 2, numComps, numDofs ); + StackArray< real64, 3, numValues *numValues > jacobianSpace( 1, numValues, numValues ); + StackArray< real64, 2, 3*numValues > residualSpace( 3, numValues ); + arraySlice1d< real64 > liquidComposition = workSpace[0]; + arraySlice1d< real64 > vapourComposition = workSpace[1]; + arraySlice1d< real64 > logLiquidFugacity = workSpace[2]; + arraySlice1d< real64 > logVapourFugacity = workSpace[3]; + arraySlice1d< real64 > kValues = workSpace[4]; + arraySlice2d< real64 > logLiquidFugacityDerivs = derivativeWorkSpace[0]; + arraySlice2d< real64 > logVapourFugacityDerivs = derivativeWorkSpace[1]; + arraySlice2d< real64 > jacobian = jacobianSpace[0]; + arraySlice1d< real64 > residual = residualSpace[0]; + arraySlice1d< real64 > displacedResidual = residualSpace[1]; + arraySlice1d< real64 > derivatives = residualSpace[2]; + + // Extract column + auto const extractColumn = [&]( integer const kc, real64 const scale = 1.0 ) + { + for( integer ic = 0; ic < numValues; ++ic ) + { + derivatives[ic] = scale * jacobian( ic, kc ); + } + }; + + KValueInitialization::computeWilsonGasLiquidKvalue( numComps, + pressure, + temperature, + componentProperties, + kValues ); + + NegativeTwoPhaseFlash::calculateResidualAndJacobian( numComps, + pressure, + temperature, + composition.toSliceConst(), + componentProperties, + flashData, + kValues.toSliceConst(), + vapourFraction, + liquidComposition, + vapourComposition, + logLiquidFugacity, + logVapourFugacity, + logLiquidFugacityDerivs, + logVapourFugacityDerivs, + residual, + jacobian ); + + // k-value derivatives + // Derivatives are wrt to log(K) + real64 constexpr deltaLogK = 1.0e-3; + for( integer kc = 0; kc < numComps; ++kc ) + { + extractColumn( kc ); + geos::testing::internal::testNumericalDerivative< numValues >( + 0.0, deltaLogK, derivatives, + [&]( real64 const dK, auto & values ) + { + real64 const kvOld = kValues[kc]; + kValues[kc] *= LvArray::math::exp( dK ); + NegativeTwoPhaseFlash::calculateResidual( numComps, + pressure, + temperature, + composition.toSliceConst(), + componentProperties, + flashData, + kValues.toSliceConst(), + vapourFraction, + liquidComposition, + vapourComposition, + logLiquidFugacity, + logVapourFugacity, + displacedResidual ); + for( integer ic = 0; ic < numValues; ++ic ) + { + values[ic] = displacedResidual[ic]; + } + kValues[kc] = kvOld; + } ); + } + + // Vapour fraction derivative + real64 constexpr deltaV = 1.0e-7; + extractColumn( numComps ); + geos::testing::internal::testNumericalDerivative< numValues >( + 0.0, deltaV, derivatives, + [&]( real64 const dV, auto & values ) + { + NegativeTwoPhaseFlash::calculateResidual( numComps, + pressure, + temperature, + composition.toSliceConst(), + componentProperties, + flashData, + kValues.toSliceConst(), + vapourFraction + dV, + liquidComposition, + vapourComposition, + logLiquidFugacity, + logVapourFugacity, + displacedResidual ); + for( integer ic = 0; ic < numValues; ++ic ) + { + values[ic] = displacedResidual[ic]; + } + } ); + } + +protected: + std::unique_ptr< TestFluid< NC > > m_fluid{}; +}; + +//using PengRobinson4 = NegativeTwoPhaseFlashJacobianTestFixture< 4, EquationOfStateType::PengRobinson, EquationOfStateType::PengRobinson >; +//using SoreideWhitson4 = NegativeTwoPhaseFlashJacobianTestFixture< 4, EquationOfStateType::SoreideWhitson, +// EquationOfStateType::PengRobinson >; +using PengRobinson9 = NegativeTwoPhaseFlashJacobianTestFixture< 9, EquationOfStateType::PengRobinson, EquationOfStateType::PengRobinson >; + +/* + TEST_P( PengRobinson4, testJacobian ) + { + testJacobian( GetParam() ); + } + + TEST_P( SoreideWhitson4, testJacobian ) + { + testJacobian( GetParam() ); + } + */ +TEST_P( PengRobinson9, testJacobian ) +{ + testJacobian( GetParam() ); +} + +//------------------------------------------------------------------------------- +// Data +//------------------------------------------------------------------------------- + +/* UNCRUSTIFY-OFF */ + +/* +INSTANTIATE_TEST_SUITE_P(NegativeTwoPhaseFlashJacobian, PengRobinson4, + ::testing::ValuesIn>({ + {1.00000e+06, 193.15, {0.0, 0.0, 0.0, 1.0}, 0.5} + }) +); + +INSTANTIATE_TEST_SUITE_P(NegativeTwoPhaseFlashJacobian, SoreideWhitson4, + ::testing::ValuesIn>({ + {1.00000e+06, 193.15, {0.0, 0.0, 0.0, 1.0}, 0.0}, + {1.00000e+06, 193.15, {0.0, 0.0, 0.0, 1.0}, 0.5}, + {1.00000e+06, 193.15, {0.0, 0.0, 0.0, 1.0}, 1.0}, + {1.00000e+06, 193.15, {0.0, 0.0, 0.0, 1.0}, 0.998}, + {1.00000e+06, 193.15, {0.0, 0.0, 0.0, 1.0}, 1.002} + }) +); +*/ + +INSTANTIATE_TEST_SUITE_P(NegativeTwoPhaseFlashJacobian, PengRobinson9, + ::testing::ValuesIn>({ + //{1.00000e+07, 278.15, {0.000363, 0.000007, 0.003471, 0.006007, 0.018423, 0.034034, 0.042565, 0.056120, 0.839010}, -0.18866442} + //{1.00000e+07, 278.15, {0.000363, 0.000007, 0.003471, 0.006007, 0.018423, 0.034034, 0.042565, 0.056120, 0.839010}, 0.18866442} + {1.00000e+07, 278.15, {0.000363, 0.000007, 0.003471, 0.006007, 0.018423, 0.034034, 0.042565, 0.056120, 0.839010}, -0.01} + }) +); + +/* UNCRUSTIFY-ON */ + +} // testing +} // geos diff --git a/src/coreComponents/constitutive/unitTests/testNegativeTwoPhaseFlashModel.cpp b/src/coreComponents/constitutive/unitTests/testNegativeTwoPhaseFlashModel.cpp index cd051c5e38d..04bc00fcbbe 100644 --- a/src/coreComponents/constitutive/unitTests/testNegativeTwoPhaseFlashModel.cpp +++ b/src/coreComponents/constitutive/unitTests/testNegativeTwoPhaseFlashModel.cpp @@ -131,6 +131,7 @@ class NegativeTwoPhaseFlashModelTestFixture : public ::testing::TestWithParam< flashParameters->m_continuousParameters[FlashParameters::STABILITY_TOLERANCE] = 1.0e-12; flashParameters->m_continuousParameters[FlashParameters::STABILITY_THRESHOLD] = -1.0e-3; flashParameters->m_discreteParameters[FlashParameters::FLASH_MAX_ITERATIONS] = 20; + flashParameters->m_continuousParameters[FlashParameters::SSI_TOLERANCE] = 1.0e-12; m_phaseTypes.emplace_back( static_cast< integer >(PhaseType::LIQUID)); m_phaseTypes.emplace_back( static_cast< integer >(PhaseType::VAPOUR)); diff --git a/src/coreComponents/constitutive/unitTests/testStabilityTest2Comp.cpp b/src/coreComponents/constitutive/unitTests/testStabilityTest2Comp.cpp index 93e9b965716..d25f59f4653 100644 --- a/src/coreComponents/constitutive/unitTests/testStabilityTest2Comp.cpp +++ b/src/coreComponents/constitutive/unitTests/testStabilityTest2Comp.cpp @@ -81,6 +81,7 @@ class StabilityTestTest2CompFixture : public ::testing::TestWithParam< Stabilit auto const parameters = FlashParameters::create( std::make_unique< ModelParameters >() ); auto const * flashParameters = parameters->get< FlashParameters >(); + flashParameters->m_discreteParameters[FlashParameters::STABILITY_MAX_ITERATIONS] = 20; bool const stabilityStatus = StabilityTest::compute( numComps, pressure, @@ -143,8 +144,8 @@ INSTANTIATE_TEST_SUITE_P(StabilityTest, PengRobinson, {5.0e+06, 353.15, 0.4000, false, 0.40000000 }, {2.0e+07, 297.15, 0.4000, false, 0.40000000 }, {2.0e+07, 353.15, 0.4000, false, 0.40000001 }, - {5.0e+06, 297.15, 0.7860, true, 0.31984391 }, - {5.0e+06, 297.15, 0.7870, true, 0.32331543 }, + {5.0e+06, 297.15, 0.7860, false, 0.31984391 }, + {5.0e+06, 297.15, 0.7870, false, 0.32331543 }, {5.0e+06, 297.15, 0.3200, true, 0.78604645 }, {5.0e+06, 297.15, 0.3220, true, 0.78662927 } }) diff --git a/src/coreComponents/constitutive/unitTests/testStabilityTest9Comp.cpp b/src/coreComponents/constitutive/unitTests/testStabilityTest9Comp.cpp index 07e3ee5bf41..33194bb05f9 100644 --- a/src/coreComponents/constitutive/unitTests/testStabilityTest9Comp.cpp +++ b/src/coreComponents/constitutive/unitTests/testStabilityTest9Comp.cpp @@ -82,6 +82,7 @@ class StabilityTestTest9CompFixture : public ::testing::TestWithParam< Stabilit auto const parameters = FlashParameters::create( std::make_unique< ModelParameters >() ); auto const * flashParameters = parameters->get< FlashParameters >(); + setMaximumNumberOfIterations( flashParameters, 20 ); bool const stabilityStatus = StabilityTest::compute( numComps, pressure, @@ -99,7 +100,7 @@ class StabilityTestTest9CompFixture : public ::testing::TestWithParam< Stabilit // Expect this to succeed ASSERT_EQ( stabilityStatus, true ); - // Check the stability labe; + // Check the stability label ASSERT_EQ( unstableMixture, expectedUnstableMixture ); // Check the incipient methane fraction @@ -138,6 +139,7 @@ class StabilityTestTest9CompFixture : public ::testing::TestWithParam< Stabilit auto const parameters = FlashParameters::create( std::make_unique< ModelParameters >() ); auto const * flashParameters = parameters->get< FlashParameters >(); + setMaximumNumberOfIterations( flashParameters, 20 ); StabilityTest::compute( numComps, pressure, @@ -167,33 +169,45 @@ class StabilityTestTest9CompFixture : public ::testing::TestWithParam< Stabilit liquidComposition, vapourComposition ); - // Expect this to succeed - ASSERT_EQ( status, true ); - - constexpr real64 epsilon = MultiFluidConstants::epsilon; - real64 const V = vapourFraction; - real64 const L = 1.0 - V; + // Expect the flash to succeed for an unstable mixture + ASSERT_EQ( status || !unstableMixture, true ); - // Check stability vs flash - if( unstableMixture ) + if( status ) { - // Unstable mixture: both L and V should be positive - ASSERT_GT( V, epsilon ); - ASSERT_GT( L, epsilon ); - } - else - { - // Stable mixture: one of L or V should be zero - ASSERT_TRUE( V < epsilon || L < epsilon ); + constexpr real64 epsilon = MultiFluidConstants::epsilon; + real64 const V = vapourFraction; + real64 const L = 1.0 - V; + + // Check stability vs flash + if( unstableMixture ) + { + // Unstable mixture: both L and V should be positive + ASSERT_GT( V, epsilon ); + ASSERT_GT( L, epsilon ); + } + else + { + // Stable mixture: one of L or V should be zero + ASSERT_TRUE( V < epsilon || L < epsilon ); + } } } protected: std::unique_ptr< TestFluid< numComps > > m_fluid{}; private: + static void setMaximumNumberOfIterations( FlashParameters const * parameters, integer const maxIterations ); static std::unique_ptr< TestFluid< numComps > > createFluid(); }; +template< EquationOfStateType EOS_TYPE > +void StabilityTestTest9CompFixture< EOS_TYPE >::setMaximumNumberOfIterations( FlashParameters const * parameters, + integer const maxIterations ) +{ + auto * flashParameters = const_cast< FlashParameters * >(parameters); + flashParameters->m_discreteParameters[FlashParameters::STABILITY_MAX_ITERATIONS] = maxIterations; +} + template< EquationOfStateType EOS_TYPE > std::unique_ptr< TestFluid< numComps > > StabilityTestTest9CompFixture< EOS_TYPE >::createFluid() { @@ -220,6 +234,7 @@ std::unique_ptr< TestFluid< numComps > > StabilityTestTest9CompFixture< EOS_TYPE using PengRobinson = StabilityTestTest9CompFixture< EquationOfStateType::PengRobinson >; using SoaveRedlichKwong = StabilityTestTest9CompFixture< EquationOfStateType::SoaveRedlichKwong >; + TEST_P( PengRobinson, testStabilityTest ) { testStability( GetParam() ); @@ -246,22 +261,118 @@ TEST_P( SoaveRedlichKwong, testStabilityWithFlash ) INSTANTIATE_TEST_SUITE_P( StabilityTest, PengRobinson, ::testing::ValuesIn({ - {1.0e+05, 288.15, { 0.0090, 0.0030, 0.5347, 0.1146, 0.0879, 0.0456, 0.0209, 0.0151, 0.1692 }, true, 0.99992981}, - {1.0e+05, 297.15, { 0.0090, 0.0030, 0.5347, 0.1146, 0.0879, 0.0456, 0.0209, 0.0151, 0.1692 }, true, 0.99989343}, - {1.0e+05, 353.15, { 0.0090, 0.0030, 0.5347, 0.1146, 0.0879, 0.0456, 0.0209, 0.0151, 0.1692 }, true, 0.99906841}, - {1.0e+05, 393.15, { 0.0090, 0.0030, 0.5347, 0.1146, 0.0879, 0.0456, 0.0209, 0.0151, 0.1692 }, true, 0.99685491}, - {1.0e+05, 573.15, { 0.0090, 0.0030, 0.5347, 0.1146, 0.0879, 0.0456, 0.0209, 0.0151, 0.1692 }, false, 0.16920000} + {1.0e+05, 288.15, {9.0000000e-03, 3.0000000e-03, 5.3470000e-01, 1.1460000e-01, 8.7900000e-02, 4.5600000e-02, 2.0900000e-02, 1.5100000e-02, 1.6920000e-01}, true, 9.99930e-01}, + {1.0e+05, 298.15, {9.0000000e-03, 3.0000000e-03, 5.3470000e-01, 1.1460000e-01, 8.7900000e-02, 4.5600000e-02, 2.0900000e-02, 1.5100000e-02, 1.6920000e-01}, true, 9.99889e-01}, + {1.0e+05, 353.15, {9.0000000e-03, 3.0000000e-03, 5.3470000e-01, 1.1460000e-01, 8.7900000e-02, 4.5600000e-02, 2.0900000e-02, 1.5100000e-02, 1.6920000e-01}, true, 9.99068e-01}, + {1.0e+05, 573.15, {9.0000000e-03, 3.0000000e-03, 5.3470000e-01, 1.1460000e-01, 8.7900000e-02, 4.5600000e-02, 2.0900000e-02, 1.5100000e-02, 1.6920000e-01}, false, 1.69200e-01}, + {5.0e+06, 288.15, {9.0000000e-03, 3.0000000e-03, 5.3470000e-01, 1.1460000e-01, 8.7900000e-02, 4.5600000e-02, 2.0900000e-02, 1.5100000e-02, 1.6920000e-01}, true, 3.31496e-06}, + {5.0e+06, 298.15, {9.0000000e-03, 3.0000000e-03, 5.3470000e-01, 1.1460000e-01, 8.7900000e-02, 4.5600000e-02, 2.0900000e-02, 1.5100000e-02, 1.6920000e-01}, true, 6.19514e-06}, + {5.0e+06, 353.15, {9.0000000e-03, 3.0000000e-03, 5.3470000e-01, 1.1460000e-01, 8.7900000e-02, 4.5600000e-02, 2.0900000e-02, 1.5100000e-02, 1.6920000e-01}, true, 1.68840e-04}, + {5.0e+06, 573.15, {9.0000000e-03, 3.0000000e-03, 5.3470000e-01, 1.1460000e-01, 8.7900000e-02, 4.5600000e-02, 2.0900000e-02, 1.5100000e-02, 1.6920000e-01}, true, 7.60017e-01}, + {1.0e+07, 288.15, {9.0000000e-03, 3.0000000e-03, 5.3470000e-01, 1.1460000e-01, 8.7900000e-02, 4.5600000e-02, 2.0900000e-02, 1.5100000e-02, 1.6920000e-01}, true, 4.69115e-05}, + {1.0e+07, 298.15, {9.0000000e-03, 3.0000000e-03, 5.3470000e-01, 1.1460000e-01, 8.7900000e-02, 4.5600000e-02, 2.0900000e-02, 1.5100000e-02, 1.6920000e-01}, true, 6.37278e-05}, + {1.0e+07, 353.15, {9.0000000e-03, 3.0000000e-03, 5.3470000e-01, 1.1460000e-01, 8.7900000e-02, 4.5600000e-02, 2.0900000e-02, 1.5100000e-02, 1.6920000e-01}, true, 4.72565e-04}, + {1.0e+07, 573.15, {9.0000000e-03, 3.0000000e-03, 5.3470000e-01, 1.1460000e-01, 8.7900000e-02, 4.5600000e-02, 2.0900000e-02, 1.5100000e-02, 1.6920000e-01}, true, 5.69337e-01}, + {1.0e+08, 288.15, {9.0000000e-03, 3.0000000e-03, 5.3470000e-01, 1.1460000e-01, 8.7900000e-02, 4.5600000e-02, 2.0900000e-02, 1.5100000e-02, 1.6920000e-01}, false, 1.69200e-01}, + {1.0e+08, 298.15, {9.0000000e-03, 3.0000000e-03, 5.3470000e-01, 1.1460000e-01, 8.7900000e-02, 4.5600000e-02, 2.0900000e-02, 1.5100000e-02, 1.6920000e-01}, false, 1.69200e-01}, + {1.0e+08, 353.15, {9.0000000e-03, 3.0000000e-03, 5.3470000e-01, 1.1460000e-01, 8.7900000e-02, 4.5600000e-02, 2.0900000e-02, 1.5100000e-02, 1.6920000e-01}, false, 1.69200e-01}, + {1.0e+08, 573.15, {9.0000000e-03, 3.0000000e-03, 5.3470000e-01, 1.1460000e-01, 8.7900000e-02, 4.5600000e-02, 2.0900000e-02, 1.5100000e-02, 1.6920000e-01}, false, 1.69200e-01}, + {1.0e+05, 288.15, {9.0000000e-03, 0.0000000e+00, 6.5230000e-01, 1.3300000e-05, 8.7900000e-02, 4.5600000e-02, 2.0900000e-02, 1.5100000e-02, 1.6918670e-01}, true, 9.99931e-01}, + {1.0e+05, 298.15, {9.0000000e-03, 0.0000000e+00, 6.5230000e-01, 1.3300000e-05, 8.7900000e-02, 4.5600000e-02, 2.0900000e-02, 1.5100000e-02, 1.6918670e-01}, true, 9.99891e-01}, + {1.0e+05, 353.15, {9.0000000e-03, 0.0000000e+00, 6.5230000e-01, 1.3300000e-05, 8.7900000e-02, 4.5600000e-02, 2.0900000e-02, 1.5100000e-02, 1.6918670e-01}, true, 9.99113e-01}, + {1.0e+05, 573.15, {9.0000000e-03, 0.0000000e+00, 6.5230000e-01, 1.3300000e-05, 8.7900000e-02, 4.5600000e-02, 2.0900000e-02, 1.5100000e-02, 1.6918670e-01}, false, 1.69187e-01}, + {5.0e+06, 288.15, {9.0000000e-03, 0.0000000e+00, 6.5230000e-01, 1.3300000e-05, 8.7900000e-02, 4.5600000e-02, 2.0900000e-02, 1.5100000e-02, 1.6918670e-01}, true, 3.10425e-06}, + {5.0e+06, 298.15, {9.0000000e-03, 0.0000000e+00, 6.5230000e-01, 1.3300000e-05, 8.7900000e-02, 4.5600000e-02, 2.0900000e-02, 1.5100000e-02, 1.6918670e-01}, true, 5.93265e-06}, + {5.0e+06, 353.15, {9.0000000e-03, 0.0000000e+00, 6.5230000e-01, 1.3300000e-05, 8.7900000e-02, 4.5600000e-02, 2.0900000e-02, 1.5100000e-02, 1.6918670e-01}, true, 2.11049e-04}, + {5.0e+06, 573.15, {9.0000000e-03, 0.0000000e+00, 6.5230000e-01, 1.3300000e-05, 8.7900000e-02, 4.5600000e-02, 2.0900000e-02, 1.5100000e-02, 1.6918670e-01}, true, 7.78486e-01}, + {1.0e+07, 288.15, {9.0000000e-03, 0.0000000e+00, 6.5230000e-01, 1.3300000e-05, 8.7900000e-02, 4.5600000e-02, 2.0900000e-02, 1.5100000e-02, 1.6918670e-01}, true, 3.30848e-05}, + {1.0e+07, 298.15, {9.0000000e-03, 0.0000000e+00, 6.5230000e-01, 1.3300000e-05, 8.7900000e-02, 4.5600000e-02, 2.0900000e-02, 1.5100000e-02, 1.6918670e-01}, true, 4.73576e-05}, + {1.0e+07, 353.15, {9.0000000e-03, 0.0000000e+00, 6.5230000e-01, 1.3300000e-05, 8.7900000e-02, 4.5600000e-02, 2.0900000e-02, 1.5100000e-02, 1.6918670e-01}, true, 4.31414e-04}, + {1.0e+07, 573.15, {9.0000000e-03, 0.0000000e+00, 6.5230000e-01, 1.3300000e-05, 8.7900000e-02, 4.5600000e-02, 2.0900000e-02, 1.5100000e-02, 1.6918670e-01}, true, 6.09508e-01}, + {1.0e+08, 288.15, {9.0000000e-03, 0.0000000e+00, 6.5230000e-01, 1.3300000e-05, 8.7900000e-02, 4.5600000e-02, 2.0900000e-02, 1.5100000e-02, 1.6918670e-01}, false, 1.69187e-01}, + {1.0e+08, 298.15, {9.0000000e-03, 0.0000000e+00, 6.5230000e-01, 1.3300000e-05, 8.7900000e-02, 4.5600000e-02, 2.0900000e-02, 1.5100000e-02, 1.6918670e-01}, false, 1.69187e-01}, + {1.0e+08, 353.15, {9.0000000e-03, 0.0000000e+00, 6.5230000e-01, 1.3300000e-05, 8.7900000e-02, 4.5600000e-02, 2.0900000e-02, 1.5100000e-02, 1.6918670e-01}, false, 1.69187e-01}, + {1.0e+08, 573.15, {9.0000000e-03, 0.0000000e+00, 6.5230000e-01, 1.3300000e-05, 8.7900000e-02, 4.5600000e-02, 2.0900000e-02, 1.5100000e-02, 1.6918670e-01}, false, 1.69187e-01}, + {1.0e+05, 288.15, {1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111200e-01}, true, 9.99380e-01}, + {1.0e+05, 298.15, {1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111200e-01}, true, 9.99048e-01}, + {1.0e+05, 353.15, {1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111200e-01}, true, 9.93558e-01}, + {1.0e+05, 573.15, {1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111200e-01}, false, 1.11112e-01}, + {5.0e+06, 288.15, {1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111200e-01}, true, 1.20037e-06}, + {5.0e+06, 298.15, {1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111200e-01}, true, 2.41217e-06}, + {5.0e+06, 353.15, {1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111200e-01}, true, 6.49045e-05}, + {5.0e+06, 573.15, {1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111200e-01}, false, 1.11112e-01}, + {1.0e+07, 288.15, {1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111200e-01}, false, 1.43584e-05}, + {1.0e+07, 298.15, {1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111200e-01}, true, 2.41543e-05}, + {1.0e+07, 353.15, {1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111200e-01}, true, 3.04002e-04}, + {1.0e+08, 288.15, {1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111200e-01}, false, 1.11112e-01}, + {1.0e+08, 298.15, {1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111200e-01}, false, 1.11112e-01}, + {1.0e+08, 353.15, {1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111200e-01}, false, 1.11112e-01}, + {1.0e+08, 573.15, {1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111200e-01}, false, 1.11112e-01}, + {1.0e+05, 573.15, {0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 1.0000000e+00}, false, 1.00000e+00}, + {1.0e+08, 288.15, {0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 1.0000000e+00}, false, 1.00000e+00}, + {1.0e+05, 573.15, {0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 1.3300000e-05, 1.3300000e-05, 9.9997340e-01}, false, 9.99973e-01}, + {1.0e+08, 288.15, {0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 1.3300000e-05, 1.3300000e-05, 9.9997340e-01}, false, 9.99973e-01}, + {1.0e+05, 573.15, {0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 1.3300000e-05, 9.9997340e-01, 1.3300000e-05}, false, 1.33000e-05}, + {1.0e+08, 288.15, {0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 1.3300000e-05, 9.9997340e-01, 1.3300000e-05}, false, 1.33000e-05} }) ); INSTANTIATE_TEST_SUITE_P( StabilityTest, SoaveRedlichKwong, ::testing::ValuesIn({ - {1.0e+05, 288.15, { 0.0090, 0.0030, 0.5347, 0.1146, 0.0879, 0.0456, 0.0209, 0.0151, 0.1692 }, true, 0.99994255}, - {1.0e+05, 297.15, { 0.0090, 0.0030, 0.5347, 0.1146, 0.0879, 0.0456, 0.0209, 0.0151, 0.1692 }, true, 0.99991141}, - {1.0e+05, 353.15, { 0.0090, 0.0030, 0.5347, 0.1146, 0.0879, 0.0456, 0.0209, 0.0151, 0.1692 }, true, 0.99916193}, - {1.0e+05, 393.15, { 0.0090, 0.0030, 0.5347, 0.1146, 0.0879, 0.0456, 0.0209, 0.0151, 0.1692 }, true, 0.99704625}, - {1.0e+05, 573.15, { 0.0090, 0.0030, 0.5347, 0.1146, 0.0879, 0.0456, 0.0209, 0.0151, 0.1692 }, false, 0.16920000} + {1.0e+05, 288.15, {9.0000000e-03, 3.0000000e-03, 5.3470000e-01, 1.1460000e-01, 8.7900000e-02, 4.5600000e-02, 2.0900000e-02, 1.5100000e-02, 1.6920000e-01}, true, 9.99943e-01}, + {1.0e+05, 298.15, {9.0000000e-03, 3.0000000e-03, 5.3470000e-01, 1.1460000e-01, 8.7900000e-02, 4.5600000e-02, 2.0900000e-02, 1.5100000e-02, 1.6920000e-01}, true, 9.99907e-01}, + {1.0e+05, 353.15, {9.0000000e-03, 3.0000000e-03, 5.3470000e-01, 1.1460000e-01, 8.7900000e-02, 4.5600000e-02, 2.0900000e-02, 1.5100000e-02, 1.6920000e-01}, true, 9.99162e-01}, + {1.0e+05, 573.15, {9.0000000e-03, 3.0000000e-03, 5.3470000e-01, 1.1460000e-01, 8.7900000e-02, 4.5600000e-02, 2.0900000e-02, 1.5100000e-02, 1.6920000e-01}, false, 1.69200e-01}, + {5.0e+06, 288.15, {9.0000000e-03, 3.0000000e-03, 5.3470000e-01, 1.1460000e-01, 8.7900000e-02, 4.5600000e-02, 2.0900000e-02, 1.5100000e-02, 1.6920000e-01}, true, 2.35604e-06}, + {5.0e+06, 298.15, {9.0000000e-03, 3.0000000e-03, 5.3470000e-01, 1.1460000e-01, 8.7900000e-02, 4.5600000e-02, 2.0900000e-02, 1.5100000e-02, 1.6920000e-01}, true, 4.52623e-06}, + {5.0e+06, 353.15, {9.0000000e-03, 3.0000000e-03, 5.3470000e-01, 1.1460000e-01, 8.7900000e-02, 4.5600000e-02, 2.0900000e-02, 1.5100000e-02, 1.6920000e-01}, true, 1.34208e-04}, + {5.0e+06, 573.15, {9.0000000e-03, 3.0000000e-03, 5.3470000e-01, 1.1460000e-01, 8.7900000e-02, 4.5600000e-02, 2.0900000e-02, 1.5100000e-02, 1.6920000e-01}, true, 7.61380e-01}, + {1.0e+07, 288.15, {9.0000000e-03, 3.0000000e-03, 5.3470000e-01, 1.1460000e-01, 8.7900000e-02, 4.5600000e-02, 2.0900000e-02, 1.5100000e-02, 1.6920000e-01}, true, 3.24422e-05}, + {1.0e+07, 298.15, {9.0000000e-03, 3.0000000e-03, 5.3470000e-01, 1.1460000e-01, 8.7900000e-02, 4.5600000e-02, 2.0900000e-02, 1.5100000e-02, 1.6920000e-01}, true, 4.49843e-05}, + {1.0e+07, 353.15, {9.0000000e-03, 3.0000000e-03, 5.3470000e-01, 1.1460000e-01, 8.7900000e-02, 4.5600000e-02, 2.0900000e-02, 1.5100000e-02, 1.6920000e-01}, true, 3.62686e-04}, + {1.0e+07, 573.15, {9.0000000e-03, 3.0000000e-03, 5.3470000e-01, 1.1460000e-01, 8.7900000e-02, 4.5600000e-02, 2.0900000e-02, 1.5100000e-02, 1.6920000e-01}, true, 5.80899e-01}, + {1.0e+08, 288.15, {9.0000000e-03, 3.0000000e-03, 5.3470000e-01, 1.1460000e-01, 8.7900000e-02, 4.5600000e-02, 2.0900000e-02, 1.5100000e-02, 1.6920000e-01}, false, 1.69200e-01}, + {1.0e+08, 298.15, {9.0000000e-03, 3.0000000e-03, 5.3470000e-01, 1.1460000e-01, 8.7900000e-02, 4.5600000e-02, 2.0900000e-02, 1.5100000e-02, 1.6920000e-01}, false, 1.69200e-01}, + {1.0e+08, 353.15, {9.0000000e-03, 3.0000000e-03, 5.3470000e-01, 1.1460000e-01, 8.7900000e-02, 4.5600000e-02, 2.0900000e-02, 1.5100000e-02, 1.6920000e-01}, false, 1.69200e-01}, + {1.0e+08, 573.15, {9.0000000e-03, 3.0000000e-03, 5.3470000e-01, 1.1460000e-01, 8.7900000e-02, 4.5600000e-02, 2.0900000e-02, 1.5100000e-02, 1.6920000e-01}, false, 1.69200e-01}, + {1.0e+05, 288.15, {9.0000000e-03, 0.0000000e+00, 6.5230000e-01, 1.3300000e-05, 8.7900000e-02, 4.5600000e-02, 2.0900000e-02, 1.5100000e-02, 1.6918670e-01}, true, 9.99944e-01}, + {1.0e+05, 298.15, {9.0000000e-03, 0.0000000e+00, 6.5230000e-01, 1.3300000e-05, 8.7900000e-02, 4.5600000e-02, 2.0900000e-02, 1.5100000e-02, 1.6918670e-01}, true, 9.99909e-01}, + {1.0e+05, 353.15, {9.0000000e-03, 0.0000000e+00, 6.5230000e-01, 1.3300000e-05, 8.7900000e-02, 4.5600000e-02, 2.0900000e-02, 1.5100000e-02, 1.6918670e-01}, true, 9.99202e-01}, + {1.0e+05, 573.15, {9.0000000e-03, 0.0000000e+00, 6.5230000e-01, 1.3300000e-05, 8.7900000e-02, 4.5600000e-02, 2.0900000e-02, 1.5100000e-02, 1.6918670e-01}, false, 1.69187e-01}, + {5.0e+06, 288.15, {9.0000000e-03, 0.0000000e+00, 6.5230000e-01, 1.3300000e-05, 8.7900000e-02, 4.5600000e-02, 2.0900000e-02, 1.5100000e-02, 1.6918670e-01}, true, 2.21223e-06}, + {5.0e+06, 298.15, {9.0000000e-03, 0.0000000e+00, 6.5230000e-01, 1.3300000e-05, 8.7900000e-02, 4.5600000e-02, 2.0900000e-02, 1.5100000e-02, 1.6918670e-01}, true, 4.33855e-06}, + {5.0e+06, 353.15, {9.0000000e-03, 0.0000000e+00, 6.5230000e-01, 1.3300000e-05, 8.7900000e-02, 4.5600000e-02, 2.0900000e-02, 1.5100000e-02, 1.6918670e-01}, true, 1.63405e-04}, + {5.0e+06, 573.15, {9.0000000e-03, 0.0000000e+00, 6.5230000e-01, 1.3300000e-05, 8.7900000e-02, 4.5600000e-02, 2.0900000e-02, 1.5100000e-02, 1.6918670e-01}, true, 7.79335e-01}, + {1.0e+07, 288.15, {9.0000000e-03, 0.0000000e+00, 6.5230000e-01, 1.3300000e-05, 8.7900000e-02, 4.5600000e-02, 2.0900000e-02, 1.5100000e-02, 1.6918670e-01}, true, 2.28707e-05}, + {1.0e+07, 298.15, {9.0000000e-03, 0.0000000e+00, 6.5230000e-01, 1.3300000e-05, 8.7900000e-02, 4.5600000e-02, 2.0900000e-02, 1.5100000e-02, 1.6918670e-01}, true, 3.33621e-05}, + {1.0e+07, 353.15, {9.0000000e-03, 0.0000000e+00, 6.5230000e-01, 1.3300000e-05, 8.7900000e-02, 4.5600000e-02, 2.0900000e-02, 1.5100000e-02, 1.6918670e-01}, true, 3.27432e-04}, + {1.0e+07, 573.15, {9.0000000e-03, 0.0000000e+00, 6.5230000e-01, 1.3300000e-05, 8.7900000e-02, 4.5600000e-02, 2.0900000e-02, 1.5100000e-02, 1.6918670e-01}, true, 6.18577e-01}, + {1.0e+08, 288.15, {9.0000000e-03, 0.0000000e+00, 6.5230000e-01, 1.3300000e-05, 8.7900000e-02, 4.5600000e-02, 2.0900000e-02, 1.5100000e-02, 1.6918670e-01}, false, 1.69187e-01}, + {1.0e+08, 298.15, {9.0000000e-03, 0.0000000e+00, 6.5230000e-01, 1.3300000e-05, 8.7900000e-02, 4.5600000e-02, 2.0900000e-02, 1.5100000e-02, 1.6918670e-01}, false, 1.69187e-01}, + {1.0e+08, 353.15, {9.0000000e-03, 0.0000000e+00, 6.5230000e-01, 1.3300000e-05, 8.7900000e-02, 4.5600000e-02, 2.0900000e-02, 1.5100000e-02, 1.6918670e-01}, false, 1.69187e-01}, + {1.0e+08, 573.15, {9.0000000e-03, 0.0000000e+00, 6.5230000e-01, 1.3300000e-05, 8.7900000e-02, 4.5600000e-02, 2.0900000e-02, 1.5100000e-02, 1.6918670e-01}, false, 1.69187e-01}, + {1.0e+05, 288.15, {1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111200e-01}, true, 9.99489e-01}, + {1.0e+05, 298.15, {1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111200e-01}, true, 9.99201e-01}, + {1.0e+05, 353.15, {1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111200e-01}, true, 9.94186e-01}, + {1.0e+05, 573.15, {1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111200e-01}, false, 1.11112e-01}, + {5.0e+06, 288.15, {1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111200e-01}, true, 8.06893e-07}, + {5.0e+06, 298.15, {1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111200e-01}, true, 1.68319e-06}, + {5.0e+06, 353.15, {1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111200e-01}, true, 5.24670e-05}, + {5.0e+06, 573.15, {1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111200e-01}, false, 1.11112e-01}, + {1.0e+07, 288.15, {1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111200e-01}, true, 8.46960e-06}, + {1.0e+07, 298.15, {1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111200e-01}, true, 1.49116e-05}, + {1.0e+07, 353.15, {1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111200e-01}, true, 2.25256e-04}, + {1.0e+08, 288.15, {1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111200e-01}, false, 1.11112e-01}, + {1.0e+08, 298.15, {1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111200e-01}, false, 1.11112e-01}, + {1.0e+08, 353.15, {1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111200e-01}, false, 1.11112e-01}, + {1.0e+08, 573.15, {1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111100e-01, 1.1111200e-01}, false, 1.11112e-01}, + {1.0e+05, 573.15, {0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 1.0000000e+00}, false, 1.00000e+00}, + {1.0e+08, 288.15, {0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 1.0000000e+00}, false, 1.00000e+00}, + {1.0e+05, 573.15, {0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 1.3300000e-05, 1.3300000e-05, 9.9997340e-01}, false, 9.99973e-01}, + {1.0e+08, 288.15, {0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 1.3300000e-05, 1.3300000e-05, 9.9997340e-01}, false, 9.99973e-01}, + {1.0e+05, 573.15, {0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 1.3300000e-05, 9.9997340e-01, 1.3300000e-05}, false, 1.33000e-05}, + {1.0e+08, 288.15, {0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 1.3300000e-05, 9.9997340e-01, 1.3300000e-05}, false, 1.33000e-05} }) ); diff --git a/src/coreComponents/constitutive/unitTests/testStabilityTestDerivatives.cpp b/src/coreComponents/constitutive/unitTests/testStabilityTestDerivatives.cpp new file mode 100644 index 00000000000..cc53cb3e004 --- /dev/null +++ b/src/coreComponents/constitutive/unitTests/testStabilityTestDerivatives.cpp @@ -0,0 +1,372 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +// Source includes +#include "codingUtilities/UnitTestUtilities.hpp" +#include "constitutive/fluid/multifluid/compositional/functions/StabilityTest.hpp" +#include "constitutive/fluid/multifluid/MultiFluidConstants.hpp" +#include "TestFluid.hpp" +#include "TestFluidUtilities.hpp" + +using namespace geos::constitutive; +using namespace geos::constitutive::compositional; + +namespace geos +{ +namespace testing +{ + +static constexpr std::array< real64, 4 > testPressures{1.0, 50.0, 100.0, 500.0}; // Pressure in bars +static constexpr std::array< real64, 4 > testTemperature{10.0, 25.0, 80.0, 150.0}; // Temperature in C + +template< integer NC > +struct FluidData +{}; + +template<> +struct FluidData< 2 > +{ + static constexpr real64 e = 1.0e-6; + + static constexpr std::array< real64, 8 > testCompositions = { + 0.5, 0.5, + 0.0, 1.0, + e, 1.0-e, + 1.0-e, e + }; + + static std::unique_ptr< TestFluid< 2 > > createFluid() + { + auto fluid = TestFluid< 2 >::create( {Fluid::CO2, Fluid::H2O} ); + fluid->setBinaryCoefficients( Feed< 1 >{0.1896} ); + return fluid; + } +}; + +template<> +struct FluidData< 4 > +{ + static constexpr real64 e = 1.0e-6; + + static constexpr std::array< real64, 16 > testCompositions = { + 0.3, 0.2, 0.2, 0.3, + 0.0, 0.0, 0.0, 1.0, + 0.0, 0.5, 0.0, 0.5, + e, e, e, 1.0-3*e + }; + + static std::unique_ptr< TestFluid< 4 > > createFluid() + { + auto fluid = TestFluid< 4 >::create( {Fluid::CH4, Fluid::CO2, Fluid::H2S, Fluid::H2O} ); + const std::array< real64 const, 6 > bics = { 0.0, 0.0, 0.0, 0.4850, 0.1896, 0.1353 }; + fluid->setBinaryCoefficients( bics ); + return fluid; + } +}; + +template<> +struct FluidData< 6 > +{ + static constexpr real64 e = 1.0e-7; + static constexpr real64 f = 0.16666; + + static constexpr std::array< real64, 30 > testCompositions = { + f, f, f, f, f, 1.0-5*f, + e, e, e, e, e, 1.0-5*e, + e, 1.0-5*e, e, e, e, e, + 0.260871, 0.0663935, 0.0906286, 0.0524321, 0.52755293, 0.00212187, + 6.233e-08, 2.99791e-09, 8.00657e-06, 8.187700789e-05, 5.10942e-08, 0.99991 + }; + + static std::unique_ptr< TestFluid< 6 > > createFluid() + { + auto fluid = TestFluid< 6 >::create( {Fluid::CH4, Fluid::C2H6, Fluid::CO2, Fluid::H2S, Fluid::N2, Fluid::H2O} ); + const std::array< real64 const, 15 > bics = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.4850, 0.4920, 0.1896, 0.1353, 0.4778 }; + fluid->setBinaryCoefficients( bics ); + return fluid; + } +}; + +template< integer NC > +using FlashData = std::tuple< + real64 const, // pressure + real64 const, // temperature + Feed< NC > const // total composition + >; + +template< integer NC, EquationOfStateType EOS_TYPE > +class StabilityTestDerivativesTestFixture : public ::testing::TestWithParam< FlashData< NC > > +{ + static constexpr real64 relTol = 1.0e-5; + static constexpr real64 absTol = 1.0e-7; + static constexpr integer numComps = NC; + static constexpr integer numDofs = NC + 2; + + using Deriv = geos::constitutive::multifluid::DerivativeOffset; + +public: + StabilityTestDerivativesTestFixture() + : m_fluid( FluidData< NC >::createFluid() ) + {} + + ~StabilityTestDerivativesTestFixture() = default; + + void testStabilityTestDerivatives( FlashData< NC > const & data ) + { + auto componentProperties = this->m_fluid->createKernelWrapper(); + + // Extract test point + real64 const pressure = std::get< 0 >( data ); + real64 const temperature = std::get< 1 >( data ); + stackArray1d< real64, numComps > composition; + TestFluid< NC >::createArray( composition, std::get< 2 >( data )); + + StackArray< real64, 2, 4*numComps > workSpace( 4, numComps ); + StackArray< real64, 3, 2*numComps *numDofs > derivativeWorkSpace( 2, numComps, numDofs ); + arraySlice1d< real64 > incipientComposition = workSpace[0]; + arraySlice1d< real64 > kValues = workSpace[1]; + arraySlice1d< real64 > derivatives = workSpace[2]; + arraySlice1d< real64 > displacedIncipientComposition = workSpace[3]; + arraySlice2d< real64 > sampleDerivatives = derivativeWorkSpace[0]; + arraySlice2d< real64 > incipientDerivatives = derivativeWorkSpace[1]; + + auto const setZero = []( real64 & value ){ value = 0.0; }; + + LvArray::forValuesInSlice( workSpace.toSlice(), setZero ); + + EquationOfStateType constexpr equationOfState = EOS_TYPE; + + constitutive::compositional::FlashData flashData; + flashData.liquidEos = equationOfState; + if constexpr (equationOfState == EquationOfStateType::SoreideWhitson) + { + flashData.vapourEos = EquationOfStateType::PengRobinson; + } + else + { + flashData.vapourEos = equationOfState; + } + + // Extract column + auto const extractColumn = [&]( integer const idof, real64 const scale = 1.0 ) + { + for( integer ic = 0; ic < numComps; ++ic ) + { + derivatives[ic] = scale * incipientDerivatives( ic, idof ); + } + }; + + if constexpr (equationOfState == EquationOfStateType::SoreideWhitson) + { + KValueInitialization::computeWilsonGasLiquidKvalue( numComps, + pressure, + temperature, + componentProperties, + kValues ); + } + else + { + KValueInitialization::computeSoreideWhitsonKvalue( numComps, + pressure, + temperature, + componentProperties, + composition.toSliceConst(), + kValues ); + } + + auto const parameters = FlashParameters::create( std::make_unique< ModelParameters >() ); + auto const * flashParameters = parameters->get< FlashParameters >(); + flashParameters->m_discreteParameters[FlashParameters::STABILITY_MAX_ITERATIONS] = 20; + + bool unstableMixture = false; + EquationOfStateType outputEos = EOS_TYPE; + bool const stabilityStatus = StabilityTest::compute( numComps, + pressure, + temperature, + composition.toSliceConst(), + componentProperties, + flashData, + kValues.toSliceConst(), + flashParameters->m_continuousParameters, + flashParameters->m_discreteParameters, + unstableMixture, + outputEos, + incipientComposition ); + EXPECT_EQ( stabilityStatus, true ); + + EquationOfStateType const sampleEos = (outputEos == flashData.vapourEos) ? flashData.liquidEos : flashData.vapourEos; + StabilityTest::computeDerivatives( numComps, + pressure, + temperature, + composition.toSliceConst(), + componentProperties, + sampleEos, + outputEos, + flashData, + incipientComposition.toSliceConst(), + incipientDerivatives, + sampleDerivatives ); + + // -- Pressure derivative + extractColumn( Deriv::dP ); + real64 const dp = 1.0e-4 * pressure; + geos::testing::internal::testNumericalDerivative< numComps >( + pressure, dp, derivatives, + [&]( real64 const p, auto & values ){ + bool displacedUnstableMixture = false; + EquationOfStateType displacedOutputEos = EOS_TYPE; + StabilityTest::compute( numComps, + p, + temperature, + composition.toSliceConst(), + componentProperties, + flashData, + kValues.toSliceConst(), + flashParameters->m_continuousParameters, + flashParameters->m_discreteParameters, + displacedUnstableMixture, + displacedOutputEos, + displacedIncipientComposition ); + for( integer ic = 0; ic < numComps; ++ic ) + { + values[ic] = displacedIncipientComposition[ic]; + } + } ); + + // -- Temperature derivative + extractColumn( Deriv::dT ); + real64 const dT = 1.0e-5 * temperature; + geos::testing::internal::testNumericalDerivative< numComps >( + temperature, dT, derivatives, + [&]( real64 const t, auto & values ){ + bool displacedUnstableMixture = false; + EquationOfStateType displacedOutputEos = EOS_TYPE; + StabilityTest::compute( numComps, + pressure, + t, + composition.toSliceConst(), + componentProperties, + flashData, + kValues.toSliceConst(), + flashParameters->m_continuousParameters, + flashParameters->m_discreteParameters, + displacedUnstableMixture, + displacedOutputEos, + displacedIncipientComposition ); + for( integer ic = 0; ic < numComps; ++ic ) + { + values[ic] = displacedIncipientComposition[ic]; + } + } ); + + // -- Composition derivatives + real64 constexpr dz = 1.0e-9; + for( integer kc = 0; kc < numComps; ++kc ) + { + // If the mixture is stable and the component doesn't exist then checking the derivative here will + // not work because the incipient phase composition is not continuous in the component dimension. + if( composition[kc] < MultiFluidConstants::minForSpeciesPresence && !unstableMixture ) + { + continue; + } + extractColumn( Deriv::dC + kc ); + geos::testing::internal::testNumericalDerivative< numComps >( + 0.0, dz, derivatives, + [&]( real64 const z, auto & values ){ + real64 const zOld = composition[kc]; + composition[kc] += z; + bool displacedUnstableMixture = false; + EquationOfStateType displacedOutputEos = EOS_TYPE; + StabilityTest::compute( numComps, + pressure, + temperature, + composition.toSliceConst(), + componentProperties, + flashData, + kValues.toSliceConst(), + flashParameters->m_continuousParameters, + flashParameters->m_discreteParameters, + displacedUnstableMixture, + displacedOutputEos, + displacedIncipientComposition ); + + for( integer ic = 0; ic < numComps; ++ic ) + { + values[ic] = displacedIncipientComposition[ic]; + } + composition[kc] = zOld; + } ); + } + } + +protected: + std::unique_ptr< TestFluid< NC > > m_fluid{}; +}; + +using SoaveRedlichKwong2 = StabilityTestDerivativesTestFixture< 2, EquationOfStateType::SoaveRedlichKwong >; +using PengRobinson4 = StabilityTestDerivativesTestFixture< 4, EquationOfStateType::PengRobinson >; +using SoreideWhitson6 = StabilityTestDerivativesTestFixture< 6, EquationOfStateType::SoreideWhitson >; + +TEST_P( SoaveRedlichKwong2, testDerivatives ) +{ + testStabilityTestDerivatives( GetParam() ); +} + +TEST_P( PengRobinson4, testDerivatives ) +{ + testStabilityTestDerivatives( GetParam() ); +} + +TEST_P( SoreideWhitson6, testDerivatives ) +{ + testStabilityTestDerivatives( GetParam() ); +} + +//------------------------------------------------------------------------------- +// Data +//------------------------------------------------------------------------------- + +template< integer NC > +stdVector< FlashData< NC > > generateData() +{ + constexpr integer dataSize = static_cast< integer >(FluidData< NC >::testCompositions.size()); + integer constexpr numSamples = dataSize / NC; + static_assert( numSamples * NC == dataSize ); + + stdVector< FlashData< NC > > data; + for( integer sampleIndex = 0; sampleIndex < numSamples; sampleIndex++ ) + { + Feed< NC > sample; + for( integer ic = 0; ic < NC; ic++ ) + { + sample[ic] = FluidData< NC >::testCompositions[sampleIndex*NC + ic]; + } + for( real64 const p : testPressures ) // Pressure in bars + { + for( real64 const t : testTemperature ) // Temperature in C + { + data.emplace_back( 1.0e5*p, t + 273.15, sample ); + } + } + } + return data; +} + +INSTANTIATE_TEST_SUITE_P( StabilityTestDerivatives, SoaveRedlichKwong2, ::testing::ValuesIn( generateData< 2 >()) ); +INSTANTIATE_TEST_SUITE_P( StabilityTestDerivatives, PengRobinson4, ::testing::ValuesIn( generateData< 4 >()) ); +INSTANTIATE_TEST_SUITE_P( StabilityTestDerivatives, SoreideWhitson6, ::testing::ValuesIn( generateData< 6 >()) ); + +} // testing +} // geos diff --git a/src/coreComponents/constitutive/unitTests/testStabilityTestJacobian.cpp b/src/coreComponents/constitutive/unitTests/testStabilityTestJacobian.cpp new file mode 100644 index 00000000000..fa86274cc94 --- /dev/null +++ b/src/coreComponents/constitutive/unitTests/testStabilityTestJacobian.cpp @@ -0,0 +1,290 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +// Source includes +#include "codingUtilities/UnitTestUtilities.hpp" +#include "constitutive/fluid/multifluid/compositional/functions/StabilityTest.hpp" +#include "constitutive/fluid/multifluid/compositional/functions/KValueInitialization.hpp" +#include "constitutive/fluid/multifluid/compositional/functions/Utilities.hpp" +#include "constitutive/fluid/multifluid/MultiFluidConstants.hpp" +#include "TestFluid.hpp" +#include "TestFluidUtilities.hpp" + +using namespace geos::constitutive; +using namespace geos::constitutive::compositional; + +namespace geos +{ +namespace testing +{ + +template< integer NC > +struct FluidData +{}; + +template<> +struct FluidData< 2 > +{ + static constexpr real64 e = 1.0e-8; + + static constexpr std::array< real64, 8 > testCompositions = { + 0.5, 0.5, + 0.0, 1.0, + e, 1.0-e, + 1.0-e, e + }; + + static std::unique_ptr< TestFluid< 2 > > createFluid() + { + auto fluid = TestFluid< 2 >::create( {Fluid::CO2, Fluid::C5H12} ); + fluid->setBinaryCoefficients( Feed< 1 >{0.025} ); + return fluid; + } +}; + +template<> +struct FluidData< 4 > +{ + static constexpr real64 e = 1.0e-7; + + static constexpr std::array< real64, 16 > testCompositions = { + 0.3, 0.2, 0.2, 0.3, + 0.0, 0.0, 0.0, 1.0, + 0.0, 0.5, 0.0, 0.5, + e, e, 0.0, 1.0-2*e + }; + + static std::unique_ptr< TestFluid< 4 > > createFluid() + { + auto fluid = TestFluid< 4 >::create( {Fluid::CH4, Fluid::CO2, Fluid::H2S, Fluid::H2O} ); + const std::array< real64 const, 6 > bics = { 0.0, 0.0, 0.0, 0.4850, 0.1896, 0.1353 }; + fluid->setBinaryCoefficients( bics ); + return fluid; + } +}; + +template< integer NC > +using FlashData = std::tuple< + real64 const, // pressure + real64 const, // temperature + Feed< NC > const, // total composition + real64 const // power to raise k-values + >; + +template< integer NC, EquationOfStateType EOS_TYPE > +class StabilityTestJacobianTestFixture : public ::testing::TestWithParam< FlashData< NC > > +{ + static constexpr real64 relTol = 1.0e-5; + static constexpr real64 absTol = 1.0e-7; + static constexpr integer numComps = NC; + static constexpr integer numDofs = NC + 2; + using Deriv = geos::constitutive::multifluid::DerivativeOffset; +public: + StabilityTestJacobianTestFixture() + : m_fluid( FluidData< NC >::createFluid() ) + {} + + ~StabilityTestJacobianTestFixture() = default; + + void testJacobian( FlashData< NC > const & data ) + { + auto componentProperties = this->m_fluid->createKernelWrapper(); + + EquationOfStateType constexpr equationOfState = EOS_TYPE; + + constitutive::compositional::FlashData flashData; + flashData.liquidEos = equationOfState; + if constexpr (equationOfState == EquationOfStateType::SoreideWhitson) + { + flashData.vapourEos = EquationOfStateType::PengRobinson; + } + else + { + flashData.vapourEos = equationOfState; + } + + real64 const pressure = std::get< 0 >( data ); + real64 const temperature = std::get< 1 >( data ); + stackArray1d< real64, numComps > composition; + TestFluid< NC >::createArray( composition, std::get< 2 >( data )); + real64 const alpha = std::get< 3 >( data ); + + StackArray< real64, 2, 7*numComps > workSpace( 7, numComps ); + StackArray< real64, 3, numComps *numDofs > derivativeWorkSpace( 1, numComps, numDofs ); + StackArray< real64, 3, numComps *numComps > jacobianSpace( 1, numComps, numComps ); + StackArray< integer, 1, numComps > nonZeroComponents( numComps ); + arraySlice1d< real64 > testComposition = workSpace[0]; + arraySlice1d< real64 > logTestComposition = workSpace[1]; + arraySlice1d< real64 > testFugacity = workSpace[2]; + arraySlice1d< real64 > kValues = workSpace[3]; + arraySlice2d< real64 > logTestFugacityDerivs = derivativeWorkSpace[0]; + arraySlice2d< real64 > jacobian = jacobianSpace[0]; + arraySlice1d< real64 > residual = workSpace[4]; + arraySlice1d< real64 > displacedResidual = workSpace[5]; + arraySlice1d< real64 > derivatives = workSpace[6]; + + auto const setZero = []( real64 & value ){ value = 0.0; }; + + LvArray::forValuesInSlice( workSpace.toSlice(), setZero ); + + // Extract column + auto const extractColumn = [&]( integer const kc, real64 const scale = 1.0 ) + { + for( integer ic = 0; ic < numComps; ++ic ) + { + derivatives[ic] = scale * jacobian( ic, kc ); + } + }; + + calculatePresentComponents( numComps, composition.toSliceConst(), nonZeroComponents ); + auto presentComponents = nonZeroComponents.toSliceConst(); + + if constexpr (equationOfState == EquationOfStateType::SoreideWhitson) + { + KValueInitialization::computeWilsonGasLiquidKvalue( numComps, + pressure, + temperature, + componentProperties, + kValues ); + } + else + { + KValueInitialization::computeSoreideWhitsonKvalue( numComps, + pressure, + temperature, + componentProperties, + composition.toSliceConst(), + kValues ); + } + + for( integer const ic: presentComponents ) + { + real64 const logK = LvArray::math::log( kValues[ic] ); + logTestComposition[ic] = LvArray::math::log( composition[ic] ) + alpha * logK; + } + + real64 tpd = 0.0; + + StabilityTest::calculateResidualAndJacobian( numComps, + presentComponents, + pressure, + temperature, + componentProperties, + equationOfState, + flashData, + logTestComposition.toSliceConst(), + tpd, + testComposition, + testFugacity, + logTestFugacityDerivs, + residual, + jacobian ); + + // composition derivatives + // Derivatives are wrt to log(Y) + real64 constexpr deltaLogY = 1.0e-3; + for( integer kc = 0; kc < numComps; ++kc ) + { + extractColumn( kc ); + geos::testing::internal::testNumericalDerivative< numComps >( + 0.0, deltaLogY, derivatives, + [&]( real64 const dY, auto & values ) + { + LvArray::forValuesInSlice( displacedResidual, setZero ); + real64 const logYOld = logTestComposition[kc]; + logTestComposition[kc] += dY; + StabilityTest::calculateResidual( numComps, + presentComponents, + pressure, + temperature, + componentProperties, + equationOfState, + flashData, + logTestComposition.toSliceConst(), + tpd, + testComposition, + testFugacity, + displacedResidual ); + for( integer ic = 0; ic < numComps; ++ic ) + { + values[ic] = displacedResidual[ic]; + } + logTestComposition[kc] = logYOld; + } ); + } + } + +protected: + std::unique_ptr< TestFluid< NC > > m_fluid{}; +}; + +using PengRobinson4 = StabilityTestJacobianTestFixture< 4, EquationOfStateType::PengRobinson >; +using SoreideWhitson4 = StabilityTestJacobianTestFixture< 4, EquationOfStateType::SoreideWhitson >; +using SoaveRedlichKwong2 = StabilityTestJacobianTestFixture< 2, EquationOfStateType::SoaveRedlichKwong >; + +TEST_P( PengRobinson4, testJacobian ) +{ + testJacobian( GetParam() ); +} + +TEST_P( SoreideWhitson4, testJacobian ) +{ + testJacobian( GetParam() ); +} + +TEST_P( SoaveRedlichKwong2, testJacobian ) +{ + testJacobian( GetParam() ); +} + +//------------------------------------------------------------------------------- +// Data +//------------------------------------------------------------------------------- + +template< integer NC > +stdVector< FlashData< NC > > generateData() +{ + constexpr integer dataSize = static_cast< integer >(FluidData< NC >::testCompositions.size()); + integer constexpr numSamples = dataSize / NC; + static_assert( numSamples * NC == dataSize ); + + stdVector< FlashData< NC > > data; + for( integer sampleIndex = 0; sampleIndex < numSamples; sampleIndex++ ) + { + Feed< NC > sample; + for( integer ic = 0; ic < NC; ic++ ) + { + sample[ic] = FluidData< NC >::testCompositions[sampleIndex*NC + ic]; + } + for( real64 const p : {1.0, 50.0, 100.0, 500.0} ) // Pressure in bars + { + for( real64 const t : {10.0, 25.0, 80.0, 150.0} ) // Temperature in C + { + for( real64 const a : {-1.0, 0.0, 1.0} ) + { + data.emplace_back( 1.0e5*p, t + 273.15, sample, a ); + } + } + } + } + return data; +} + +INSTANTIATE_TEST_SUITE_P( StabilityTestJacobian, PengRobinson4, ::testing::ValuesIn( generateData< 4 >()) ); +INSTANTIATE_TEST_SUITE_P( StabilityTestJacobian, SoreideWhitson4, ::testing::ValuesIn( generateData< 4 >()) ); +INSTANTIATE_TEST_SUITE_P( StabilityTestJacobian, SoaveRedlichKwong2, ::testing::ValuesIn( generateData< 2 >()) ); + +} // testing +} // geos diff --git a/src/coreComponents/schema/schema.xsd b/src/coreComponents/schema/schema.xsd index 5be40a17911..6e68e66858f 100644 --- a/src/coreComponents/schema/schema.xsd +++ b/src/coreComponents/schema/schema.xsd @@ -7493,6 +7493,8 @@ if the table is requested to be output in the log, and it is too large, a CSV fi + + @@ -7555,6 +7557,8 @@ if the table is requested to be output in the log, and it is too large, a CSV fi + + @@ -7594,6 +7598,8 @@ if the table is requested to be output in the log, and it is too large, a CSV fi + + @@ -7644,6 +7650,8 @@ if the table is requested to be output in the log, and it is too large, a CSV fi + +