<commit_before>/**
 *  @file ThermoPhase.cpp
 * Definition file for class ThermoPhase, the base class for phases with
 * thermodynamic properties
 * (see class \link Cantera::ThermoPhase ThermoPhase\endlink).
 */

//  Copyright 2002 California Institute of Technology

#include "cantera/thermo/ThermoPhase.h"
#include "cantera/base/stringUtils.h"
#include "cantera/thermo/ThermoFactory.h"
#include "cantera/base/ctml.h"
#include "cantera/base/vec_functions.h"

#include <iomanip>
#include <fstream>

using namespace std;
using namespace ctml;

namespace Cantera
{

ThermoPhase::ThermoPhase() :
    m_spthermo(0), m_speciesData(0),
    m_phi(0.0),
    m_hasElementPotentials(false),
    m_chargeNeutralityNecessary(false),
    m_ssConvention(cSS_CONVENTION_TEMPERATURE),
    m_tlast(0.0)
{
}

ThermoPhase::~ThermoPhase()
{
    for (size_t k = 0; k < m_speciesData.size(); k++) {
        delete m_speciesData[k];
    }
    delete m_spthermo;
}

ThermoPhase::ThermoPhase(const ThermoPhase& right)  :
    m_spthermo(0),
    m_speciesData(0),
    m_phi(0.0),
    m_hasElementPotentials(false),
    m_chargeNeutralityNecessary(false),
    m_ssConvention(cSS_CONVENTION_TEMPERATURE)
{
    /*
     * Call the assignment operator
     */
    *this = operator=(right);
}

ThermoPhase& ThermoPhase::operator=(const ThermoPhase& right)
{
    /*
     * Check for self assignment.
     */
    if (this == &right) {
        return *this;
    }

    /*
     * We need to destruct first
     */
    for (size_t k = 0; k < m_speciesData.size(); k++) {
        delete m_speciesData[k];
    }
    delete m_spthermo;

    /*
     * Call the base class assignment operator
     */
    (void)Phase::operator=(right);

    /*
     * Pointer to the species thermodynamic property manager
     * We own this, so we need to do a deep copy
     */
    m_spthermo = (right.m_spthermo)->duplMyselfAsSpeciesThermo();

    /*
     * Do a deep copy of species Data, because we own this
     */
    m_speciesData.resize(m_kk);
    for (size_t k = 0; k < m_kk; k++) {
        m_speciesData[k] = new XML_Node(*(right.m_speciesData[k]));
    }

    m_phi = right.m_phi;
    m_lambdaRRT = right.m_lambdaRRT;
    m_hasElementPotentials = right.m_hasElementPotentials;
    m_chargeNeutralityNecessary = right.m_chargeNeutralityNecessary;
    m_ssConvention = right.m_ssConvention;
    m_tlast = right.m_tlast;
    return *this;
}

ThermoPhase* ThermoPhase::duplMyselfAsThermoPhase() const
{
    return new ThermoPhase(*this);
}

int ThermoPhase::activityConvention() const
{
    return cAC_CONVENTION_MOLAR;
}

int ThermoPhase::standardStateConvention() const
{
    return m_ssConvention;
}

doublereal ThermoPhase::logStandardConc(size_t k) const
{
    return log(standardConcentration(k));
}

void ThermoPhase::getActivities(doublereal* a) const
{
    getActivityConcentrations(a);
    for (size_t k = 0; k < nSpecies(); k++) {
        a[k] /= standardConcentration(k);
    }
}

void ThermoPhase::getLnActivityCoefficients(doublereal* lnac) const
{
    getActivityCoefficients(lnac);
    for (size_t k = 0; k < m_kk; k++) {
        lnac[k] = std::log(lnac[k]);
    }
}

void ThermoPhase::setState_TPX(doublereal t, doublereal p, const doublereal* x)
{
    setMoleFractions(x);
    setState_TP(t,p);
}

void ThermoPhase::setState_TPX(doublereal t, doublereal p, const compositionMap& x)
{
    setMoleFractionsByName(x);
    setState_TP(t,p);
}

void ThermoPhase::setState_TPX(doublereal t, doublereal p, const std::string& x)
{
    setMoleFractionsByName(x);
    setState_TP(t,p);
}

void ThermoPhase::setState_TPY(doublereal t, doublereal p, const doublereal* y)
{
    setMassFractions(y);
    setState_TP(t,p);
}

void ThermoPhase::setState_TPY(doublereal t, doublereal p, const compositionMap& y)
{
    setMassFractionsByName(y);
    setState_TP(t,p);
}

void ThermoPhase::setState_TPY(doublereal t, doublereal p, const std::string& y)
{
    setMassFractionsByName(y);
    setState_TP(t,p);
}

void ThermoPhase::setState_TP(doublereal t, doublereal p)
{
    setTemperature(t);
    setPressure(p);
}

void ThermoPhase::setState_PX(doublereal p, doublereal* x)
{
    setMoleFractions(x);
    setPressure(p);
}

void ThermoPhase::setState_PY(doublereal p, doublereal* y)
{
    setMassFractions(y);
    setPressure(p);
}

void ThermoPhase::setState_HP(doublereal Htarget, doublereal p,
                              doublereal dTtol)
{
    setState_HPorUV(Htarget, p, dTtol, false);
}

void ThermoPhase::setState_UV(doublereal u, doublereal v, doublereal dTtol)
{
    setState_HPorUV(u, v, dTtol, true);
}

void ThermoPhase::setState_conditional_TP(doublereal t, doublereal p, bool set_p)
{
    setTemperature(t);
    if (set_p) {
        setPressure(p);
    }
}

void ThermoPhase::setState_HPorUV(doublereal Htarget, doublereal p,
                                  doublereal dTtol, bool doUV)
{
    doublereal dt;
    doublereal v = 0.0;

    // Assign the specific volume or pressure and make sure it's positive
    if (doUV) {
        doublereal v = p;
        if (v < 1.0E-300) {
            throw CanteraError("setState_HPorUV (UV)",
                               "Input specific volume is too small or negative. v = " + fp2str(v));
        }
        setDensity(1.0/v);
    } else {
        if (p < 1.0E-300) {
            throw CanteraError("setState_HPorUV (HP)",
                               "Input pressure is too small or negative. p = " + fp2str(p));
        }
        setPressure(p);
    }
    double Tmax = maxTemp() + 0.1;
    double Tmin = minTemp() - 0.1;

    // Make sure we are within the temperature bounds at the start
    // of the iteration
    double Tnew = temperature();
    double Tinit = Tnew;
    if (Tnew > Tmax) {
        Tnew = Tmax - 1.0;
    } else if (Tnew < Tmin) {
        Tnew = Tmin + 1.0;
    }
    if (Tnew != Tinit) {
        setState_conditional_TP(Tnew, p, !doUV);
    }

    double Hnew = (doUV) ? intEnergy_mass() : enthalpy_mass();
    double Cpnew = (doUV) ? cv_mass() : cp_mass();
    double Htop = Hnew;
    double Ttop = Tnew;
    double Hbot = Hnew;
    double Tbot = Tnew;
    double Told = Tnew;
    double Hold = Hnew;

    bool ignoreBounds = false;
    // Unstable phases are those for which
    // cp < 0.0. These are possible for cases where
    // we have passed the spinodal curve.
    bool unstablePhase = false;
    // Counter indicating the last temperature point where the
    // phase was unstable
    double Tunstable = -1.0;
    bool unstablePhaseNew = false;

    // Newton iteration
    for (int n = 0; n < 500; n++) {
        Told = Tnew;
        Hold = Hnew;
        double cpd = Cpnew;
        if (cpd < 0.0) {
            unstablePhase = true;
            Tunstable = Tnew;
        }
        // limit step size to 100 K
        dt = clip((Htarget - Hold)/cpd, -100.0, 100.0);

        // Calculate the new T
        Tnew = Told + dt;

        // Limit the step size so that we are convergent
        // This is the step that makes it different from a
        // Newton's algorithm
        if ((dt > 0.0 && unstablePhase) || (dt <= 0.0 && !unstablePhase)) {
            if (Hbot < Htarget && Tnew < (0.75 * Tbot + 0.25 * Told)) {
                dt = 0.75 * (Tbot - Told);
                Tnew = Told + dt;
            }
        } else if (Htop > Htarget && Tnew > (0.75 * Ttop + 0.25 * Told)) {
            dt = 0.75 * (Ttop - Told);
            Tnew = Told + dt;
        }

        // Check Max and Min values
        if (Tnew > Tmax && !ignoreBounds) {
            setState_conditional_TP(Tmax, p, !doUV);
            double Hmax = (doUV) ? intEnergy_mass() : enthalpy_mass();
            if (Hmax >= Htarget) {
                if (Htop < Htarget) {
                    Ttop = Tmax;
                    Htop = Hmax;
                }
            } else {
                Tnew = Tmax + 1.0;
                ignoreBounds = true;
            }
        }
        if (Tnew < Tmin && !ignoreBounds) {
            setState_conditional_TP(Tmin, p, !doUV);
            double Hmin = (doUV) ? intEnergy_mass() : enthalpy_mass();
            if (Hmin <= Htarget) {
                if (Hbot > Htarget) {
                    Tbot = Tmin;
                    Hbot = Hmin;
                }
            } else {
                Tnew = Tmin - 1.0;
                ignoreBounds = true;
            }
        }

        // Try to keep phase within its region of stability
        // -> Could do a lot better if I calculate the
        //    spinodal value of H.
        for (int its = 0; its < 10; its++) {
            Tnew = Told + dt;
            if (Tnew < Told / 3.0) {
                Tnew = Told / 3.0;
                dt = -2.0 * Told / 3.0;
            }
            setState_conditional_TP(Tnew, p, !doUV);
            if (doUV) {
                Hnew = intEnergy_mass();
                Cpnew = cv_mass();
            } else {
                Hnew = enthalpy_mass();
                Cpnew = cp_mass();
            }
            if (Cpnew < 0.0) {
                unstablePhaseNew = true;
                Tunstable = Tnew;
            } else {
                unstablePhaseNew = false;
                break;
            }
            if (unstablePhase == false && unstablePhaseNew == true) {
                dt *= 0.25;
            }
        }

        if (Hnew == Htarget) {
            return;
        } else if (Hnew > Htarget && (Htop < Htarget || Hnew < Htop)) {
            Htop = Hnew;
            Ttop = Tnew;
        } else if (Hnew < Htarget && (Hbot > Htarget || Hnew > Hbot)) {
            Hbot = Hnew;
            Tbot = Tnew;
        }
        // Convergence in H
        double Herr = Htarget - Hnew;
        double acpd = std::max(fabs(cpd), 1.0E-5);
        double denom = std::max(fabs(Htarget), acpd * dTtol);
        double HConvErr = fabs((Herr)/denom);
        if (HConvErr < 0.00001 *dTtol || fabs(dt) < dTtol) {
            return;
        }
    }
    // We are here when there hasn't been convergence
    /*
     * Formulate a detailed error message, since questions seem to
     * arise often about the lack of convergence.
     */
    string ErrString =  "No convergence in 500 iterations\n";
    if (doUV) {
        ErrString += "\tTarget Internal Energy  = " + fp2str(Htarget) + "\n";
        ErrString += "\tCurrent Specific Volume = " + fp2str(v) + "\n";
        ErrString += "\tStarting Temperature    = " + fp2str(Tinit) + "\n";
        ErrString += "\tCurrent Temperature     = " + fp2str(Tnew) + "\n";
        ErrString += "\tCurrent Internal Energy = " + fp2str(Hnew) + "\n";
        ErrString += "\tCurrent Delta T         = " + fp2str(dt) + "\n";
    } else {
        ErrString += "\tTarget Enthalpy         = " + fp2str(Htarget) + "\n";
        ErrString += "\tCurrent Pressure        = " + fp2str(p) + "\n";
        ErrString += "\tStarting Temperature    = " + fp2str(Tinit) + "\n";
        ErrString += "\tCurrent Temperature     = " + fp2str(Tnew) + "\n";
        ErrString += "\tCurrent Enthalpy        = " + fp2str(Hnew) + "\n";
        ErrString += "\tCurrent Delta T         = " + fp2str(dt) + "\n";
    }
    if (unstablePhase) {
        ErrString += "\t  - The phase became unstable (Cp < 0) T_unstable_last = "
                     + fp2str(Tunstable) + "\n";
    }
    if (doUV) {
        throw CanteraError("setState_HPorUV (UV)", ErrString);
    } else {
        throw CanteraError("setState_HPorUV (HP)", ErrString);
    }
}

void ThermoPhase::setState_SP(doublereal Starget, doublereal p,
                              doublereal dTtol)
{
    setState_SPorSV(Starget, p, dTtol, false);
}

void ThermoPhase::setState_SV(doublereal Starget, doublereal v,
                              doublereal dTtol)
{
    setState_SPorSV(Starget, v, dTtol, true);
}

void ThermoPhase::setState_SPorSV(doublereal Starget, doublereal p,
                                  doublereal dTtol, bool doSV)
{
    doublereal v = 0.0;
    doublereal dt;
    if (doSV) {
        v = p;
        if (v < 1.0E-300) {
            throw CanteraError("setState_SPorSV (SV)",
                               "Input specific volume is too small or negative. v = " + fp2str(v));
        }
        setDensity(1.0/v);
    } else {
        if (p < 1.0E-300) {
            throw CanteraError("setState_SPorSV (SP)",
                               "Input pressure is too small or negative. p = " + fp2str(p));
        }
        setPressure(p);
    }
    double Tmax = maxTemp() + 0.1;
    double Tmin = minTemp() - 0.1;

    // Make sure we are within the temperature bounds at the start
    // of the iteration
    double Tnew = temperature();
    double Tinit = Tnew;
    if (Tnew > Tmax) {
        Tnew = Tmax - 1.0;
    } else if (Tnew < Tmin) {
        Tnew = Tmin + 1.0;
    }
    if (Tnew != Tinit) {
        setState_conditional_TP(Tnew, p, !doSV);
    }

    double Snew = entropy_mass();
    double Cpnew = (doSV) ? cv_mass() : cp_mass();

    double Stop = Snew;
    double Ttop = Tnew;
    double Sbot = Snew;
    double Tbot = Tnew;
    double Told = Tnew;
    double Sold = Snew;

    bool ignoreBounds = false;
    // Unstable phases are those for which
    // Cp < 0.0. These are possible for cases where
    // we have passed the spinodal curve.
    bool unstablePhase = false;
    double Tunstable = -1.0;
    bool unstablePhaseNew = false;

    // Newton iteration
    for (int n = 0; n < 500; n++) {
        Told = Tnew;
        Sold = Snew;
        double cpd = Cpnew;
        if (cpd < 0.0) {
            unstablePhase = true;
            Tunstable = Tnew;
        }
        // limit step size to 100 K
        dt = clip((Starget - Sold)*Told/cpd, -100.0, 100.0);
        Tnew = Told + dt;

        // Limit the step size so that we are convergent
        if ((dt > 0.0 && unstablePhase) || (dt <= 0.0 && !unstablePhase)) {
            if (Sbot < Starget && Tnew < Tbot) {
                dt = 0.75 * (Tbot - Told);
                Tnew = Told + dt;
            }
        } else if (Stop > Starget && Tnew > Ttop) {
            dt = 0.75 * (Ttop - Told);
            Tnew = Told + dt;
        }

        // Check Max and Min values
        if (Tnew > Tmax && !ignoreBounds) {
            setState_conditional_TP(Tmax, p, !doSV);
            double Smax = entropy_mass();
            if (Smax >= Starget) {
                if (Stop < Starget) {
                    Ttop = Tmax;
                    Stop = Smax;
                }
            } else {
                Tnew = Tmax + 1.0;
                ignoreBounds = true;
            }
        } else if (Tnew < Tmin && !ignoreBounds) {
            setState_conditional_TP(Tmin, p, !doSV);
            double Smin = entropy_mass();
            if (Smin <= Starget) {
                if (Sbot > Starget) {
                    Tbot = Tmin;
                    Sbot = Smin;
                }
            } else {
                Tnew = Tmin - 1.0;
                ignoreBounds = true;
            }
        }

        // Try to keep phase within its region of stability
        // -> Could do a lot better if I calculate the
        //    spinodal value of H.
        for (int its = 0; its < 10; its++) {
            Tnew = Told + dt;
            setState_conditional_TP(Tnew, p, !doSV);
            Cpnew = (doSV) ? cv_mass() : cp_mass();
            Snew = entropy_mass();
            if (Cpnew < 0.0) {
                unstablePhaseNew = true;
                Tunstable = Tnew;
            } else {
                unstablePhaseNew = false;
                break;
            }
            if (unstablePhase == false && unstablePhaseNew == true) {
                dt *= 0.25;
            }
        }

        if (Snew == Starget) {
            return;
        } else if (Snew > Starget && (Stop < Starget || Snew < Stop)) {
            Stop = Snew;
            Ttop = Tnew;
        } else if (Snew < Starget && (Sbot > Starget || Snew > Sbot)) {
            Sbot = Snew;
            Tbot = Tnew;
        }
        // Convergence in S
        double Serr = Starget - Snew;
        double acpd = std::max(fabs(cpd), 1.0E-5);
        double denom = std::max(fabs(Starget), acpd * dTtol);
        double SConvErr = fabs((Serr * Tnew)/denom);
        if (SConvErr < 0.00001 *dTtol || fabs(dt) < dTtol) {
            return;
        }
    }
    // We are here when there hasn't been convergence
    /*
     * Formulate a detailed error message, since questions seem to
     * arise often about the lack of convergence.
     */
    string ErrString =  "No convergence in 500 iterations\n";
    if (doSV) {
        ErrString += "\tTarget Entropy          = " + fp2str(Starget) + "\n";
        ErrString += "\tCurrent Specific Volume = " + fp2str(v) + "\n";
        ErrString += "\tStarting Temperature    = " + fp2str(Tinit) + "\n";
        ErrString += "\tCurrent Temperature     = " + fp2str(Tnew) + "\n";
        ErrString += "\tCurrent Entropy          = " + fp2str(Snew) + "\n";
        ErrString += "\tCurrent Delta T         = " + fp2str(dt) + "\n";
    } else {
        ErrString += "\tTarget Entropy          = " + fp2str(Starget) + "\n";
        ErrString += "\tCurrent Pressure        = " + fp2str(p) + "\n";
        ErrString += "\tStarting Temperature    = " + fp2str(Tinit) + "\n";
        ErrString += "\tCurrent Temperature     = " + fp2str(Tnew) + "\n";
        ErrString += "\tCurrent Entropy         = " + fp2str(Snew) + "\n";
        ErrString += "\tCurrent Delta T         = " + fp2str(dt) + "\n";
    }
    if (unstablePhase) {
        ErrString += "\t  - The phase became unstable (Cp < 0) T_unstable_last = "
                     + fp2str(Tunstable) + "\n";
    }
    if (doSV) {
        throw CanteraError("setState_SPorSV (SV)", ErrString);
    } else {
        throw CanteraError("setState_SPorSV (SP)", ErrString);
    }
}

void ThermoPhase::getUnitsStandardConc(double* uA, int k, int sizeUA) const
{
    for (int i = 0; i < sizeUA; i++) {
        if (i == 0) {
            uA[0] = 1.0;
        }
        if (i == 1) {
            uA[1] = -int(nDim());
        }
        if (i == 2) {
            uA[2] = 0.0;
        }
        if (i == 3) {
            uA[3] = 0.0;
        }
        if (i == 4) {
            uA[4] = 0.0;
        }
        if (i == 5) {
            uA[5] = 0.0;
        }
    }
}

void ThermoPhase::setSpeciesThermo(SpeciesThermo* spthermo)
{
    if (m_spthermo) {
        if (m_spthermo != spthermo) {
            delete m_spthermo;
        }
    }
    m_spthermo = spthermo;
}

SpeciesThermo& ThermoPhase::speciesThermo(int k)
{
    if (!m_spthermo) {
        throw CanteraError("ThermoPhase::speciesThermo()",
                           "species reference state thermo manager was not set");
    }
    return *m_spthermo;
}

void ThermoPhase::initThermoFile(const std::string& inputFile,
                                 const std::string& id)
{
    if (inputFile.size() == 0) {
        throw CanteraError("ThermoPhase::initThermoFile",
                           "input file is null");
    }
    string path = findInputFile(inputFile);
    ifstream fin(path.c_str());
    if (!fin) {
        throw CanteraError("initThermoFile","could not open "
                           +path+" for reading.");
    }
    /*
     * The phase object automatically constructs an XML object.
     * Use this object to store information.
     */
    XML_Node* fxml = new XML_Node();
    fxml->build(fin);
    XML_Node* fxml_phase = findXMLPhase(fxml, id);
    if (!fxml_phase) {
        throw CanteraError("ThermoPhase::initThermo",
                           "ERROR: Can not find phase named " +
                           id + " in file named " + inputFile);
    }
    bool m_ok = importPhase(*fxml_phase, this);
    if (!m_ok) {
        throw CanteraError("ThermoPhase::initThermoFile","importPhase failed ");
    }
    delete fxml;
}

void ThermoPhase::initThermoXML(XML_Node& phaseNode, const std::string& id)
{
    if (phaseNode.hasChild("state")) {
        XML_Node& stateNode = phaseNode.child("state");
        setStateFromXML(stateNode);
    }
    setReferenceComposition(0);
}<commit_msg>[Thermo] Use findXMLPhase in initThermoFile<commit_after>/**
 *  @file ThermoPhase.cpp
 * Definition file for class ThermoPhase, the base class for phases with
 * thermodynamic properties
 * (see class \link Cantera::ThermoPhase ThermoPhase\endlink).
 */

//  Copyright 2002 California Institute of Technology

#include "cantera/thermo/ThermoPhase.h"
#include "cantera/base/stringUtils.h"
#include "cantera/thermo/ThermoFactory.h"
#include "cantera/base/ctml.h"
#include "cantera/base/vec_functions.h"

#include <iomanip>
#include <fstream>

using namespace std;
using namespace ctml;

namespace Cantera
{

ThermoPhase::ThermoPhase() :
    m_spthermo(0), m_speciesData(0),
    m_phi(0.0),
    m_hasElementPotentials(false),
    m_chargeNeutralityNecessary(false),
    m_ssConvention(cSS_CONVENTION_TEMPERATURE),
    m_tlast(0.0)
{
}

ThermoPhase::~ThermoPhase()
{
    for (size_t k = 0; k < m_speciesData.size(); k++) {
        delete m_speciesData[k];
    }
    delete m_spthermo;
}

ThermoPhase::ThermoPhase(const ThermoPhase& right)  :
    m_spthermo(0),
    m_speciesData(0),
    m_phi(0.0),
    m_hasElementPotentials(false),
    m_chargeNeutralityNecessary(false),
    m_ssConvention(cSS_CONVENTION_TEMPERATURE)
{
    /*
     * Call the assignment operator
     */
    *this = operator=(right);
}

ThermoPhase& ThermoPhase::operator=(const ThermoPhase& right)
{
    /*
     * Check for self assignment.
     */
    if (this == &right) {
        return *this;
    }

    /*
     * We need to destruct first
     */
    for (size_t k = 0; k < m_speciesData.size(); k++) {
        delete m_speciesData[k];
    }
    delete m_spthermo;

    /*
     * Call the base class assignment operator
     */
    (void)Phase::operator=(right);

    /*
     * Pointer to the species thermodynamic property manager
     * We own this, so we need to do a deep copy
     */
    m_spthermo = (right.m_spthermo)->duplMyselfAsSpeciesThermo();

    /*
     * Do a deep copy of species Data, because we own this
     */
    m_speciesData.resize(m_kk);
    for (size_t k = 0; k < m_kk; k++) {
        m_speciesData[k] = new XML_Node(*(right.m_speciesData[k]));
    }

    m_phi = right.m_phi;
    m_lambdaRRT = right.m_lambdaRRT;
    m_hasElementPotentials = right.m_hasElementPotentials;
    m_chargeNeutralityNecessary = right.m_chargeNeutralityNecessary;
    m_ssConvention = right.m_ssConvention;
    m_tlast = right.m_tlast;
    return *this;
}

ThermoPhase* ThermoPhase::duplMyselfAsThermoPhase() const
{
    return new ThermoPhase(*this);
}

int ThermoPhase::activityConvention() const
{
    return cAC_CONVENTION_MOLAR;
}

int ThermoPhase::standardStateConvention() const
{
    return m_ssConvention;
}

doublereal ThermoPhase::logStandardConc(size_t k) const
{
    return log(standardConcentration(k));
}

void ThermoPhase::getActivities(doublereal* a) const
{
    getActivityConcentrations(a);
    for (size_t k = 0; k < nSpecies(); k++) {
        a[k] /= standardConcentration(k);
    }
}

void ThermoPhase::getLnActivityCoefficients(doublereal* lnac) const
{
    getActivityCoefficients(lnac);
    for (size_t k = 0; k < m_kk; k++) {
        lnac[k] = std::log(lnac[k]);
    }
}

void ThermoPhase::setState_TPX(doublereal t, doublereal p, const doublereal* x)
{
    setMoleFractions(x);
    setState_TP(t,p);
}

void ThermoPhase::setState_TPX(doublereal t, doublereal p, const compositionMap& x)
{
    setMoleFractionsByName(x);
    setState_TP(t,p);
}

void ThermoPhase::setState_TPX(doublereal t, doublereal p, const std::string& x)
{
    setMoleFractionsByName(x);
    setState_TP(t,p);
}

void ThermoPhase::setState_TPY(doublereal t, doublereal p, const doublereal* y)
{
    setMassFractions(y);
    setState_TP(t,p);
}

void ThermoPhase::setState_TPY(doublereal t, doublereal p, const compositionMap& y)
{
    setMassFractionsByName(y);
    setState_TP(t,p);
}

void ThermoPhase::setState_TPY(doublereal t, doublereal p, const std::string& y)
{
    setMassFractionsByName(y);
    setState_TP(t,p);
}

void ThermoPhase::setState_TP(doublereal t, doublereal p)
{
    setTemperature(t);
    setPressure(p);
}

void ThermoPhase::setState_PX(doublereal p, doublereal* x)
{
    setMoleFractions(x);
    setPressure(p);
}

void ThermoPhase::setState_PY(doublereal p, doublereal* y)
{
    setMassFractions(y);
    setPressure(p);
}

void ThermoPhase::setState_HP(doublereal Htarget, doublereal p,
                              doublereal dTtol)
{
    setState_HPorUV(Htarget, p, dTtol, false);
}

void ThermoPhase::setState_UV(doublereal u, doublereal v, doublereal dTtol)
{
    setState_HPorUV(u, v, dTtol, true);
}

void ThermoPhase::setState_conditional_TP(doublereal t, doublereal p, bool set_p)
{
    setTemperature(t);
    if (set_p) {
        setPressure(p);
    }
}

void ThermoPhase::setState_HPorUV(doublereal Htarget, doublereal p,
                                  doublereal dTtol, bool doUV)
{
    doublereal dt;
    doublereal v = 0.0;

    // Assign the specific volume or pressure and make sure it's positive
    if (doUV) {
        doublereal v = p;
        if (v < 1.0E-300) {
            throw CanteraError("setState_HPorUV (UV)",
                               "Input specific volume is too small or negative. v = " + fp2str(v));
        }
        setDensity(1.0/v);
    } else {
        if (p < 1.0E-300) {
            throw CanteraError("setState_HPorUV (HP)",
                               "Input pressure is too small or negative. p = " + fp2str(p));
        }
        setPressure(p);
    }
    double Tmax = maxTemp() + 0.1;
    double Tmin = minTemp() - 0.1;

    // Make sure we are within the temperature bounds at the start
    // of the iteration
    double Tnew = temperature();
    double Tinit = Tnew;
    if (Tnew > Tmax) {
        Tnew = Tmax - 1.0;
    } else if (Tnew < Tmin) {
        Tnew = Tmin + 1.0;
    }
    if (Tnew != Tinit) {
        setState_conditional_TP(Tnew, p, !doUV);
    }

    double Hnew = (doUV) ? intEnergy_mass() : enthalpy_mass();
    double Cpnew = (doUV) ? cv_mass() : cp_mass();
    double Htop = Hnew;
    double Ttop = Tnew;
    double Hbot = Hnew;
    double Tbot = Tnew;
    double Told = Tnew;
    double Hold = Hnew;

    bool ignoreBounds = false;
    // Unstable phases are those for which
    // cp < 0.0. These are possible for cases where
    // we have passed the spinodal curve.
    bool unstablePhase = false;
    // Counter indicating the last temperature point where the
    // phase was unstable
    double Tunstable = -1.0;
    bool unstablePhaseNew = false;

    // Newton iteration
    for (int n = 0; n < 500; n++) {
        Told = Tnew;
        Hold = Hnew;
        double cpd = Cpnew;
        if (cpd < 0.0) {
            unstablePhase = true;
            Tunstable = Tnew;
        }
        // limit step size to 100 K
        dt = clip((Htarget - Hold)/cpd, -100.0, 100.0);

        // Calculate the new T
        Tnew = Told + dt;

        // Limit the step size so that we are convergent
        // This is the step that makes it different from a
        // Newton's algorithm
        if ((dt > 0.0 && unstablePhase) || (dt <= 0.0 && !unstablePhase)) {
            if (Hbot < Htarget && Tnew < (0.75 * Tbot + 0.25 * Told)) {
                dt = 0.75 * (Tbot - Told);
                Tnew = Told + dt;
            }
        } else if (Htop > Htarget && Tnew > (0.75 * Ttop + 0.25 * Told)) {
            dt = 0.75 * (Ttop - Told);
            Tnew = Told + dt;
        }

        // Check Max and Min values
        if (Tnew > Tmax && !ignoreBounds) {
            setState_conditional_TP(Tmax, p, !doUV);
            double Hmax = (doUV) ? intEnergy_mass() : enthalpy_mass();
            if (Hmax >= Htarget) {
                if (Htop < Htarget) {
                    Ttop = Tmax;
                    Htop = Hmax;
                }
            } else {
                Tnew = Tmax + 1.0;
                ignoreBounds = true;
            }
        }
        if (Tnew < Tmin && !ignoreBounds) {
            setState_conditional_TP(Tmin, p, !doUV);
            double Hmin = (doUV) ? intEnergy_mass() : enthalpy_mass();
            if (Hmin <= Htarget) {
                if (Hbot > Htarget) {
                    Tbot = Tmin;
                    Hbot = Hmin;
                }
            } else {
                Tnew = Tmin - 1.0;
                ignoreBounds = true;
            }
        }

        // Try to keep phase within its region of stability
        // -> Could do a lot better if I calculate the
        //    spinodal value of H.
        for (int its = 0; its < 10; its++) {
            Tnew = Told + dt;
            if (Tnew < Told / 3.0) {
                Tnew = Told / 3.0;
                dt = -2.0 * Told / 3.0;
            }
            setState_conditional_TP(Tnew, p, !doUV);
            if (doUV) {
                Hnew = intEnergy_mass();
                Cpnew = cv_mass();
            } else {
                Hnew = enthalpy_mass();
                Cpnew = cp_mass();
            }
            if (Cpnew < 0.0) {
                unstablePhaseNew = true;
                Tunstable = Tnew;
            } else {
                unstablePhaseNew = false;
                break;
            }
            if (unstablePhase == false && unstablePhaseNew == true) {
                dt *= 0.25;
            }
        }

        if (Hnew == Htarget) {
            return;
        } else if (Hnew > Htarget && (Htop < Htarget || Hnew < Htop)) {
            Htop = Hnew;
            Ttop = Tnew;
        } else if (Hnew < Htarget && (Hbot > Htarget || Hnew > Hbot)) {
            Hbot = Hnew;
            Tbot = Tnew;
        }
        // Convergence in H
        double Herr = Htarget - Hnew;
        double acpd = std::max(fabs(cpd), 1.0E-5);
        double denom = std::max(fabs(Htarget), acpd * dTtol);
        double HConvErr = fabs((Herr)/denom);
        if (HConvErr < 0.00001 *dTtol || fabs(dt) < dTtol) {
            return;
        }
    }
    // We are here when there hasn't been convergence
    /*
     * Formulate a detailed error message, since questions seem to
     * arise often about the lack of convergence.
     */
    string ErrString =  "No convergence in 500 iterations\n";
    if (doUV) {
        ErrString += "\tTarget Internal Energy  = " + fp2str(Htarget) + "\n";
        ErrString += "\tCurrent Specific Volume = " + fp2str(v) + "\n";
        ErrString += "\tStarting Temperature    = " + fp2str(Tinit) + "\n";
        ErrString += "\tCurrent Temperature     = " + fp2str(Tnew) + "\n";
        ErrString += "\tCurrent Internal Energy = " + fp2str(Hnew) + "\n";
        ErrString += "\tCurrent Delta T         = " + fp2str(dt) + "\n";
    } else {
        ErrString += "\tTarget Enthalpy         = " + fp2str(Htarget) + "\n";
        ErrString += "\tCurrent Pressure        = " + fp2str(p) + "\n";
        ErrString += "\tStarting Temperature    = " + fp2str(Tinit) + "\n";
        ErrString += "\tCurrent Temperature     = " + fp2str(Tnew) + "\n";
        ErrString += "\tCurrent Enthalpy        = " + fp2str(Hnew) + "\n";
        ErrString += "\tCurrent Delta T         = " + fp2str(dt) + "\n";
    }
    if (unstablePhase) {
        ErrString += "\t  - The phase became unstable (Cp < 0) T_unstable_last = "
                     + fp2str(Tunstable) + "\n";
    }
    if (doUV) {
        throw CanteraError("setState_HPorUV (UV)", ErrString);
    } else {
        throw CanteraError("setState_HPorUV (HP)", ErrString);
    }
}

void ThermoPhase::setState_SP(doublereal Starget, doublereal p,
                              doublereal dTtol)
{
    setState_SPorSV(Starget, p, dTtol, false);
}

void ThermoPhase::setState_SV(doublereal Starget, doublereal v,
                              doublereal dTtol)
{
    setState_SPorSV(Starget, v, dTtol, true);
}

void ThermoPhase::setState_SPorSV(doublereal Starget, doublereal p,
                                  doublereal dTtol, bool doSV)
{
    doublereal v = 0.0;
    doublereal dt;
    if (doSV) {
        v = p;
        if (v < 1.0E-300) {
            throw CanteraError("setState_SPorSV (SV)",
                               "Input specific volume is too small or negative. v = " + fp2str(v));
        }
        setDensity(1.0/v);
    } else {
        if (p < 1.0E-300) {
            throw CanteraError("setState_SPorSV (SP)",
                               "Input pressure is too small or negative. p = " + fp2str(p));
        }
        setPressure(p);
    }
    double Tmax = maxTemp() + 0.1;
    double Tmin = minTemp() - 0.1;

    // Make sure we are within the temperature bounds at the start
    // of the iteration
    double Tnew = temperature();
    double Tinit = Tnew;
    if (Tnew > Tmax) {
        Tnew = Tmax - 1.0;
    } else if (Tnew < Tmin) {
        Tnew = Tmin + 1.0;
    }
    if (Tnew != Tinit) {
        setState_conditional_TP(Tnew, p, !doSV);
    }

    double Snew = entropy_mass();
    double Cpnew = (doSV) ? cv_mass() : cp_mass();

    double Stop = Snew;
    double Ttop = Tnew;
    double Sbot = Snew;
    double Tbot = Tnew;
    double Told = Tnew;
    double Sold = Snew;

    bool ignoreBounds = false;
    // Unstable phases are those for which
    // Cp < 0.0. These are possible for cases where
    // we have passed the spinodal curve.
    bool unstablePhase = false;
    double Tunstable = -1.0;
    bool unstablePhaseNew = false;

    // Newton iteration
    for (int n = 0; n < 500; n++) {
        Told = Tnew;
        Sold = Snew;
        double cpd = Cpnew;
        if (cpd < 0.0) {
            unstablePhase = true;
            Tunstable = Tnew;
        }
        // limit step size to 100 K
        dt = clip((Starget - Sold)*Told/cpd, -100.0, 100.0);
        Tnew = Told + dt;

        // Limit the step size so that we are convergent
        if ((dt > 0.0 && unstablePhase) || (dt <= 0.0 && !unstablePhase)) {
            if (Sbot < Starget && Tnew < Tbot) {
                dt = 0.75 * (Tbot - Told);
                Tnew = Told + dt;
            }
        } else if (Stop > Starget && Tnew > Ttop) {
            dt = 0.75 * (Ttop - Told);
            Tnew = Told + dt;
        }

        // Check Max and Min values
        if (Tnew > Tmax && !ignoreBounds) {
            setState_conditional_TP(Tmax, p, !doSV);
            double Smax = entropy_mass();
            if (Smax >= Starget) {
                if (Stop < Starget) {
                    Ttop = Tmax;
                    Stop = Smax;
                }
            } else {
                Tnew = Tmax + 1.0;
                ignoreBounds = true;
            }
        } else if (Tnew < Tmin && !ignoreBounds) {
            setState_conditional_TP(Tmin, p, !doSV);
            double Smin = entropy_mass();
            if (Smin <= Starget) {
                if (Sbot > Starget) {
                    Tbot = Tmin;
                    Sbot = Smin;
                }
            } else {
                Tnew = Tmin - 1.0;
                ignoreBounds = true;
            }
        }

        // Try to keep phase within its region of stability
        // -> Could do a lot better if I calculate the
        //    spinodal value of H.
        for (int its = 0; its < 10; its++) {
            Tnew = Told + dt;
            setState_conditional_TP(Tnew, p, !doSV);
            Cpnew = (doSV) ? cv_mass() : cp_mass();
            Snew = entropy_mass();
            if (Cpnew < 0.0) {
                unstablePhaseNew = true;
                Tunstable = Tnew;
            } else {
                unstablePhaseNew = false;
                break;
            }
            if (unstablePhase == false && unstablePhaseNew == true) {
                dt *= 0.25;
            }
        }

        if (Snew == Starget) {
            return;
        } else if (Snew > Starget && (Stop < Starget || Snew < Stop)) {
            Stop = Snew;
            Ttop = Tnew;
        } else if (Snew < Starget && (Sbot > Starget || Snew > Sbot)) {
            Sbot = Snew;
            Tbot = Tnew;
        }
        // Convergence in S
        double Serr = Starget - Snew;
        double acpd = std::max(fabs(cpd), 1.0E-5);
        double denom = std::max(fabs(Starget), acpd * dTtol);
        double SConvErr = fabs((Serr * Tnew)/denom);
        if (SConvErr < 0.00001 *dTtol || fabs(dt) < dTtol) {
            return;
        }
    }
    // We are here when there hasn't been convergence
    /*
     * Formulate a detailed error message, since questions seem to
     * arise often about the lack of convergence.
     */
    string ErrString =  "No convergence in 500 iterations\n";
    if (doSV) {
        ErrString += "\tTarget Entropy          = " + fp2str(Starget) + "\n";
        ErrString += "\tCurrent Specific Volume = " + fp2str(v) + "\n";
        ErrString += "\tStarting Temperature    = " + fp2str(Tinit) + "\n";
        ErrString += "\tCurrent Temperature     = " + fp2str(Tnew) + "\n";
        ErrString += "\tCurrent Entropy          = " + fp2str(Snew) + "\n";
        ErrString += "\tCurrent Delta T         = " + fp2str(dt) + "\n";
    } else {
        ErrString += "\tTarget Entropy          = " + fp2str(Starget) + "\n";
        ErrString += "\tCurrent Pressure        = " + fp2str(p) + "\n";
        ErrString += "\tStarting Temperature    = " + fp2str(Tinit) + "\n";
        ErrString += "\tCurrent Temperature     = " + fp2str(Tnew) + "\n";
        ErrString += "\tCurrent Entropy         = " + fp2str(Snew) + "\n";
        ErrString += "\tCurrent Delta T         = " + fp2str(dt) + "\n";
    }
    if (unstablePhase) {
        ErrString += "\t  - The phase became unstable (Cp < 0) T_unstable_last = "
                     + fp2str(Tunstable) + "\n";
    }
    if (doSV) {
        throw CanteraError("setState_SPorSV (SV)", ErrString);
    } else {
        throw CanteraError("setState_SPorSV (SP)", ErrString);
    }
}

void ThermoPhase::getUnitsStandardConc(double* uA, int k, int sizeUA) const
{
    for (int i = 0; i < sizeUA; i++) {
        if (i == 0) {
            uA[0] = 1.0;
        }
        if (i == 1) {
            uA[1] = -int(nDim());
        }
        if (i == 2) {
            uA[2] = 0.0;
        }
        if (i == 3) {
            uA[3] = 0.0;
        }
        if (i == 4) {
            uA[4] = 0.0;
        }
        if (i == 5) {
            uA[5] = 0.0;
        }
    }
}

void ThermoPhase::setSpeciesThermo(SpeciesThermo* spthermo)
{
    if (m_spthermo) {
        if (m_spthermo != spthermo) {
            delete m_spthermo;
        }
    }
    m_spthermo = spthermo;
}

SpeciesThermo& ThermoPhase::speciesThermo(int k)
{
    if (!m_spthermo) {
        throw CanteraError("ThermoPhase::speciesThermo()",
                           "species reference state thermo manager was not set");
    }
    return *m_spthermo;
}

void ThermoPhase::initThermoFile(const std::string& inputFile,
                                 const std::string& id)
{
    XML_Node* fxml = get_XML_File(inputFile);
    XML_Node* fxml_phase = findXMLPhase(fxml, id);
    if (!fxml_phase) {
        throw CanteraError("ThermoPhase::initThermoFile",
                           "ERROR: Can not find phase named " +
                           id + " in file named " + inputFile);
    }
    importPhase(*fxml_phase, this);
}

void ThermoPhase::initThermoXML(XML_Node& phaseNode, const std::string& id)
{
    if (phaseNode.hasChild("state")) {
        XML_Node& stateNode = phaseNode.child("state");
        setStateFromXML(stateNode);
    }
    setReferenceComposition(0);
}
{'fraction_non_alphanumeric': 0.12280301879147039, 'fraction_numerical': 0.01263739046750747, 'mean_word_length': 1.8070661832657995, 'pattern_counts': {'":': 0, '<': 101, '<?xml version=': 0, '>': 54, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 12}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '5294551', 'n_tokens_mistral': 13773, 'n_tokens_neox': 12790, 'n_words': 3258}