CanteraThermo.cpp

Go to the documentation of this file.

#include "thermo/CanteraThermo.h"
 
#include "thermo/PlasmaState.h"
 
#include "cantera/base/global.h"
#include "cantera/thermo/ThermoFactory.h"
 
#include <algorithm>
#include <cmath>
#include <stdexcept>
#include <string>
 
namespace rizer
{
 
    // ------------------------------------------------------------------------------ //
    //   Constructor                                                                  //
    // ------------------------------------------------------------------------------ //
 
    CanteraThermo::CanteraThermo(const std::string &mech, const std::string &phase)
    {
        m_phase = Cantera::newThermo(mech, phase);
        m_thermo = m_phase.get();
        // Resolved once, by dynamic_cast: every two-temperature call below is a
        // PlasmaPhase member, and a plain IdealGasPhase would carry no electron
        // temperature at all. A static_cast would reach them anyway and read a
        // phase that does not have them.
        m_plasma = dynamic_cast<Cantera::PlasmaPhase *>(m_thermo);
        if (m_plasma == nullptr)
        {
            throw std::invalid_argument(
                "CanteraThermo: phase \"" + phase + "\" in \"" + mech +
                "\" is a " + m_thermo->type() +
                ", not a PlasmaPhase, so it carries no electron temperature. "
                "Declare it with \"thermo: plasma\" in the mechanism.");
        }
        m_nsp = m_thermo->nSpecies();
        m_ie = m_plasma->electronSpeciesIndex();
        m_W_electron = m_thermo->molecularWeight(m_ie);
        m_Yc.assign(m_nsp, 0.0);
        m_work.assign(m_nsp, 0.0);
    }
 
    // ------------------------------------------------------------------------------ //
    //   Public methods                                                               //
    // ------------------------------------------------------------------------------ //
 
    void CanteraThermo::pressures(double /*r*/, double /*t*/, double rho, double Tg,
                                  double Te, const double *Y, double &p_heavy,
                                  double &p_electron) const
    {
        setState(rho, Tg, Te, Y);
        // electronPressure() is n_e k_B T_e and pressure() is n_tot k_B T_mean
        // with T_mean the mole-fraction-weighted mean, so their difference is
        // n_H k_B T_H exactly: the split of PHYSICS.md (G6), taken from the
        // phase rather than recomputed from it.
        const TwoTemperatureState eos = readTwoTemperatureState(*m_plasma);
        p_electron = eos.electron_pressure;
        p_heavy = eos.heavy_pressure;
    }
 
    double CanteraThermo::soundSpeed(double /*r*/, double /*t*/, double rho,
                                     double Tg, double Te, const double *Y) const
    {
        // PHYSICS.md appendix A.4, frozen two-pool speed: the pools' isentropic
        // bulk moduli add, so c = sqrt((gamma_H p_H + gamma_e p_e)/rho) with the
        // same partial-pressure split pressures() uses.
        //
        // gamma_H is summed over the heavy species alone rather than taken as
        // the mixture's cp_mass/cv_mass. The electron's heat capacity per kg of
        // ELECTRONS is enormous (1.5 R / W_e), so even at Y_e ~ 1e-7 it carries
        // a percent of the mixture cp and the mixture ratio is not the heavy
        // pool's. The electron pool enters through ElectronGamma instead.
        double p_heavy = 0.0;
        double p_electron = 0.0;
        pressures(0.0, 0.0, rho, Tg, Te, Y, p_heavy, p_electron);
        return soundSpeedAtSetState(p_heavy, p_electron, Tg);
    }
 
    double CanteraThermo::heavyEnergy(double Tg, const double *Y) const
    {
        // An ideal gas's internal energy is set by temperature and composition
        // alone, so the density only has to keep the state valid, and the
        // electron temperature cancels: PlasmaPhase evaluates the electron term
        // at whatever Te is set and the subtraction below takes the same term
        // back out, which is why Tg is passed for both.
        setState(/*rho=*/1.0, Tg, /*Te=*/Tg, Y);
        return heavyEnergyAtSetState();
    }
 
    double CanteraThermo::electronCv() const
    {
        return monatomicElectronCv(m_W_electron);
    }
 
    double CanteraThermo::densityFromPressure(double pressure, double Tg,
                                              double Te, const double *Y) const
    {
        // PlasmaPhase::setPressure() is the inverse of the pressure() that
        // pressures() splits: it sets rho = p Wbar / (R T_mean) with the same
        // mean temperature, so there is nothing here to invert by hand.
        setState(/*rho=*/1.0, Tg, Te, Y);
        m_thermo->setPressure(pressure);
        const double density = m_thermo->density();
        if (!(density > 0.0))
        {
            throw std::invalid_argument(
                "CanteraThermo::densityFromPressure: a target pressure of " +
                std::to_string(pressure) + " Pa gives a density of " +
                std::to_string(density) +
                " kg/m^3; the pressure, the temperatures or the mass fractions "
                "do not describe a gas.");
        }
        return density;
    }
 
    double CanteraThermo::heavyTemperature(double /*rho*/, const double *Y,
                                           double heavy_internal_energy_target,
                                           double Tg_seed) const
    {
        // Newton on e_H(Tg) = target. e_H is monotone in Tg (every c_v > 0), so
        // the only way to leave the bracket is a bad step; damping to the bracket
        // rather than trusting the step keeps a stale seed from diverging.
        constexpr double T_low = TemperatureFloor;
        constexpr double T_high = TemperatureCeiling;
        constexpr int max_iterations = 100;
        constexpr double tolerance = 1.0e-10;
 
        if (!std::isfinite(Tg_seed))
        {
            throw std::runtime_error(
                "CanteraThermo::heavyTemperature: the Tg seed is " +
                std::to_string(Tg_seed) + "; a non-finite seed means the state "
                "handed in is already corrupted.");
        }
        double Tg = std::min(std::max(Tg_seed, T_low), T_high);
 
        for (int iteration = 0; iteration < max_iterations; iteration++)
        {
            double internal_energy, slope; // [J/kg of mixture], [J/kg/K]
            heavyEnergyAndCv(Tg, Y, internal_energy, slope);
            const double residual = internal_energy - heavy_internal_energy_target;
            if (!(slope > 0.0))
            {
                throw std::runtime_error(
                    "CanteraThermo::heavyTemperature: non-positive heavy heat "
                    "capacity at Tg = " +
                    std::to_string(Tg) + " K.");
            }
            const double step = residual / slope;
            const double Tg_next = std::min(std::max(Tg - step, T_low), T_high);
            if (std::abs(Tg_next - Tg) <= tolerance * std::max(1.0, Tg_next))
            {
                // A step that stops moving is only convergence if the residual went
                // with it. Sitting against a bracket edge also stops moving -- the
                // clamp returns the same value every iteration - and reporting that
                // as a converged temperature would hand the caller the floor
                // dressed up as an answer.
                const double scale =
                    std::max(std::abs(heavy_internal_energy_target), 1.0);
                const double final_residual =
                    heavyEnergy(Tg_next, Y) - heavy_internal_energy_target;
                if (std::abs(final_residual) > 1.0e-6 * scale)
                {
                    throw std::runtime_error(
                        "CanteraThermo::heavyTemperature: stalled at Tg = " +
                        std::to_string(Tg_next) + " K with a relative energy residual of " +
                        std::to_string(std::abs(final_residual) / scale) +
                        ". The target internal energy of " +
                        std::to_string(heavy_internal_energy_target) +
                        " J/kg is outside the bracket [" + std::to_string(T_low) + ", " +
                        std::to_string(T_high) + "] K for this composition.");
                }
                return Tg_next;
            }
            Tg = Tg_next;
        }
        throw std::runtime_error(
            "CanteraThermo::heavyTemperature: Newton did not converge in " +
            std::to_string(max_iterations) + " iterations from seed " +
            std::to_string(Tg_seed) + " K.");
    }
 
    StateProperties CanteraThermo::stateProperties(double /*r*/, double /*t*/,
                                                   double rho, double Tg, double Te,
                                                   const double *Y) const
    {
        setState(rho, Tg, Te, Y);
        const TwoTemperatureState eos = readTwoTemperatureState(*m_plasma);
        StateProperties properties;
        properties.p_heavy = eos.heavy_pressure;
        properties.p_electron = eos.electron_pressure;
        properties.sound_speed =
            soundSpeedAtSetState(eos.heavy_pressure, eos.electron_pressure, Tg);
        properties.heavy_energy = heavyEnergyAtSetState();
        return properties;
    }
 
    // ------------------------------------------------------------------------------ //
    //   Private methods                                                              //
    // ------------------------------------------------------------------------------ //
 
    void CanteraThermo::setState(double rho, double Tg, double Te,
                                 const double *Y) const
    {
        // Rejected rather than floored. Cantera's NASA polynomials are fitted
        // from 200 K up and a hydrodynamic state that arrives colder is broken,
        // not merely cold: flooring it would hand every caller below a pressure
        // and an energy for a state the solver never asked about. The one place
        // a sub-floor temperature is legitimate is the heavy-pool inversion,
        // and heavyTemperature() keeps its own iterates inside the bracket.
        if (!(Tg >= TemperatureFloor) || !(Te >= TemperatureFloor) ||
            !(Tg <= TemperatureCeiling))
        {
            throw std::runtime_error(
                "CanteraThermo::setState: Tg = " + std::to_string(Tg) +
                " K, Te = " + std::to_string(Te) + " K; Tg must lie within [" +
                std::to_string(TemperatureFloor) + ", " +
                std::to_string(TemperatureCeiling) +
                "] K and Te must be at least " + std::to_string(TemperatureFloor) +
                " K, where the mechanism's thermo fits hold.");
        }
        // Negative mass fractions are clipped, as PlasmaCell::setState clips
        // them (see there for why, and for the measured sizes): the hydro step
        // reads Y after the reaction sub-step, and measured inputs reach -3.5e-8
        // (C2H6+) on the Python split solver. A NaN is a broken state, not
        // round-off, so it raises instead of being clipped to zero.
        for (std::size_t k = 0; k < m_nsp; k++)
        {
            if (std::isnan(Y[k]))
            {
                throw std::runtime_error(
                    "CanteraThermo::setState: the mass fraction of species " +
                    m_thermo->speciesName(k) + " is NaN at Tg = " +
                    std::to_string(Tg) + " K, Te = " + std::to_string(Te) + " K.");
            }
            m_Yc[k] = (Y[k] > 0.0) ? Y[k] : 0.0;
        }
        m_thermo->setMassFractions_NoNorm(
            Cantera::span<const double>(m_Yc.data(), m_nsp));
        m_plasma->setElectronTemperature(Te);
        m_thermo->setState_TD(Tg, rho);
    }
 
    void CanteraThermo::heavyEnergyAndCv(double Tg, const double *Y,
                                         double &internal_energy,
                                         double &cv) const
    {
        // heavyEnergy() leaves the state it set on the phase, which is the state
        // heavyCpAndCv() then reads: one state-setting for both, and the slope
        // is taken at the same state as the energy it is the slope of.
        internal_energy = heavyEnergy(Tg, Y);
        double cp_heavy = 0.0; // [J/kg of mixture/K], not needed here
        heavyCpAndCv(cp_heavy, cv);
    }
 
    void CanteraThermo::heavyCpAndCv(double &cp, double &cv) const
    {
        // Mayer per species, c_v = (c_p^molar - R) / W, summed over the heavy
        // species alone. Cantera offers no heavy-pool-only capacity: cp_mass()
        // carries the electron at its own temperature.
        m_thermo->getPartialMolarCp(Cantera::span<double>(m_work.data(), m_nsp));
        cp = 0.0;
        cv = 0.0;
        for (std::size_t k = 0; k < m_nsp; k++)
        {
            if (k == m_ie)
                continue;
            const double moles =
                m_thermo->massFraction(k) / m_thermo->molecularWeight(k);
            cp += moles * m_work[k];
            cv += moles * (m_work[k] - Cantera::GasConstant);
        }
    }
 
    double CanteraThermo::soundSpeedAtSetState(double p_heavy, double p_electron,
                                               double Tg) const
    {
        double cp_heavy = 0.0; // [J/kg of mixture/K]
        double cv_heavy = 0.0; // [J/kg of mixture/K]
        heavyCpAndCv(cp_heavy, cv_heavy);
        if (!(cv_heavy > 0.0))
        {
            throw std::runtime_error(
                "CanteraThermo::soundSpeed: non-positive heavy heat capacity at "
                "Tg = " +
                std::to_string(Tg) + " K; the mixture has no heavy pool to carry "
                                     "an acoustic wave.");
        }
        return std::sqrt(((cp_heavy / cv_heavy) * p_heavy +
                          ElectronGamma * p_electron) /
                         m_thermo->density());
    }
 
    double CanteraThermo::heavyEnergyAtSetState() const
    {
        // intEnergy_mass() is sum_k Y_k u_k / W_k over the whole mixture, with
        // PlasmaPhase evaluating the electron's u at Te and every heavy
        // species' at Tg. Removing the electron term leaves the heavy pool on
        // the per kg of MIXTURE basis that PHYSICS.md (4b) conserves.
        m_thermo->getPartialMolarIntEnergies(
            Cantera::span<double>(m_work.data(), m_nsp));
        return m_thermo->intEnergy_mass() -
               m_thermo->massFraction(m_ie) * m_work[m_ie] / m_W_electron;
    }
 
} // namespace rizer