VLERoutines.h

#ifndef VLEROUTINES_H
#define VLEROUTINES_H

#include "HelmholtzEOSMixtureBackend.h"

namespace CoolProp {

namespace SaturationSolvers {
struct saturation_T_pure_Akasaka_options
{
    bool use_guesses;  
    CoolPropDbl omega, rhoL, rhoV, pL, pV;
    saturation_T_pure_Akasaka_options(bool use_guesses = false)
      : use_guesses(use_guesses), omega(_HUGE), rhoL(_HUGE), rhoV(_HUGE), pL(_HUGE), pV(_HUGE) {}
};
struct saturation_T_pure_options
{
    CoolPropDbl omega, rhoL, rhoV, pL, pV, p, T;
    saturation_T_pure_options() : omega(_HUGE), rhoL(_HUGE), rhoV(_HUGE), pL(_HUGE), pV(_HUGE), p(_HUGE), T(_HUGE) {}
};

struct saturation_D_pure_options
{
    enum imposed_rho_options
    {
        IMPOSED_RHOL,
        IMPOSED_RHOV
    };
    bool use_guesses,  
      use_logdelta;    
    CoolPropDbl omega, rhoL, rhoV, pL, pV;
    int imposed_rho;
    int max_iterations;
    saturation_D_pure_options() : use_guesses(false), rhoL(_HUGE), rhoV(_HUGE), pL(_HUGE), pV(_HUGE), imposed_rho(0), max_iterations(200) {
        use_logdelta = true;
        omega = 1.0;
    }  // Defaults
};

enum sstype_enum
{
    imposed_T,
    imposed_p
};
struct mixture_VLE_IO
{
    sstype_enum sstype;
    int Nstep_max;
    CoolPropDbl rhomolar_liq, rhomolar_vap, p, T, beta;
    std::vector<CoolPropDbl> x, y, K;
};

static CoolPropDbl Wilson_lnK_factor(const HelmholtzEOSMixtureBackend& HEOS, CoolPropDbl T, CoolPropDbl p, std::size_t i) {
    const double pci = HEOS.get_fluid_constant(i, iP_critical);
    const double Tci = HEOS.get_fluid_constant(i, iT_critical);
    const double omegai = HEOS.get_fluid_constant(i, iacentric_factor);
    return log(pci / p) + 5.373 * (1 + omegai) * (1 - Tci / T);
};

void saturation_D_pure(HelmholtzEOSMixtureBackend& HEOS, CoolPropDbl rhomolar, saturation_D_pure_options& options);
void saturation_T_pure(HelmholtzEOSMixtureBackend& HEOS, CoolPropDbl T, saturation_T_pure_options& options);
void saturation_T_pure_Akasaka(HelmholtzEOSMixtureBackend& HEOS, CoolPropDbl T, saturation_T_pure_Akasaka_options& options);
void saturation_T_pure_Maxwell(HelmholtzEOSMixtureBackend& HEOS, CoolPropDbl T, saturation_T_pure_Akasaka_options& options);

struct saturation_PHSU_pure_options
{
    enum specified_variable_options
    {
        IMPOSED_HL,
        IMPOSED_HV,
        IMPOSED_PL,
        IMPOSED_PV,
        IMPOSED_SL,
        IMPOSED_SV,
        IMPOSED_UL,
        IMPOSED_UV,
        IMPOSED_INVALID_INPUT
    };
    bool use_guesses,  
      use_logdelta;    
    specified_variable_options specified_variable;
    CoolPropDbl omega, rhoL, rhoV, pL, pV, T, p;
    saturation_PHSU_pure_options() : use_logdelta(true), rhoL(_HUGE), rhoV(_HUGE), pL(_HUGE), pV(_HUGE), T(_HUGE), p(_HUGE) {
        specified_variable = IMPOSED_INVALID_INPUT;
        use_guesses = true;
        omega = 1.0;
    }
};

#if !defined(NO_FMTLIB) && FMT_VERSION >= 90000
inline int format_as(saturation_PHSU_pure_options::specified_variable_options option) {
    return fmt::underlying(option);
}
#endif

void saturation_PHSU_pure(HelmholtzEOSMixtureBackend& HEOS, CoolPropDbl specified_value, saturation_PHSU_pure_options& options);

/* \brief This is a backup saturation_p solver for the case where the Newton solver cannot approach closely enough the solution
     *
     * This is especially a problem at low pressures where catastrophic truncation error occurs, especially in the saturated vapor side
     *
     * @param HEOS The Helmholtz EOS backend instance to be used
     * @param p Imposed pressure in kPa
     * @param options Options to be passed to the function (at least T, rhoL and rhoV must be provided)
     */
void saturation_P_pure_1D_T(HelmholtzEOSMixtureBackend& HEOS, CoolPropDbl p, saturation_PHSU_pure_options& options);

/* \brief This is a backup saturation_T solver for the case where the Newton solver cannot approach closely enough the solution
     *
     * This is especially a problem at low pressures where catastrophic truncation error occurs, especially in the saturated vapor side
     *
     * @param HEOS The Helmholtz EOS backend instance to be used
     * @param T Imposed temperature in K
     * @param options Options to be passed to the function (at least p, rhoL and rhoV must be provided)
     */
void saturation_T_pure_1D_P(HelmholtzEOSMixtureBackend& HEOS, CoolPropDbl T, saturation_T_pure_options& options);

/* \brief A robust but slow solver in the very-near-critical region
     *
     * This solver operates in the following fashion:
     * 1. Using a bounded interval for rho'':[rhoc, rhoc-??], guess a value for rho''
     * 2. For guessed value of rho'' and given value of T, calculate p
     * 3. Using a Brent solver on the other co-existing phase (rho'), calculate the (bounded) value of rho' that yields the same pressure
     * 4. Use another outer Brent solver on rho'' to enforce the same Gibbs function between liquid and vapor
     * 5. Fin.
     *
     * @param HEOS The Helmholtz EOS backend instance to be used
     * @param ykey The CoolProp::parameters key to be imposed - one of iT or iP
     * @param y The value for the imposed variable
     */
void saturation_critical(HelmholtzEOSMixtureBackend& HEOS, CoolProp::parameters ykey, CoolPropDbl y);

void successive_substitution(HelmholtzEOSMixtureBackend& HEOS, const CoolPropDbl beta, CoolPropDbl T, CoolPropDbl p,
                             const std::vector<CoolPropDbl>& z, std::vector<CoolPropDbl>& K, mixture_VLE_IO& options);
void x_and_y_from_K(CoolPropDbl beta, const std::vector<CoolPropDbl>& K, const std::vector<CoolPropDbl>& z, std::vector<CoolPropDbl>& x,
                    std::vector<CoolPropDbl>& y);

class WilsonK_resid : public FuncWrapper1D
{
   public:
    sstype_enum input_type;
    double T, p, beta;
    const std::vector<CoolPropDbl>& z;
    std::vector<CoolPropDbl>& K;
    const HelmholtzEOSMixtureBackend& HEOS;

    WilsonK_resid(const HelmholtzEOSMixtureBackend& HEOS, double beta, double imposed_value, sstype_enum input_type,
                  const std::vector<CoolPropDbl>& z, std::vector<CoolPropDbl>& K)
      : input_type(input_type),
        T(imposed_value),
        p(imposed_value),
        beta(beta),
        z(z),
        K(K),
        HEOS(HEOS) {}  // if input_type == imposed_T -> use T, else use p; init both
    double call(double input_value) {
        double summer = 0;
        if (input_type == imposed_T) {
            p = input_value;  // Iterate on pressure
        } else {
            T = input_value;  // Iterate on temperature, pressure imposed
        }
        for (unsigned int i = 0; i < z.size(); i++) {
            K[i] = exp(Wilson_lnK_factor(HEOS, T, p, i));
            summer += z[i] * (K[i] - 1) / (1 - beta + beta * K[i]);
        }
        return summer;
    }
};
inline double saturation_preconditioner(HelmholtzEOSMixtureBackend& HEOS, double input_value, sstype_enum input_type,
                                        const std::vector<CoolPropDbl>& z) {
    double ptriple = 0, pcrit = 0, Ttriple = 0, Tcrit = 0;

    if (HEOS.get_components().empty()) {
        return -1;
    }

    for (unsigned int i = 0; i < z.size(); i++) {
        Tcrit += HEOS.get_fluid_constant(i, iT_critical) * z[i];
        pcrit += HEOS.get_fluid_constant(i, iP_critical) * z[i];
        Ttriple += HEOS.get_fluid_constant(i, iT_triple) * z[i];
        ptriple += HEOS.get_fluid_constant(i, iP_triple) * z[i];
    }
    // Return an invalid number if either triple point temperature or pressure are not available
    if (!ValidNumber(Ttriple) || !ValidNumber(ptriple)) {
        return _HUGE;
    }

    if (input_type == imposed_T) {
        return exp(log(pcrit / ptriple) / (Tcrit - Ttriple) * (input_value - Ttriple) + log(ptriple));
    } else if (input_type == imposed_p) {
        return 1 / (1 / Tcrit - (1 / Ttriple - 1 / Tcrit) / log(pcrit / ptriple) * log(input_value / pcrit));
    } else {
        throw ValueError();
    }
}
inline double saturation_Wilson(HelmholtzEOSMixtureBackend& HEOS, double beta, double input_value, sstype_enum input_type,
                                const std::vector<CoolPropDbl>& z, double guess) {
    double out = 0;
    std::string errstr;

    // If T is input and beta = 0 or beta = 1, explicit solution for p is possible
    if (input_type == imposed_T && (std::abs(beta) < 1e-12 || std::abs(beta - 1) < 1e-12)) {
        const std::vector<double> z = HEOS.get_mole_fractions_ref();
        bool beta0 = std::abs(beta) < 1e-12;  // True is beta is approx. zero
        for (int i = 0; i < static_cast<int>(z.size()); ++i) {
            double pci = HEOS.get_fluid_constant(i, iP_critical);
            double Tci = HEOS.get_fluid_constant(i, iT_critical);
            double omegai = HEOS.get_fluid_constant(i, iacentric_factor);
            if (beta0) {
                out += z[i] * pci * exp(5.373 * (1 + omegai) * (1 - Tci / input_value));
            } else {
                out += z[i] / (pci * exp(5.373 * (1 + omegai) * (1 - Tci / input_value)));
            }
        }
        if (!beta0) {       // beta = 1
            out = 1 / out;  // summation is for 1/p, take reciprocal to get p
        }
        std::vector<CoolPropDbl>& K = HEOS.get_K();
        for (int i = 0; i < static_cast<int>(z.size()); ++i) {
            double pci = HEOS.get_fluid_constant(i, iP_critical);
            double Tci = HEOS.get_fluid_constant(i, iT_critical);
            double omegai = HEOS.get_fluid_constant(i, iacentric_factor);
            K[i] = pci / out * exp(5.373 * (1 + omegai) * (1 - Tci / input_value));
        }
    } else {
        // Find first guess for output variable using Wilson K-factors
        WilsonK_resid Resid(HEOS, beta, input_value, input_type, z, HEOS.get_K());
        if (guess < 0 || !ValidNumber(guess))
            out = Brent(Resid, 50, 10000, 1e-10, 1e-10, 100);
        else
            out = Secant(Resid, guess, 0.001, 1e-10, 100);
        if (!ValidNumber(out)) {
            throw ValueError("saturation_p_Wilson failed to get good output value");
        }
    }
    return out;
}
struct SuccessiveSubstitutionStep
{
    CoolPropDbl T, p;
};

struct newton_raphson_twophase_options
{
    enum imposed_variable_options
    {
        NO_VARIABLE_IMPOSED = 0,
        P_IMPOSED,
        T_IMPOSED
    };
    int Nstep_max;
    std::size_t Nsteps;
    CoolPropDbl beta, omega, rhomolar_liq, rhomolar_vap, pL, pV, p, T, hmolar_liq, hmolar_vap, smolar_liq, smolar_vap;
    imposed_variable_options imposed_variable;
    std::vector<CoolPropDbl> x, y, z;
    newton_raphson_twophase_options()
      : Nstep_max(30),
        Nsteps(0),
        beta(-1),
        omega(1),
        rhomolar_liq(_HUGE),
        rhomolar_vap(_HUGE),
        pL(_HUGE),
        pV(_HUGE),
        p(_HUGE),
        T(_HUGE),
        hmolar_liq(_HUGE),
        hmolar_vap(_HUGE),
        smolar_liq(_HUGE),
        smolar_vap(_HUGE),
        imposed_variable(NO_VARIABLE_IMPOSED) {}  // Defaults
};

class newton_raphson_twophase
{
   public:
    HelmholtzEOSMixtureBackend* HEOS;
    newton_raphson_twophase_options::imposed_variable_options imposed_variable;
    CoolPropDbl error_rms, rhomolar_liq, rhomolar_vap, T, p, min_rel_change, beta;
    std::size_t N;
    bool logging;
    int Nsteps;
    Eigen::MatrixXd J;
    Eigen::Vector2d r, err_rel;
    std::vector<CoolPropDbl> K, x, y, z;
    std::vector<SuccessiveSubstitutionStep> step_logger;

    newton_raphson_twophase()
      : HEOS(NULL),
        imposed_variable(newton_raphson_twophase_options::NO_VARIABLE_IMPOSED),
        error_rms(_HUGE),
        rhomolar_liq(_HUGE),
        rhomolar_vap(_HUGE),
        T(_HUGE),
        p(_HUGE),
        min_rel_change(_HUGE),
        beta(_HUGE),
        N(0),
        logging(false),
        Nsteps(0){};

    void resize(unsigned int N);

    // Reset the state of all the internal variables
    void pre_call() {
        K.clear();
        x.clear();
        y.clear();
        step_logger.clear();
        error_rms = 1e99;
        Nsteps = 0;
        rhomolar_liq = _HUGE;
        rhomolar_vap = _HUGE;
        T = _HUGE;
        p = _HUGE;
    };

    void call(HelmholtzEOSMixtureBackend& HEOS, newton_raphson_twophase_options& IO);

    /* \brief Build the arrays for the Newton-Raphson solve
         *
         */
    void build_arrays();
};

struct newton_raphson_saturation_options
{
    enum imposed_variable_options
    {
        NO_VARIABLE_IMPOSED = 0,
        P_IMPOSED,
        RHOV_IMPOSED,
        T_IMPOSED
    };
    int Nstep_max;
    bool bubble_point;
    std::size_t Nsteps;
    CoolPropDbl omega, rhomolar_liq, rhomolar_vap, pL, pV, p, T, hmolar_liq, hmolar_vap, smolar_liq, smolar_vap;
    imposed_variable_options imposed_variable;
    std::vector<CoolPropDbl> x, y;
    newton_raphson_saturation_options()
      : bubble_point(false),
        omega(_HUGE),
        rhomolar_liq(_HUGE),
        rhomolar_vap(_HUGE),
        pL(_HUGE),
        pV(_HUGE),
        p(_HUGE),
        T(_HUGE),
        hmolar_liq(_HUGE),
        hmolar_vap(_HUGE),
        smolar_liq(_HUGE),
        smolar_vap(_HUGE),
        imposed_variable(NO_VARIABLE_IMPOSED) {
        Nstep_max = 30;
        Nsteps = 0;
    }  // Defaults
};

class newton_raphson_saturation
{
   public:
    newton_raphson_saturation_options::imposed_variable_options imposed_variable;
    CoolPropDbl error_rms, rhomolar_liq, rhomolar_vap, T, p, min_rel_change;
    std::size_t N;
    bool logging;
    bool bubble_point;
    int Nsteps;
    Eigen::MatrixXd J;
    HelmholtzEOSMixtureBackend* HEOS;
    CoolPropDbl dTsat_dPsat, dPsat_dTsat;
    std::vector<CoolPropDbl> K, x, y;
    Eigen::VectorXd r, err_rel;
    std::vector<SuccessiveSubstitutionStep> step_logger;

    newton_raphson_saturation(){};

    void resize(std::size_t N);

    // Reset the state of all the internal variables
    void pre_call() {
        step_logger.clear();
        error_rms = 1e99;
        Nsteps = 0;
        rhomolar_liq = _HUGE;
        rhomolar_vap = _HUGE;
        T = _HUGE;
        p = _HUGE;
    };

    void call(HelmholtzEOSMixtureBackend& HEOS, const std::vector<CoolPropDbl>& z, std::vector<CoolPropDbl>& z_incipient,
              newton_raphson_saturation_options& IO);

    void build_arrays();

    void check_Jacobian();
};

struct PTflash_twophase_options
{
    int Nstep_max;
    std::size_t Nsteps;
    CoolPropDbl omega, rhomolar_liq, rhomolar_vap, pL, pV, p, T;
    std::vector<CoolPropDbl> x,  
      y,                         
      z;                         
    PTflash_twophase_options() : omega(_HUGE), rhomolar_liq(_HUGE), rhomolar_vap(_HUGE), pL(_HUGE), pV(_HUGE), p(_HUGE), T(_HUGE) {
        Nstep_max = 30;
        Nsteps = 0;  // Defaults
    }
};

class PTflash_twophase
{
   public:
    double error_rms;
    bool logging;
    int Nsteps;
    Eigen::MatrixXd J;
    Eigen::VectorXd r, err_rel;
    HelmholtzEOSMixtureBackend& HEOS;
    PTflash_twophase_options& IO;
    std::vector<SuccessiveSubstitutionStep> step_logger;

    PTflash_twophase(HelmholtzEOSMixtureBackend& HEOS, PTflash_twophase_options& IO) : HEOS(HEOS), IO(IO){};

    void solve();

    void build_arrays();
};
};  // namespace SaturationSolvers

namespace StabilityRoutines {

class StabilityEvaluationClass
{
   protected:
    HelmholtzEOSMixtureBackend& HEOS;
    std::vector<double> lnK, K, K0, x, y, xL, xH;
    const std::vector<double>& z;
    double rhomolar_liq, rhomolar_vap, beta, tpd_liq, tpd_vap, DELTAG_nRT;
    double m_T,  
      m_p;       
   private:
    bool _stable;
    bool debug;

   public:
    StabilityEvaluationClass(HelmholtzEOSMixtureBackend& HEOS)
      : HEOS(HEOS),
        z(HEOS.get_mole_fractions_doubleref()),
        rhomolar_liq(-1),
        rhomolar_vap(-1),
        beta(-1),
        tpd_liq(10000),
        tpd_vap(100000),
        DELTAG_nRT(10000),
        m_T(-1),
        m_p(-1),
        _stable(false),
        debug(false){};
    void set_TP(double T, double p) {
        m_T = T;
        m_p = p;
    };
    void rho_TP_w_guesses();
    void rho_TP_global();
    void rho_TP_SRK_translated();

    void trial_compositions();
    void successive_substitution(int num_steps);
    void check_stability();
    bool is_stable() {
        trial_compositions();
        successive_substitution(3);
        check_stability();
        return _stable;
    }
    void get_liq(std::vector<double>& x, double& rhomolar) {
        x = this->x;
        rhomolar = rhomolar_liq;
    }
    void get_vap(std::vector<double>& y, double& rhomolar) {
        y = this->y;
        rhomolar = rhomolar_vap;
    }
};

}; /* namespace StabilityRoutines*/

}; /* namespace CoolProp*/

#endif