Phase.cpp

Go to the documentation of this file.

/**
 *  @file Phase.cpp
 *   Definition file for class Phase.
 */
 
// This file is part of Cantera. See License.txt in the top-level directory or
// at https://cantera.org/license.txt for license and copyright information.
 
#include "cantera/thermo/Phase.h"
#include "cantera/base/utilities.h"
#include "cantera/base/stringUtils.h"
#include "cantera/thermo/Species.h"
#include "cantera/thermo/ThermoFactory.h"
 
using namespace std;
 
namespace Cantera
{
 
string Phase::name() const
{
    return m_name;
}
 
void Phase::setName(const string& name)
{
    m_name = name;
}
 
size_t Phase::nElements() const
{
    return m_mm;
}
 
size_t Phase::checkElementIndex(size_t m) const
{
    if (m < m_mm) {
        return m;
    }
    throw IndexError("Phase::checkElementIndex", "elements", m, m_mm);
}
 
string Phase::elementName(size_t m) const
{
    if (m < m_mm) {
        return m_elementNames[m];
    }
    throw IndexError("Phase::elementName", "element", m, m_mm);
}
 
size_t Phase::elementIndex(const string& elementName, bool raise) const
{
    for (size_t i = 0; i < m_mm; i++) {
        if (m_elementNames[i] == elementName) {
            return i;
        }
    }
    if (!raise) {
        return npos;
    }
    throw CanteraError("Phase::elementIndex", "Element '{}' not found", elementName);
}
 
const vector<string>& Phase::elementNames() const
{
    return m_elementNames;
}
 
double Phase::atomicWeight(size_t m) const
{
    return m_atomicWeights[m];
}
 
double Phase::entropyElement298(size_t m) const
{
    return m_entropy298[checkElementIndex(m)];
}
 
span<const double> Phase::atomicWeights() const
{
    return m_atomicWeights;
}
 
int Phase::atomicNumber(size_t m) const
{
    return m_atomicNumbers[m];
}
 
int Phase::elementType(size_t m) const
{
    return m_elem_type[m];
}
 
int Phase::changeElementType(int m, int elem_type)
{
    int old = m_elem_type[m];
    m_elem_type[m] = elem_type;
    return old;
}
 
double Phase::nAtoms(size_t k, size_t m) const
{
    checkElementIndex(m);
    checkSpeciesIndex(k);
    return m_speciesComp[m_mm * k + m];
}
 
size_t Phase::findSpeciesLower(const string& name) const
{
    string nLower = toLowerCopy(name);
    size_t loc = npos;
    auto it = m_speciesLower.find(nLower);
    if (it == m_speciesLower.end()) {
        return npos;
    } else {
        loc = it->second;
        if (loc == npos) {
            throw CanteraError("Phase::findSpeciesLower",
                               "Lowercase species name '{}' is not unique. "
                               "Set Phase::caseSensitiveSpecies to true to "
                               "enforce case sensitive species names", nLower);
        }
    }
    return loc;
}
 
size_t Phase::speciesIndex(const string& name, bool raise) const
{
    size_t loc = npos;
    auto it = m_speciesIndices.find(name);
    if (it != m_speciesIndices.end()) {
        return it->second;
    } else if (!m_caseSensitiveSpecies) {
        loc = findSpeciesLower(name);
    }
    if (loc==npos && raise) {
        throw CanteraError("Phase::speciesIndex", "Species '{}' not found", name);
    }
 
    return loc;
}
 
string Phase::speciesName(size_t k) const
{
    if (k < m_kk) {
        return m_speciesNames[k];
    }
    throw IndexError("Phase::speciesName", "species", k, m_kk);
}
 
const vector<string>& Phase::speciesNames() const
{
    return m_speciesNames;
}
 
size_t Phase::checkSpeciesIndex(size_t k) const
{
    if (k < m_kk) {
        return k;
    }
    throw IndexError("Phase::checkSpeciesIndex", "species", k, m_kk);
}
 
map<string, size_t> Phase::nativeState() const
{
    if (isPure()) {
        if (isCompressible()) {
            return { {"T", 0}, {"D", 1} };
        } else {
            return { {"T", 0}, {"P", 1} };
        }
    } else {
        if (isCompressible()) {
            return { {"T", 0}, {"D", 1}, {"Y", 2} };
        } else {
            return { {"T", 0}, {"P", 1}, {"Y", 2} };
        }
    }
}
 
string Phase::nativeMode() const
{
    map<size_t, string> states; // reverse lookup
    for (auto& [name, value] : nativeState()) {
        states[value] = name;
    }
    string out;
    for (auto& [value, name] : states) {
        out += name;
    }
    return out;
}
 
vector<string> Phase::fullStates() const
{
    if (isPure()) {
        if (isCompressible()) {
            return {"TD", "TP", "UV", "DP", "HP", "SP", "SV"};
        } else {
            return {"TP", "HP", "SP"};
        }
    } else {
        if (isCompressible()) {
            return {"TDX", "TDY", "TPX", "TPY", "UVX", "UVY", "DPX", "DPY",
                    "HPX", "HPY", "SPX", "SPY", "SVX", "SVY"};
        } else {
            return {"TPX", "TPY", "HPX", "HPY", "SPX", "SPY"};
        }
    }
}
 
vector<string> Phase::partialStates() const
{
    if (isPure()) {
        return {};
    } else {
        if (isCompressible()) {
            return {"TD", "TP", "UV", "DP", "HP", "SP", "SV"};
        } else {
            return {"TP", "HP", "SP"};
        }
    }
}
 
void Phase::savePartialState(span<double> state) const
{
    checkArraySize("Phase::savePartialState", state.size(), partialStateSize());
    auto native = nativeState();
    state[native.at("T")] = temperature();
    if (isCompressible()) {
        state[native.at("D")] = density();
    } else {
        state[native.at("P")] = pressure();
    }
}
 
void Phase::restorePartialState(span<const double> state)
{
    checkArraySize("Phase::restorePartialState", state.size(), partialStateSize());
    auto native = nativeState();
    setTemperature(state[native.at("T")]);
    if (isCompressible()) {
        setDensity(state[native.at("D")]);
    } else {
        setPressure(state[native.at("P")]);
    }
}
 
size_t Phase::stateSize() const {
    if (isPure()) {
        return partialStateSize();
    } else {
        return partialStateSize() + nSpecies();
    }
}
 
void Phase::saveState(span<double> state) const
{
    checkArraySize("Phase::saveState", state.size(), stateSize());
    savePartialState(state);
    auto native = nativeState();
    if (native.count("X")) {
        getMoleFractions(state.subspan(native["X"], nSpecies()));
    } else if (native.count("Y")) {
        getMassFractions(state.subspan(native["Y"], nSpecies()));
    }
}
 
void Phase::restoreState(span<const double> state)
{
    checkArraySize("Phase::restoreState", state.size(), stateSize());
    restorePartialState(state);
    auto native = nativeState();
 
    if (native.count("X")) {
        setMoleFractions_NoNorm(state.subspan(native["X"], nSpecies()));
    } else if (native.count("Y")) {
        setMassFractions_NoNorm(state.subspan(native["Y"], nSpecies()));
    }
    compositionChanged();
}
 
void Phase::setMoleFractions(span<const double> x)
{
    checkArraySize("Phase::setMoleFractions", x.size(), m_kk);
    // Use m_y as a temporary work vector for the non-negative mole fractions
    double norm = 0.0;
    // sum is calculated below as the unnormalized molecular weight
    double sum = 0;
    for (size_t k = 0; k < m_kk; k++) {
        double xk = std::max(x[k], 0.0); // Ignore negative mole fractions
        m_y[k] = xk;
        norm += xk;
        sum += m_molwts[k] * xk;
    }
 
    // Set m_ym to the normalized mole fractions divided by the normalized mean
    // molecular weight:
    //     m_ym_k = X_k / (sum_k X_k M_k)
    const double invSum = 1.0/sum;
    if (!std::isfinite(invSum) && m_kk != 0) {
        throw CanteraError("Phase::setMoleFractions",
            "The sum of the species mole fractions is zero, too small to normalize,"
            " or not-a-number.");
    }
    for (size_t k=0; k < m_kk; k++) {
        m_ym[k] = m_y[k]*invSum;
    }
 
    // Now set m_y to the normalized mass fractions:
    //     m_y =  X_k M_k / (sum_k X_k M_k)
    for (size_t k=0; k < m_kk; k++) {
        m_y[k] = m_ym[k] * m_molwts[k];
    }
 
    // Calculate the normalized molecular weight
    m_mmw = sum/norm;
    compositionChanged();
}
 
void Phase::setMoleFractions_NoNorm(span<const double> x)
{
    checkArraySize("Phase::setMoleFractions_NoNorm", x.size(), m_kk);
    m_mmw = dot(x.begin(), x.begin() + m_kk, m_molwts.begin());
    scale(x.begin(), x.begin() + m_kk, m_ym.begin(), 1.0/m_mmw);
    transform(m_ym.begin(), m_ym.begin() + m_kk, m_molwts.begin(),
              m_y.begin(), multiplies<double>());
    compositionChanged();
}
 
void Phase::setMoleFractionsByName(const Composition& xMap)
{
    vector<double> mf = getCompositionFromMap(xMap);
    setMoleFractions(mf);
}
 
void Phase::setMoleFractionsByName(const string& x)
{
    setMoleFractionsByName(parseCompString(x));
}
 
void Phase::setMassFractions(span<const double> y)
{
    checkArraySize("Phase::setMassFractions", y.size(), m_kk);
    for (size_t k = 0; k < m_kk; k++) {
        m_y[k] = std::max(y[k], 0.0); // Ignore negative mass fractions
    }
    double invSum = 1.0 / accumulate(m_y.begin(), m_y.end(), 0.0);
    if (!std::isfinite(invSum) && m_kk != 0) {
        throw CanteraError("Phase::setMassFractions",
            "The sum of the species mass fractions is zero, too small to normalize,"
            " or not-a-number.");
    }
    scale(m_y.begin(), m_y.end(), m_y.begin(), invSum);
 
    transform(m_y.begin(), m_y.end(), m_rmolwts.begin(),
              m_ym.begin(), multiplies<double>());
    m_mmw = 1.0 / accumulate(m_ym.begin(), m_ym.end(), 0.0);
    compositionChanged();
}
 
void Phase::setMassFractions_NoNorm(span<const double> y)
{
    checkArraySize("Phase::setMassFractions_NoNorm", y.size(), m_kk);
    double sum = 0.0;
    copy(y.begin(), y.begin() + m_kk, m_y.begin());
    transform(m_y.begin(), m_y.end(), m_rmolwts.begin(), m_ym.begin(),
              multiplies<double>());
    sum = accumulate(m_ym.begin(), m_ym.end(), 0.0);
    m_mmw = 1.0/sum;
    compositionChanged();
}
 
void Phase::setMassFractionsByName(const Composition& yMap)
{
    vector<double> mf = getCompositionFromMap(yMap);
    setMassFractions(mf);
}
 
void Phase::setMassFractionsByName(const string& y)
{
    setMassFractionsByName(parseCompString(y));
}
 
void Phase::setState_TD(double t, double rho)
{
    vector<double> state(partialStateSize());
    savePartialState(state);
    try {
        setTemperature(t);
        setDensity(rho);
    } catch (std::exception&) {
        restorePartialState(state);
        throw;
    }
}
 
double Phase::molecularWeight(size_t k) const
{
    checkSpeciesIndex(k);
    return m_molwts[k];
}
 
void Phase::getMolecularWeights(span<double> weights) const
{
    auto mw = molecularWeights();
    checkArraySize("Phase::getMolecularWeights", weights.size(), mw.size());
    copy(mw.begin(), mw.end(), weights.begin());
}
 
span<const double> Phase::molecularWeights() const
{
    return m_molwts;
}
 
span<const double> Phase::inverseMolecularWeights() const
{
    return m_rmolwts;
}
 
void Phase::getCharges(span<double> charges) const
{
    checkArraySize("Phase::getCharges", charges.size(), m_speciesCharge.size());
    copy(m_speciesCharge.begin(), m_speciesCharge.end(), charges.begin());
}
 
Composition Phase::getMoleFractionsByName(double threshold) const
{
    Composition comp;
    for (size_t k = 0; k < m_kk; k++) {
        double x = moleFraction(k);
        if (x > threshold) {
            comp[speciesName(k)] = x;
        }
    }
    return comp;
}
 
Composition Phase::getMassFractionsByName(double threshold) const
{
    Composition comp;
    for (size_t k = 0; k < m_kk; k++) {
        double x = massFraction(k);
        if (x > threshold) {
            comp[speciesName(k)] = x;
        }
    }
    return comp;
}
 
void Phase::getMoleFractions(span<double> x) const
{
    checkArraySize("Phase::getMoleFractions", x.size(), m_kk);
    scale(m_ym.begin(), m_ym.end(), x.begin(), m_mmw);
}
 
double Phase::moleFraction(size_t k) const
{
    checkSpeciesIndex(k);
    return m_ym[k] * m_mmw;
}
 
double Phase::moleFraction(const string& nameSpec) const
{
    size_t iloc = speciesIndex(nameSpec, false);
    if (iloc != npos) {
        return moleFraction(iloc);
    } else {
        return 0.0;
    }
}
 
double Phase::massFraction(size_t k) const
{
    checkSpeciesIndex(k);
    return m_y[k];
}
 
double Phase::massFraction(const string& nameSpec) const
{
    size_t iloc = speciesIndex(nameSpec, false);
    if (iloc != npos) {
        return massFractions()[iloc];
    } else {
        return 0.0;
    }
}
 
void Phase::getMassFractions(span<double> y) const
{
    checkArraySize("Phase::getMassFractions", y.size(), m_kk);
    copy(m_y.begin(), m_y.end(), y.begin());
}
 
double Phase::concentration(const size_t k) const
{
    checkSpeciesIndex(k);
    return m_y[k] * m_dens * m_rmolwts[k];
}
 
void Phase::getConcentrations(span<double> c) const
{
    checkArraySize("Phase::getConcentrations", c.size(), m_kk);
    scale(m_ym.begin(), m_ym.end(), c.begin(), m_dens);
}
 
void Phase::setConcentrations(span<const double> conc)
{
    assertCompressible("setConcentrations");
    checkArraySize("Phase::setConcentrations", conc.size(), m_kk);
 
    // Use m_y as temporary storage for non-negative concentrations
    double sum = 0.0, norm = 0.0;
    for (size_t k = 0; k != m_kk; ++k) {
        double ck = std::max(conc[k], 0.0); // Ignore negative concentrations
        m_y[k] = ck;
        sum += ck * m_molwts[k];
        norm += ck;
    }
    m_mmw = sum/norm;
    setDensity(sum);
    double rsum = 1.0/sum;
    for (size_t k = 0; k != m_kk; ++k) {
        m_ym[k] = m_y[k] * rsum;
        m_y[k] = m_ym[k] * m_molwts[k]; // m_y is now the mass fraction
    }
    compositionChanged();
}
 
void Phase::setConcentrationsNoNorm(span<const double> conc)
{
    assertCompressible("setConcentrationsNoNorm");
    checkArraySize("Phase::setConcentrationsNoNorm", conc.size(), m_kk);
 
    double sum = 0.0, norm = 0.0;
    for (size_t k = 0; k != m_kk; ++k) {
        sum += conc[k] * m_molwts[k];
        norm += conc[k];
    }
    m_mmw = sum/norm;
    setDensity(sum);
    double rsum = 1.0/sum;
    for (size_t k = 0; k != m_kk; ++k) {
        m_ym[k] = conc[k] * rsum;
        m_y[k] = m_ym[k] * m_molwts[k];
    }
    compositionChanged();
}
 
void Phase::setMolesNoTruncate(span<const double> N)
{
    checkArraySize("Phase::setMolesNoTruncate", N.size(), m_kk);
    // get total moles
    copy(N.begin(), N.begin() + m_kk, m_ym.begin());
    double totalMoles = accumulate(m_ym.begin(), m_ym.end(), 0.0);
    // get total mass
    copy(N.begin(), N.begin() + m_kk, m_y.begin());
    transform(m_y.begin(), m_y.end(), m_molwts.begin(), m_y.begin(), multiplies<double>());
    double totalMass = accumulate(m_y.begin(), m_y.end(), 0.0);
    // mean molecular weight
    m_mmw = totalMass/totalMoles;
    // mass fractions
    scale(m_y.begin(), m_y.end(), m_y.begin(), 1/totalMass);
    // moles fractions/m_mmw
    scale(m_ym.begin(), m_ym.end(), m_ym.begin(), 1/(m_mmw * totalMoles));
    // composition has changed
    compositionChanged();
}
 
double Phase::elementalMassFraction(const size_t m) const
{
    double Z_m = 0.0;
    for (size_t k = 0; k != m_kk; ++k) {
        Z_m += nAtoms(k, m) * atomicWeight(m) / molecularWeight(k)
            * massFraction(k);
    }
    return Z_m;
}
 
double Phase::elementalMoleFraction(const size_t m) const
{
    double denom = 0;
    for (size_t k = 0; k < m_kk; k++) {
        double atoms = 0;
        for (size_t j = 0; j < nElements(); j++) {
            atoms += nAtoms(k, j);
        }
        denom += atoms * moleFraction(k);
    }
    double numerator = 0.0;
    for (size_t k = 0; k != m_kk; ++k) {
        numerator += nAtoms(k, m) * moleFraction(k);
    }
    return numerator / denom;
}
 
double Phase::molarDensity() const
{
    return density()/meanMolecularWeight();
}
 
double Phase::molarVolume() const
{
    return 1.0/molarDensity();
}
 
void Phase::setDensity(const double density_)
{
    assertCompressible("setDensity");
    if (density_ > 0.0) {
        m_dens = density_;
    } else {
        throw CanteraError("Phase::setDensity",
            "density must be positive. density = {}", density_);
    }
}
 
void Phase::assignDensity(const double density_)
{
    if (density_ > 0.0) {
        m_dens = density_;
    } else {
        throw CanteraError("Phase::assignDensity",
            "density must be positive. density = {}", density_);
    }
}
 
double Phase::chargeDensity() const
{
    double cdens = 0.0;
    for (size_t k = 0; k < m_kk; k++) {
        cdens += charge(k)*moleFraction(k);
    }
    return cdens * Faraday;
}
 
double Phase::mean_X(span<const double> Q) const
{
    checkArraySize("Phase::mean_X", Q.size(), m_kk);
    return m_mmw*std::inner_product(m_ym.begin(), m_ym.end(), Q.begin(), 0.0);
}
 
double Phase::sum_xlogx() const
{
    double sumxlogx = 0;
    for (size_t k = 0; k < m_kk; k++) {
        sumxlogx += m_ym[k] * std::log(std::max(m_ym[k], SmallNumber));
    }
    return m_mmw * sumxlogx + std::log(m_mmw);
}
 
size_t Phase::addElement(const string& symbol, double weight, int atomic_number,
                         double entropy298, int elem_type)
{
    // Look up the atomic weight if not given
    if (weight == 0.0) {
        try {
            weight = getElementWeight(symbol);
        } catch (CanteraError&) {
            // assume this is just a custom element with zero atomic weight
        }
    } else if (weight == -12345.0) {
        weight = getElementWeight(symbol);
    }
 
    // Try to look up the standard entropy if not given. Fail silently.
    if (entropy298 == ENTROPY298_UNKNOWN) {
        try {
            const static AnyMap db = AnyMap::fromYamlFile(
                "element-standard-entropies.yaml");
            const AnyMap& elem = db["elements"].getMapWhere("symbol", symbol);
            entropy298 = elem.convert("entropy298", "J/kmol/K", ENTROPY298_UNKNOWN);
        } catch (CanteraError&) {
        }
    }
 
    // Check for duplicates
    auto iter = find(m_elementNames.begin(), m_elementNames.end(), symbol);
    if (iter != m_elementNames.end()) {
        size_t m = iter - m_elementNames.begin();
        if (m_atomicWeights[m] != weight) {
            throw CanteraError("Phase::addElement",
                "Duplicate elements ({}) have different weights", symbol);
        } else {
            // Ignore attempt to add duplicate element with the same weight
            return m;
        }
    }
 
    // Add the new element
    m_atomicWeights.push_back(weight);
    m_elementNames.push_back(symbol);
    m_atomicNumbers.push_back(atomic_number);
    m_entropy298.push_back(entropy298);
    if (symbol == "E") {
        m_elem_type.push_back(CT_ELEM_TYPE_ELECTRONCHARGE);
    } else {
        m_elem_type.push_back(elem_type);
    }
    m_mm++;
 
    // Update species compositions
    if (m_kk) {
        vector<double> old(m_speciesComp);
        m_speciesComp.resize(m_kk*m_mm, 0.0);
        for (size_t k = 0; k < m_kk; k++) {
            size_t m_old = m_mm - 1;
            for (size_t m = 0; m < m_old; m++) {
                m_speciesComp[k * m_mm + m] = old[k * (m_old) + m];
            }
            m_speciesComp[k * (m_mm) + (m_mm-1)] = 0.0;
        }
    }
 
    return m_mm-1;
}
 
vector<AnyMap> Phase::elementDefinitions() const
{
    vector<AnyMap> definitions;
    definitions.reserve(nElements());
    const static AnyMap db = AnyMap::fromYamlFile("element-standard-entropies.yaml");
    const static auto standardEntropies = db["elements"].asMap("symbol");
    const auto& weights = elementWeights();
 
    for (size_t m = 0; m < nElements(); m++) {
        const string& symbol = m_elementNames[m];
        bool exportDefinition = false;
        AnyMap element;
        element["symbol"] = symbol;
 
        auto weightIter = weights.find(symbol);
        if (weightIter == weights.end()) {
            exportDefinition = true;
            element["atomic-weight"] = m_atomicWeights[m];
        } else if (m_atomicWeights[m] != weightIter->second) {
            exportDefinition = true;
            element["atomic-weight"] = m_atomicWeights[m];
        }
 
        if (m_atomicNumbers[m] != 0) {
            exportDefinition = true;
            element["atomic-number"] = m_atomicNumbers[m];
        }
 
        if (m_entropy298[m] != ENTROPY298_UNKNOWN) {
            auto entropyIter = standardEntropies.find(symbol);
            if (entropyIter == standardEntropies.end()
                || !entropyIter->second->hasKey("entropy298"))
            {
                exportDefinition = true;
                element["entropy298"].setQuantity(m_entropy298[m], "J/kmol/K");
            } else {
                double defaultEntropy =
                    entropyIter->second->convert("entropy298", "J/kmol/K");
                if (m_entropy298[m] != defaultEntropy) {
                    exportDefinition = true;
                    element["entropy298"].setQuantity(m_entropy298[m], "J/kmol/K");
                }
            }
        }
 
        if (exportDefinition) {
            definitions.push_back(std::move(element));
        }
    }
    return definitions;
}
 
bool Phase::addSpecies(shared_ptr<Species> spec)
{
    if (m_nSpeciesLocks) {
        throw CanteraError("Phase::addSpecies",
                "Cannot add species to ThermoPhase '{}' because it is being "
                "used by another object,\nsuch as a Reactor, Domain1D (flame), "
                "SolutionArray, or MultiPhase (Mixture) object.", m_name);
    }
 
    if (m_species.find(spec->name) != m_species.end()) {
        throw CanteraError("Phase::addSpecies",
            "Phase '{}' already contains a species named '{}'.",
            m_name, spec->name);
    }
 
    vector<double> comp(nElements());
    for (const auto& [eName, stoich] : spec->composition) {
        size_t m = elementIndex(eName, false);
        if (m == npos) { // Element doesn't exist in this phase
            switch (m_undefinedElementBehavior) {
            case UndefElement::ignore:
                return false;
 
            case UndefElement::add:
                addElement(eName);
                comp.resize(nElements());
                m = elementIndex(eName, true);
                break;
 
            case UndefElement::error:
            default:
                throw CanteraError("Phase::addSpecies",
                    "Species '{}' contains an undefined element '{}'.",
                    spec->name, eName);
            }
        }
        comp[m] = stoich;
    }
 
    size_t ne = nElements();
    auto aw = atomicWeights();
    if (spec->charge != 0.0) {
        size_t eindex = elementIndex("E", false);
        if (eindex != npos) {
            double ecomp = comp[eindex];
            if (fabs(spec->charge + ecomp) > 0.001) {
                if (ecomp != 0.0) {
                    throw CanteraError("Phase::addSpecies",
                                       "Input charge and element E compositions differ "
                                       "for species " + spec->name);
                } else {
                    // Just fix up the element E composition based on the input
                    // species charge
                    comp[eindex] = -spec->charge;
                }
            }
        } else {
            addElement("E", 0.000545, 0, 0.0, CT_ELEM_TYPE_ELECTRONCHARGE);
            ne = nElements();
            eindex = elementIndex("E", true);
            comp.resize(ne);
            comp[ne - 1] = - spec->charge;
        }
    }
 
    double wt = 0.0;
    for (size_t m = 0; m < ne; m++) {
        wt += comp[m] * aw[m];
    }
 
    // Some surface phases may define species representing empty sites
    // that have zero molecular weight. Give them a very small molecular
    // weight to avoid dividing by zero.
    wt = std::max(wt, Tiny);
 
    spec->setMolecularWeight(wt);
 
    m_molwts.push_back(wt);
    m_rmolwts.push_back(1.0/wt);
 
    for (size_t m = 0; m < ne; m++) {
        m_speciesComp.push_back(comp[m]);
    }
 
    m_speciesNames.push_back(spec->name);
    m_species[spec->name] = spec;
    m_speciesIndices[spec->name] = m_kk;
    m_speciesCharge.push_back(spec->charge);
 
    string nLower = toLowerCopy(spec->name);
    if (m_speciesLower.find(nLower) == m_speciesLower.end()) {
        m_speciesLower[nLower] = m_kk;
    } else {
        m_speciesLower[nLower] = npos;
    }
    m_kk++;
 
    // Ensure that the Phase has a valid mass fraction vector that sums to
    // one. We will assume that species 0 has a mass fraction of 1.0 and mass
    // fraction of all other species is 0.0.
    if (m_kk == 1) {
        m_y.push_back(1.0);
        m_ym.push_back(m_rmolwts[0]);
        m_mmw = 1.0 / m_ym[0];
    } else {
        m_y.push_back(0.0);
        m_ym.push_back(0.0);
    }
    m_workS.push_back(0.0);
    invalidateCache();
    return true;
}
 
void Phase::modifySpecies(size_t k, shared_ptr<Species> spec)
{
    if (speciesName(k) != spec->name) {
        throw CanteraError("Phase::modifySpecies",
            "New species name '{}' does not match existing name '{}'",
                           spec->name, speciesName(k));
    }
    const shared_ptr<Species>& old = m_species[spec->name];
    if (spec->composition != old->composition) {
        throw CanteraError("Phase::modifySpecies",
            "New composition for '{}' does not match existing composition",
            spec->name);
    }
    m_species[spec->name] = spec;
    invalidateCache();
}
 
void Phase::addSpeciesAlias(const string& name, const string& alias)
{
    if (speciesIndex(alias, false) != npos) {
        throw CanteraError("Phase::addSpeciesAlias",
            "Invalid alias '{}': species already exists", alias);
    }
    size_t k = speciesIndex(name, false);
    if (k != npos) {
        m_speciesIndices[alias] = k;
    } else {
        throw CanteraError("Phase::addSpeciesAlias",
            "Unable to add alias '{}' "
            "(original species '{}' not found).", alias, name);
    }
}
 
void Phase::removeSpeciesLock()
{
    if (!m_nSpeciesLocks) {
        throw CanteraError("Phase::removeSpeciesLock",
                "ThermoPhase '{}' has no current species locks.", m_name);
    }
    m_nSpeciesLocks--;
}
 
vector<string> Phase::findIsomers(const Composition& compMap) const
{
    vector<string> isomerNames;
 
    for (const auto& [name, species] : m_species) {
        if (species->composition == compMap) {
            isomerNames.emplace_back(name);
        }
    }
 
    return isomerNames;
}
 
vector<string> Phase::findIsomers(const string& comp) const
{
    return findIsomers(parseCompString(comp));
}
 
shared_ptr<Species> Phase::species(const string& name) const
{
    size_t k = speciesIndex(name, true);
    return m_species.at(speciesName(k));
}
 
shared_ptr<Species> Phase::species(size_t k) const
{
    checkSpeciesIndex(k);
    return m_species.at(m_speciesNames[k]);
}
 
void Phase::ignoreUndefinedElements() {
    m_undefinedElementBehavior = UndefElement::ignore;
}
 
void Phase::addUndefinedElements() {
    m_undefinedElementBehavior = UndefElement::add;
}
 
void Phase::throwUndefinedElements() {
    m_undefinedElementBehavior = UndefElement::error;
}
 
bool Phase::ready() const
{
    return (m_kk > 0);
}
 
void Phase::invalidateCache() {
    m_stateNum++;
    m_cache.clear();
}
 
void Phase::setMolecularWeight(const int k, const double mw)
{
    m_molwts[k] = mw;
    m_rmolwts[k] = 1.0/mw;
 
    transform(m_y.begin(), m_y.end(), m_rmolwts.begin(), m_ym.begin(),
              multiplies<double>());
    m_mmw = 1.0 / accumulate(m_ym.begin(), m_ym.end(), 0.0);
}
 
void Phase::compositionChanged() {
    m_stateNum++;
}
 
vector<double> Phase::getCompositionFromMap(const Composition& comp) const
{
    vector<double> X(m_kk);
    for (const auto& [name, value] : comp) {
        size_t loc = speciesIndex(name, true);
        X[loc] = value;
    }
    return X;
}
 
void Phase::massFractionsToMoleFractions(span<const double> Y, span<double> X) const
{
    checkArraySize("Phase::massFractionsToMoleFractions", Y.size(), m_kk);
    checkArraySize("Phase::massFractionsToMoleFractions", X.size(), m_kk);
    double rmmw = 0.0;
    for (size_t k = 0; k != m_kk; ++k) {
        rmmw += Y[k]/m_molwts[k];
    }
    if (rmmw == 0.0) {
        throw CanteraError("Phase::massFractionsToMoleFractions",
                           "no input composition given");
    }
    for (size_t k = 0; k != m_kk; ++k) {
        X[k] = Y[k]/(rmmw*m_molwts[k]);
    }
}
 
void Phase::moleFractionsToMassFractions(span<const double> X, span<double> Y) const
{
    checkArraySize("Phase::moleFractionsToMassFractions", X.size(), m_kk);
    checkArraySize("Phase::moleFractionsToMassFractions", Y.size(), m_kk);
    double mmw = dot(m_molwts.begin(), m_molwts.end(), X.begin());
    if (mmw == 0.0) {
        throw CanteraError("Phase::moleFractionsToMassFractions",
                           "no input composition given");
    }
    double rmmw = 1.0/mmw;
    for (size_t k = 0; k != m_kk; ++k) {
        Y[k] = X[k]*m_molwts[k]*rmmw;
    }
}
 
} // namespace Cantera