Remove penalization in fugacity computations and improve Gibbs energy calculation

.dependencies/python.devenv.yml

@@ -3,6 +3,7 @@
 dependencies:
   - numpy
   - pandas
+  - matplotlib
   - pip
   - tabulate
   - colorama

Reaktoro/Core/ChemicalProperties.cpp

@@ -321,6 +321,13 @@ auto ChemicalProperties::phaseSpecificVolumes() const -> ChemicalVector
     return phaseAmounts()/phaseMasses() % phaseMolarVolumes();
 }
 
+auto ChemicalProperties::phaseCompressibilityFactors() const -> ChemicalVector
+{
+    const auto& R = universalGasConstant;
+    const auto& v = phaseMolarVolumes();
+    return P * v / (R * T);
+}
+
 auto ChemicalProperties::phaseSpecificEntropies() const -> ChemicalVector
 {
     return phaseAmounts()/phaseMasses() % phaseMolarEntropies();

Reaktoro/Core/ChemicalProperties.hpp

@@ -151,6 +151,9 @@ class ChemicalProperties
     /// Return the specific volumes of the phases (in units of m3/kg).
     auto phaseSpecificVolumes() const -> ChemicalVector;
 
+    /// Return the compressibility factor of the phases.
+    auto phaseCompressibilityFactors() const -> ChemicalVector;
+
     /// Return the specific entropies of the phases (in units of J/(kg*K)).
     auto phaseSpecificEntropies() const -> ChemicalVector;
 

Reaktoro/Thermodynamics/EOS/CubicEOS.cpp

@@ -156,6 +156,131 @@ auto Psi(CubicEOS::Model type) -> double
 
 } // namespace internal
 
+auto computeSpeciesFugacity(
+    const ThermoScalar& P,
+    const ThermoScalar& T,
+    const ChemicalScalar& xi,
+    const ChemicalScalar& ai,
+    const double& bi,
+    const ChemicalScalar& aiT,
+    const ChemicalScalar& amix,
+    const ChemicalScalar& amixT,
+    const ChemicalScalar& bmix,
+    const ChemicalScalar& A,
+    const ChemicalScalar& B,
+    const ChemicalScalar& C,
+    const ChemicalScalar& Z,
+    const double epsilon,
+    const double sigma) -> ChemicalScalar
+{
+    const double R = universalGasConstant;
+
+    const double almost_zero = 1e-20;
+
+    const ChemicalScalar q = amix/(bmix*R*T);
+    const ChemicalScalar qT = q*(amixT/amix - 1.0/T);
+
+    const ChemicalScalar beta = P*bmix/(R*T);
+    const ThermoScalar betai = P*bi/(R*T);
+
+    const ChemicalScalar qi = q*(1 + ai/amix - bi/bmix);
+    const ChemicalScalar qiT = qi*qT/q + q*(aiT - ai*amixT/amix)/amix;
+
+    const ThermoScalar Ai = (epsilon + sigma - 1.0)*betai - 1.0;
+    const ChemicalScalar Bi = (epsilon*sigma - epsilon - sigma)*(2*beta*betai - beta*beta) - (epsilon + sigma - q)*(betai - beta) - (epsilon + sigma - qi)*beta;
+    const ChemicalScalar Ci = -3*sigma*epsilon*beta*beta*betai + 2*epsilon*sigma*beta*beta*beta - (epsilon*sigma + qi)*beta*beta - 2*(epsilon*sigma + q)*(beta*betai - beta*beta);
+    const ChemicalScalar Zi = -(Ai*Z*Z + (Bi + B)*Z + Ci + 2*C)/(3*Z*Z + 2*A*Z + B);
+
+    // Calculate the integration factor I
+    ChemicalScalar I;
+    if(std::abs(epsilon - sigma) > almost_zero) I = log((Z + sigma*beta)/(Z + epsilon*beta))/(sigma - epsilon);
+                    else I = beta/(Z + epsilon*beta);
+
+    // Calculate the integration factor I for each component
+    ChemicalScalar Ii;
+    if(std::abs(epsilon - sigma) > almost_zero) Ii = I + ((Zi + sigma*betai)/(Z + sigma*beta) - (Zi + epsilon*betai)/(Z + epsilon*beta))/(sigma - epsilon);
+                    else Ii = I * (1 + betai/beta - (Zi + epsilon*betai)/(Z + epsilon*beta));
+
+    auto Gi_res = R*T*(Zi - (Zi - betai)/(Z - beta) - log(Z - beta) - qi*I - q*Ii + q*I);
+    auto ln_fi = Gi_res/(R*T) + log(xi * P);
+    auto fi = exp(ln_fi);
+    return fi;
+}
+
+auto calculateNormalizedPhaseGibbsEnergy(
+    const unsigned& nspecies,
+    const ThermoScalar& P,
+    const ThermoScalar& T,
+    const ChemicalVector& x,
+    const ChemicalVector& abar,
+    const Vector& bbar,
+    const ChemicalVector& abarT,
+    const ChemicalScalar& amix,
+    const ChemicalScalar& amixT,
+    const ChemicalScalar& bmix,
+    const ChemicalScalar& A,
+    const ChemicalScalar& B,
+    const ChemicalScalar& C,
+    const ChemicalScalar& Z,
+    const double epsilon,
+    const double sigma) -> ChemicalScalar
+{
+    ChemicalScalar G_normalized;
+    for(unsigned i = 0; i < nspecies; ++i)
+    {
+        auto xi = x[i];
+        auto fi = computeSpeciesFugacity(P, T, xi, abar[i], bbar[i], abarT[i], amix, amixT, bmix, A, B, C, Z, epsilon, sigma);
+        G_normalized += xi * log(fi);
+    }
+    return G_normalized;
+}
+
+auto selectCompressibilityFactorByGibbsEnergy(
+    const unsigned& nspecies,
+    std::vector<ChemicalScalar> Zs,
+    const ThermoScalar& P,
+    const ThermoScalar& T,
+    const ChemicalVector& x,
+    const ChemicalVector& abar,
+    const Vector& bbar,
+    const ChemicalVector& abarT,
+    const ChemicalScalar& amix,
+    const ChemicalScalar& amixT,
+    const ChemicalScalar& bmix,
+    const ChemicalScalar& A,
+    const ChemicalScalar& B,
+    const ChemicalScalar& C,
+    const double epsilon,
+    const double sigma) -> ChemicalScalar
+{
+    ChemicalScalar Z(nspecies);
+    if (Zs.size() == 1)
+    {
+        Z = Zs[0];
+        return Z;
+    }
+
+    if (Zs.size() != 2) {
+        Exception exception;
+        exception.error << "selectCompressibilityByGibbsEnergy received invalid input";
+        exception.reason << "Zs should have size 1 or 2 in selectCompressibilityByGibbsEnergy, "
+            << "but has a size of " << Zs.size();
+        RaiseError(exception);
+    }
+
+    std::vector<ChemicalScalar> normalized_Gs;
+    normalized_Gs.push_back(
+        calculateNormalizedPhaseGibbsEnergy(nspecies, P, T, x, abar, bbar, abarT, amix, amixT, bmix, A, B, C, Zs[0], epsilon, sigma)
+    );
+    normalized_Gs.push_back(
+        calculateNormalizedPhaseGibbsEnergy(nspecies, P, T, x, abar, bbar, abarT, amix, amixT, bmix, A, B, C, Zs[1], epsilon, sigma)
+    );
+
+    Z = normalized_Gs[0] < normalized_Gs[1] ? Zs[0] : Zs[1];
+
+    return Z;
+}
+
 struct CubicEOS::Impl
 {
     /// The number of species in the phase.
@@ -340,11 +465,15 @@ struct CubicEOS::Impl
         const auto cubic_size = cubicEOS_roots.size();
 
         std::vector<ChemicalScalar> Zs;
-        if (cubic_size == 1 || cubic_size == 2)
+        if (cubic_size == 1)
         {
-            //even if cubicEOS_roots has 2 roots, assume that the smallest does not have physical meaning
             Zs.push_back(ChemicalScalar(nspecies, cubicEOS_roots[0]));
         }
+        else if (cubic_size == 2)
+        {
+            Zs.push_back(ChemicalScalar(nspecies, cubicEOS_roots[0]));
+            Zs.push_back(ChemicalScalar(nspecies, cubicEOS_roots[1]));
+        }
         else
         {
             if (cubic_size != 3) {
@@ -357,16 +486,18 @@ struct CubicEOS::Impl
             Zs.push_back(ChemicalScalar(nspecies, cubicEOS_roots[2]));  // Z_min
         }
 
-        // Selecting compressibility factor - Z_liq < Z_gas
         ChemicalScalar Z(nspecies);
+
+        // Selecting compressibility factor - Z_liq < Z_gas
+        //TODO: study the possibilty to use selectCompressibilityFactorByGibbsEnergy for Z selection
         if (isvapor)
             Z.val = *std::max_element(cubicEOS_roots.begin(), cubicEOS_roots.end());
         else
             Z.val = *std::min_element(cubicEOS_roots.begin(), cubicEOS_roots.end());
 
         auto input_phase_type = isvapor ? PhaseType::Gas : PhaseType::Liquid;
-        auto identified_phase_type = input_phase_type;
 
+        auto identified_phase_type = input_phase_type;
         switch (phase_identification_method)
         {
         case PhaseIdentificationMethod::None:
@@ -390,8 +521,7 @@ struct CubicEOS::Impl
             throw std::logic_error("CubicEOS received an unexpected phaseIdentificationMethod");
         }
 
-        if (identified_phase_type != input_phase_type)
-        {
+        if (identified_phase_type != input_phase_type) {
             // Since the phase is identified as different than the expect input phase type, it is
             // deemed inappropriate. Artificially high values are configured for fugacities, so that
             // this condition is "removed" by the optimizer.
@@ -403,7 +533,7 @@ struct CubicEOS::Impl
             result.partial_molar_volumes.fill(0.0);
             result.residual_partial_molar_gibbs_energies.fill(0.0);
             result.residual_partial_molar_enthalpies.fill(0.0);
-            result.ln_fugacity_coefficients.fill(100.0);
+            result.ln_fugacity_coefficients.fill(0.0);
             return result;
         }
 
@@ -493,7 +623,7 @@ CubicEOS::Result::Result(unsigned nspecies)
 {}
 
 /// Sanity check free function to verify if BIPs matrices have proper dimensions. Considering that the phase has
-/// n species, the BIP matricies k, kT and kTT should have (n, n) as dimensions.
+/// n species, the BIP matrices k, kT and kTT should have (n, n) as dimensions.
 /// @see CubicEOS::setInteractionParamsFunction
 auto sanityCheckInteractionParamsFunction(const unsigned& nspecies, const CubicEOS::InteractionParamsFunction& func) -> void
 {

Reaktoro/Thermodynamics/EOS/CubicEOS.hpp

@@ -50,10 +50,13 @@ class CubicEOS
     /// Note that the BIPs can depend on the temperature, thus kT and kTT should be also provided.
     struct InteractionParamsResult
     {
+        /// The BIPs matrix. The size must be (n, n), where n is the number of species
         MatrixXd k;
 
+        /// The derivative of each k entry w.r.t T.
         MatrixXd kT;
 
+        /// The derivative of each kT entry w.r.t T.
         MatrixXd kTT;
     };
 

Reaktoro/Thermodynamics/EOS/PhaseIdentification.cpp

@@ -107,6 +107,38 @@ auto pressureComparison(const ThermoScalar& Pressure, const ThermoScalar& Temper
     RaiseError(exception);
 }
 
+auto computeGibbsResidualEnergy(
+    const ThermoScalar& pressure,
+    const ThermoScalar& temperature,
+    const ChemicalScalar& amix,
+    const ChemicalScalar& bmix,
+    const ChemicalScalar& A,
+    const ChemicalScalar& B,
+    const ChemicalScalar& Z,
+    const double epsilon,
+    const double sigma) -> ChemicalScalar
+{
+    auto const almost_zero = 1e-20;
+    auto constexpr R = universalGasConstant;
+    auto const& T = temperature;
+
+    // Computing the values of residual Gibbs energy for all Zs
+    ChemicalScalar Gs;
+    const double factor = -1.0 / (3 * Z.val*Z.val + 2 * A.val*Z.val + B.val);
+    const ChemicalScalar beta = pressure * bmix / (R * T);
+    const ChemicalScalar q = amix / (bmix * R * T);
+
+    // Calculate the integration factor I
+    ChemicalScalar I;
+    if (std::abs(epsilon - sigma) > almost_zero)
+        I = log((Z + sigma * beta) / (Z + epsilon * beta)) / (sigma - epsilon);
+    else 
+        I = beta / (Z + epsilon * beta);        
+
+    Gs = R * temperature*(Z - 1 - log(Z - beta) - q * I);
+
+    return Gs;
+}
 
 auto gibbsResidualEnergyComparison(
     const ThermoScalar& pressure,
@@ -127,24 +159,24 @@ auto gibbsResidualEnergyComparison(
     std::vector<ChemicalScalar> Gs;
     for (const auto Z : {Z_max, Z_min})
     {
-        const double factor = -1.0 / (3 * Z.val*Z.val + 2 * A.val*Z.val + B.val);
-        const ChemicalScalar beta = pressure * bmix / (R * T);
-        const ChemicalScalar q = amix / (bmix * R * T);
-
-        // Calculate the integration factor I and its temperature derivative IT
-        ChemicalScalar I;
-        if (epsilon != sigma) 
-            I = log((Z + sigma * beta) / (Z + epsilon * beta)) / (sigma - epsilon);
-        else 
-            I = beta / (Z + epsilon * beta);        
-
-        Gs.push_back(R * temperature*(Z - 1 - log(Z - beta) - q * I));
+        auto current_G = computeGibbsResidualEnergy(
+            pressure,
+            temperature,
+            amix,
+            bmix,
+            A,
+            B,
+            Z,
+            epsilon,
+            sigma
+        );
+
+        Gs.push_back(current_G);
     }
 
     return (Gs[0].val < Gs[1].val) ? PhaseType::Gas : PhaseType::Liquid;
 }
 
-
 auto identifyPhaseUsingGibbsEnergyAndEos(
     const ThermoScalar& pressure,
     const ThermoScalar& temperature,

Reaktoro/Thermodynamics/EOS/PhaseIdentification.hpp

@@ -65,7 +65,7 @@ auto identifyPhaseUsingIsothermalCompressibility(
 /// @param temperature Phase temperature
 /// @param amix attractive parameter
 /// @param Zs Z-roots, as calculated by the cubic EOS. Must have size 1 or 2 here.
-///     If size(Z) == 2, the values od Gibbs residual energy are compared. It is a liquid phase if
+///     If size(Z) == 2, the values of Gibbs residual energy are compared. It is a liquid phase if
 ///     Gibbs residual energy of Z_min is the smallest and gaseous if Gibbs residual energy of Z_max
 ///     is the smallest.
 ///     If size(Z) == 1, the pressue is compared with the local P_min and local P_max of the EoS.

python/PyReaktoro/Core/PyChemicalProperties.cpp

@@ -67,6 +67,7 @@ void exportChemicalProperties(py::module& m)
         .def("phaseSpecificGibbsEnergies", &ChemicalProperties::phaseSpecificGibbsEnergies)
         .def("phaseSpecificEnthalpies", &ChemicalProperties::phaseSpecificEnthalpies)
         .def("phaseSpecificVolumes", &ChemicalProperties::phaseSpecificVolumes)
+        .def("phaseCompressibilityFactors", &ChemicalProperties::phaseCompressibilityFactors)
         .def("phaseSpecificEntropies", &ChemicalProperties::phaseSpecificEntropies)
         .def("phaseSpecificInternalEnergies", &ChemicalProperties::phaseSpecificInternalEnergies)
         .def("phaseSpecificHelmholtzEnergies", &ChemicalProperties::phaseSpecificHelmholtzEnergies)

python/reaktoro/CMakeLists.txt

@@ -21,7 +21,7 @@ add_custom_target(reaktoro ALL
     COMMAND ${PYTHON_EXECUTABLE} -m pip install
         reaktoro --prefix=${CMAKE_BINARY_DIR}
         --find-links=${CMAKE_CURRENT_BINARY_DIR}/dist
-        --no-index --no-deps --no-cache-dir
+        --no-index --no-deps --no-cache-dir --force-reinstall
     WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
 
 # Set dependencies of reaktoro target
@@ -66,7 +66,7 @@ install(CODE
         COMMAND ${PYTHON_EXECUTABLE} -m pip
             install reaktoro --prefix=\${REAKTORO_PYTHON_INSTALL_PREFIX_NATIVE}
             --find-links=\${REAKTORO_PYTHON_DIST_DIR_NATIVE}
-	    --no-index --no-deps --no-cache-dir
+	    --no-index --no-deps --no-cache-dir --force-reinstall
         WORKING_DIRECTORY ${REAKTORO_PYTHON_INSTALL_PREFIX})
 "
 )