CubicBackend.h

/*

The goal of this backend is to allow the Helmholtz-based structure for cubics and to replace the entire
multi-fluid model with a one-fluid model.  The other changes are relatively trivial.  The primary
change is to replace the core residual Helmholtz energy derivatives from HelmholtzEOSMixtureBackend
with the derivatives from this class.

The core code for the Helmholtz translations is from the publication
"Helmholtz energy translations for common cubic equations of state for use in one-fluid and multi-fluid mixture models"
by Ian H. Bell and Andreas Jaeger, J. Res. NIST, 2016

*/

#ifndef CUBICBACKEND_H_
#define CUBICBACKEND_H_

#include "CoolPropTools.h"
#include "DataStructures.h"
#include "GeneralizedCubic.h"
#include "CubicsLibrary.h"
#include "Configuration.h"
#include "AbstractState.h"
#include "Backends/Helmholtz/HelmholtzEOSMixtureBackend.h"
#include "Exceptions.h"
#include <vector>

namespace CoolProp {

// Forward declaration for use in initialization of AbstractCubicBackend
class CubicResidualHelmholtz;

class AbstractCubicBackend : public HelmholtzEOSMixtureBackend
{
   protected:
    shared_ptr<AbstractCubic> cubic;
    std::vector<CubicLibrary::CubicsValues> components;  
   public:
    void setup(bool generate_SatL_and_SatV = true);

    void set_alpha_from_components();

    void set_alpha0_from_components();

    shared_ptr<AbstractCubic>& get_cubic() {
        return cubic;
    };

    std::vector<std::string> calc_fluid_names(void);

    bool using_mole_fractions(void) {
        return true;
    };
    bool using_mass_fractions(void) {
        return false;
    };
    bool using_volu_fractions(void) {
        return false;
    };

    void set_mass_fractions(const std::vector<CoolPropDbl>& mass_fractions) {
        throw NotImplementedError("Mass composition has not been implemented.");
    };
    void set_volu_fractions(const std::vector<CoolPropDbl>& volu_fractions) {
        throw NotImplementedError("Volume composition has not been implemented.");
    };
    const std::vector<CoolPropDbl>& get_mole_fractions(void) {
        return this->mole_fractions;
    };

    const double get_fluid_constant(std::size_t i, parameters param) const {
        switch (param) {
            case iP_critical:
                return cubic->get_pc()[i];
            case iT_reducing:
            case iT_critical:
                return cubic->get_Tc()[i];
            case iacentric_factor:
                return cubic->get_acentric()[i];
            case imolar_mass:
                return components[i].molemass;
            case iT_triple:
                return HelmholtzEOSMixtureBackend::get_components()[i].EOS().sat_min_liquid.T;  // From the base class data structure
            case iP_triple:
                return HelmholtzEOSMixtureBackend::get_components()[i].EOS().sat_min_liquid.p;  // From the base class data structure
            case irhomolar_reducing:
            case irhomolar_critical:
                return components[i].rhomolarc;
            case igas_constant:
                return get_config_double(R_U_CODATA);
            default:
                throw ValueError(format("I don't know what to do with this fluid constant: %s", get_parameter_information(param, "short").c_str()));
        }
    }
    
    std::string fluid_param_string(const std::string&);
    CoolPropDbl calc_gas_constant(void) {
        return cubic->get_R_u();
    };
    SimpleState calc_reducing_state_nocache(const std::vector<CoolPropDbl>& mole_fractions) {
        SimpleState reducing;
        reducing.T = cubic->get_Tr();
        reducing.rhomolar = cubic->get_rhor();
        return reducing;
    };
    CoolPropDbl calc_reduced_density(void) {
        return _rhomolar / get_cubic()->get_rhor();
    };
    CoolPropDbl calc_reciprocal_reduced_temperature(void) {
        return get_cubic()->get_Tr() / _T;
    };
    std::vector<double> spinodal_densities();

    CoolPropDbl calc_T_critical(void) {
        if (is_pure_or_pseudopure) {
            return cubic->get_Tc()[0];
        } else {
            return HelmholtzEOSMixtureBackend::calc_T_critical();
        }
    };
    CoolPropDbl calc_p_critical(void) {
        if (is_pure_or_pseudopure) {
            return cubic->get_pc()[0];
        } else {
            return HelmholtzEOSMixtureBackend::calc_p_critical();
        }
    };
    CoolPropDbl calc_acentric_factor(void) {
        if (is_pure_or_pseudopure) {
            return cubic->get_acentric()[0];
        } else {
            throw ValueError("acentric factor cannot be calculated for mixtures");
        }
    }
    CoolPropDbl calc_rhomolar_critical(void) {
        if (is_pure_or_pseudopure) {
            // Curve fit from all the pure fluids in CoolProp (thanks to recommendation of A. Kazakov)
            double v_c_Lmol = 2.14107171795 * (cubic->get_Tc()[0] / cubic->get_pc()[0] * 1000) + 0.00773144012514;  // [L/mol]
            return 1 / (v_c_Lmol / 1000.0);
        } else {
            return HelmholtzEOSMixtureBackend::calc_rhomolar_critical();
        }
    };

    void get_linear_reducing_parameters(double& rhomolar, double& T);

    void get_critical_point_starting_values(double& delta0, double& tau0);

    void get_critical_point_search_radii(double& R_delta, double& R_tau);

    bool get_critical_is_terminated(double& delta, double& tau);

    CoolPropDbl calc_alphar_deriv_nocache(const int nTau, const int nDelta, const std::vector<CoolPropDbl>& mole_fractions, const CoolPropDbl& tau,
                                          const CoolPropDbl& delta);

    CoolPropDbl calc_pressure_nocache(CoolPropDbl T, CoolPropDbl rhomolar);

    virtual void update_DmolarT();

    virtual void update(CoolProp::input_pairs input_pair, double value1, double value2);

    void rho_Tp_cubic(CoolPropDbl T, CoolPropDbl p, int& Nsolns, double& rho0, double& rho1, double& rho2);

    CoolPropDbl solver_rho_Tp_SRK(CoolPropDbl T, CoolPropDbl p, phases phase) {
        return solver_rho_Tp(T, p);
    };
    CoolPropDbl solver_rho_Tp(CoolPropDbl T, CoolPropDbl p, CoolPropDbl rho_guess = -1);

    CoolPropDbl solver_rho_Tp_global(CoolPropDbl T, CoolPropDbl p, CoolPropDbl rhomax);

    void update_TPD_state() {
        TPD_state.reset(get_copy());
    };

    void saturation(CoolProp::input_pairs inputs);

    CoolPropDbl calc_molar_mass(void);

    void set_binary_interaction_double(const std::size_t i1, const std::size_t i2, const std::string& parameter, const double value);
    double get_binary_interaction_double(const std::size_t i1, const std::size_t i2, const std::string& parameter);

    void set_binary_interaction_double(const std::string& CAS1, const std::string& CAS2, const std::string& parameter, const double value) {
        throw ValueError("set_binary_interaction_double not defined for AbstractCubic not defined for CAS #");
    }
    double get_binary_interaction_double(const std::string& CAS1, const std::string& CAS2, const std::string& parameter) {
        throw ValueError("get_binary_interaction_double not defined for AbstractCubic not defined for CAS #");
    };

    // Return a 1-1 copy of this class
    virtual HelmholtzEOSMixtureBackend* get_copy(bool generate_SatL_and_SatV = true) = 0;

    // Copy the entire kij matrix from another instance in one shot
    void copy_k(AbstractCubicBackend* donor);

    //
    void copy_all_alpha_functions(AbstractCubicBackend* donor);

    void copy_internals(AbstractCubicBackend& donor);

    // Set the cubic alpha function's constants:
    void set_cubic_alpha_C(const size_t i, const std::string& parameter, const double c1, const double c2, const double c3);

    // Set fluid parameter (currently the volume translation parameter)
    void set_fluid_parameter_double(const size_t i, const std::string& parameter, const double value);

    // Get fluid parameter (currently the volume translation parameter)
    double get_fluid_parameter_double(const size_t i, const std::string& parameter);
};

class SRKBackend : public AbstractCubicBackend
{

   public:
    SRKBackend(const std::vector<double>& Tc, const std::vector<double>& pc, const std::vector<double>& acentric, double R_u,
               bool generate_SatL_and_SatV = true) {
        cubic.reset(new SRK(Tc, pc, acentric, R_u));
        setup(generate_SatL_and_SatV);
    };
    SRKBackend(double Tc, double pc, double acentric, double R_u, bool generate_SatL_and_SatV = true) {
        cubic.reset(new SRK(Tc, pc, acentric, R_u));
        setup(generate_SatL_and_SatV);
    }
    SRKBackend(const std::vector<std::string> fluid_identifiers, const double R_u = get_config_double(R_U_CODATA),
               bool generate_SatL_and_SatV = true) {
        std::vector<double> Tc, pc, acentric;
        N = fluid_identifiers.size();
        components.resize(N);
        for (std::size_t i = 0; i < fluid_identifiers.size(); ++i) {
            components[i] = CubicLibrary::get_cubic_values(fluid_identifiers[i]);
            Tc.push_back(components[i].Tc);
            pc.push_back(components[i].pc);
            acentric.push_back(components[i].acentric);
        }
        cubic.reset(new SRK(Tc, pc, acentric, R_u));
        setup(generate_SatL_and_SatV);
    }
    HelmholtzEOSMixtureBackend* get_copy(bool generate_SatL_and_SatV = true) {
        AbstractCubicBackend* ACB = new SRKBackend(cubic->get_Tc(), cubic->get_pc(), cubic->get_acentric(), cubic->get_R_u(), generate_SatL_and_SatV);
        ACB->copy_internals(*this);
        return static_cast<HelmholtzEOSMixtureBackend*>(ACB);
    }
    std::string backend_name(void) {
        return get_backend_string(SRK_BACKEND);
    }
};

class PengRobinsonBackend : public AbstractCubicBackend
{

   public:
    PengRobinsonBackend(){};  // Default constructor (make sure you know what you are doing)
    PengRobinsonBackend(const std::vector<double>& Tc, const std::vector<double>& pc, const std::vector<double>& acentric, double R_u,
                        bool generate_SatL_and_SatV = true) {
        cubic.reset(new PengRobinson(Tc, pc, acentric, R_u));
        setup(generate_SatL_and_SatV);
    };
    PengRobinsonBackend(double Tc, double pc, double acentric, double R_u, bool generate_SatL_and_SatV = true) {
        cubic.reset(new PengRobinson(Tc, pc, acentric, R_u));
        setup(generate_SatL_and_SatV);
    };
    PengRobinsonBackend(const std::vector<std::string> fluid_identifiers, const double R_u = get_config_double(R_U_CODATA),
                        bool generate_SatL_and_SatV = true) {
        std::vector<double> Tc, pc, acentric;
        N = fluid_identifiers.size();
        components.resize(N);
        for (std::size_t i = 0; i < fluid_identifiers.size(); ++i) {
            components[i] = CubicLibrary::get_cubic_values(fluid_identifiers[i]);
            Tc.push_back(components[i].Tc);
            pc.push_back(components[i].pc);
            acentric.push_back(components[i].acentric);
        }
        cubic.reset(new PengRobinson(Tc, pc, acentric, R_u));
        setup(generate_SatL_and_SatV);
    };
    HelmholtzEOSMixtureBackend* get_copy(bool generate_SatL_and_SatV = true) {
        AbstractCubicBackend* ACB =
          new PengRobinsonBackend(cubic->get_Tc(), cubic->get_pc(), cubic->get_acentric(), cubic->get_R_u(), generate_SatL_and_SatV);
        ACB->copy_internals(*this);
        return static_cast<HelmholtzEOSMixtureBackend*>(ACB);
    }
    std::string backend_name(void) {
        return get_backend_string(PR_BACKEND);
    }
};

class CubicResidualHelmholtz : public ResidualHelmholtz
{

   protected:
    AbstractCubicBackend* ACB;

   public:
    CubicResidualHelmholtz() {
        ACB = NULL;
    };
    CubicResidualHelmholtz(AbstractCubicBackend* ACB) : ACB(ACB){};

    // copy assignment
    CubicResidualHelmholtz& operator=(CubicResidualHelmholtz& other) {
        ACB = other.ACB;
        return *this;
    }

    virtual HelmholtzDerivatives all(HelmholtzEOSMixtureBackend& HEOS, const std::vector<CoolPropDbl>& mole_fractions, double tau, double delta,
                                     bool cache_values = false) {
        HelmholtzDerivatives a;
        std::vector<double> z = std::vector<double>(mole_fractions.begin(), mole_fractions.end());
        shared_ptr<AbstractCubic>& cubic = ACB->get_cubic();
        a.alphar = cubic->alphar(tau, delta, z, 0, 0);
        a.dalphar_dtau = cubic->alphar(tau, delta, z, 1, 0);
        a.dalphar_ddelta = cubic->alphar(tau, delta, z, 0, 1);
        a.d2alphar_dtau2 = cubic->alphar(tau, delta, z, 2, 0);
        a.d2alphar_ddelta_dtau = cubic->alphar(tau, delta, z, 1, 1);
        a.d2alphar_ddelta2 = cubic->alphar(tau, delta, z, 0, 2);
        a.d3alphar_dtau3 = cubic->alphar(tau, delta, z, 3, 0);
        a.d3alphar_ddelta_dtau2 = cubic->alphar(tau, delta, z, 2, 1);
        a.d3alphar_ddelta2_dtau = cubic->alphar(tau, delta, z, 1, 2);
        a.d3alphar_ddelta3 = cubic->alphar(tau, delta, z, 0, 3);
        a.d4alphar_dtau4 = cubic->alphar(tau, delta, z, 4, 0);
        a.d4alphar_ddelta_dtau3 = cubic->alphar(tau, delta, z, 3, 1);
        a.d4alphar_ddelta2_dtau2 = cubic->alphar(tau, delta, z, 2, 2);
        a.d4alphar_ddelta3_dtau = cubic->alphar(tau, delta, z, 1, 3);
        a.d4alphar_ddelta4 = cubic->alphar(tau, delta, z, 0, 4);
        return a;
    }
    virtual CoolPropDbl dalphar_dxi(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) {
        return ACB->get_cubic()->d_alphar_dxi(HEOS.tau(), HEOS.delta(), HEOS.get_mole_fractions_doubleref(), 0, 0, i, xN_flag == XN_INDEPENDENT);
    }
    virtual CoolPropDbl d2alphar_dxi_dTau(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) {
        return ACB->get_cubic()->d_alphar_dxi(HEOS.tau(), HEOS.delta(), HEOS.get_mole_fractions_doubleref(), 1, 0, i, xN_flag == XN_INDEPENDENT);
    }
    virtual CoolPropDbl d2alphar_dxi_dDelta(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) {
        return ACB->get_cubic()->d_alphar_dxi(HEOS.tau(), HEOS.delta(), HEOS.get_mole_fractions_doubleref(), 0, 1, i, xN_flag == XN_INDEPENDENT);
    }
    virtual CoolPropDbl d3alphar_dxi_dTau2(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) {
        return ACB->get_cubic()->d_alphar_dxi(HEOS.tau(), HEOS.delta(), HEOS.get_mole_fractions_doubleref(), 2, 0, i, xN_flag == XN_INDEPENDENT);
    }
    virtual CoolPropDbl d3alphar_dxi_dDelta_dTau(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) {
        return ACB->get_cubic()->d_alphar_dxi(HEOS.tau(), HEOS.delta(), HEOS.get_mole_fractions_doubleref(), 1, 1, i, xN_flag == XN_INDEPENDENT);
    }
    virtual CoolPropDbl d3alphar_dxi_dDelta2(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) {
        return ACB->get_cubic()->d_alphar_dxi(HEOS.tau(), HEOS.delta(), HEOS.get_mole_fractions_doubleref(), 0, 2, i, xN_flag == XN_INDEPENDENT);
    }

    virtual CoolPropDbl d2alphardxidxj(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag) {
        return ACB->get_cubic()->d2_alphar_dxidxj(HEOS.tau(), HEOS.delta(), HEOS.get_mole_fractions_doubleref(), 0, 0, i, j,
                                                  xN_flag == XN_INDEPENDENT);
    }
    virtual CoolPropDbl d3alphar_dxi_dxj_dTau(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag) {
        return ACB->get_cubic()->d2_alphar_dxidxj(HEOS.tau(), HEOS.delta(), HEOS.get_mole_fractions_doubleref(), 1, 0, i, j,
                                                  xN_flag == XN_INDEPENDENT);
    }
    virtual CoolPropDbl d3alphar_dxi_dxj_dDelta(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag) {
        return ACB->get_cubic()->d2_alphar_dxidxj(HEOS.tau(), HEOS.delta(), HEOS.get_mole_fractions_doubleref(), 0, 1, i, j,
                                                  xN_flag == XN_INDEPENDENT);
    }

    virtual CoolPropDbl d4alphar_dxi_dTau3(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) {
        return ACB->get_cubic()->d_alphar_dxi(HEOS.tau(), HEOS.delta(), HEOS.get_mole_fractions_doubleref(), 3, 0, i, xN_flag == XN_INDEPENDENT);
    }
    virtual CoolPropDbl d4alphar_dxi_dDelta2_dTau(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) {
        return ACB->get_cubic()->d_alphar_dxi(HEOS.tau(), HEOS.delta(), HEOS.get_mole_fractions_doubleref(), 1, 2, i, xN_flag == XN_INDEPENDENT);
    }
    virtual CoolPropDbl d4alphar_dxi_dDelta_dTau2(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) {
        return ACB->get_cubic()->d_alphar_dxi(HEOS.tau(), HEOS.delta(), HEOS.get_mole_fractions_doubleref(), 2, 1, i, xN_flag == XN_INDEPENDENT);
    }
    virtual CoolPropDbl d4alphar_dxi_dDelta3(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) {
        return ACB->get_cubic()->d_alphar_dxi(HEOS.tau(), HEOS.delta(), HEOS.get_mole_fractions_doubleref(), 0, 3, i, xN_flag == XN_INDEPENDENT);
    }
    virtual CoolPropDbl d4alphar_dxi_dxj_dTau2(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag) {
        return ACB->get_cubic()->d2_alphar_dxidxj(HEOS.tau(), HEOS.delta(), HEOS.get_mole_fractions_doubleref(), 2, 0, i, j,
                                                  xN_flag == XN_INDEPENDENT);
    }
    virtual CoolPropDbl d4alphar_dxi_dxj_dDelta_dTau(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag) {
        return ACB->get_cubic()->d2_alphar_dxidxj(HEOS.tau(), HEOS.delta(), HEOS.get_mole_fractions_doubleref(), 1, 1, i, j,
                                                  xN_flag == XN_INDEPENDENT);
    }
    virtual CoolPropDbl d4alphar_dxi_dxj_dDelta2(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag) {
        return ACB->get_cubic()->d2_alphar_dxidxj(HEOS.tau(), HEOS.delta(), HEOS.get_mole_fractions_doubleref(), 0, 2, i, j,
                                                  xN_flag == XN_INDEPENDENT);
    }
    virtual CoolPropDbl d3alphardxidxjdxk(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j, std::size_t k,
                                          x_N_dependency_flag xN_flag) {
        return ACB->get_cubic()->d3_alphar_dxidxjdxk(HEOS.tau(), HEOS.delta(), HEOS.get_mole_fractions_doubleref(), 0, 0, i, j, k,
                                                     xN_flag == XN_INDEPENDENT);
    }
    virtual CoolPropDbl d4alphar_dxi_dxj_dxk_dTau(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j, std::size_t k,
                                                  x_N_dependency_flag xN_flag) {
        return ACB->get_cubic()->d3_alphar_dxidxjdxk(HEOS.tau(), HEOS.delta(), HEOS.get_mole_fractions_doubleref(), 1, 0, i, j, k,
                                                     xN_flag == XN_INDEPENDENT);
    }
    virtual CoolPropDbl d4alphar_dxi_dxj_dxk_dDelta(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j, std::size_t k,
                                                    x_N_dependency_flag xN_flag) {
        return ACB->get_cubic()->d3_alphar_dxidxjdxk(HEOS.tau(), HEOS.delta(), HEOS.get_mole_fractions_doubleref(), 0, 1, i, j, k,
                                                     xN_flag == XN_INDEPENDENT);
    }
};

} /* namespace CoolProp */
#endif /* CUBICBACKEND_H_ */