GammaLawThermo.cpp

Go to the documentation of this file.

#include "thermo/GammaLawThermo.h"
 
#include "thermo/PlasmaState.h"
 
#include <cmath>
#include <stdexcept>
#include <string>
#include <utility>
#include <vector>
 
namespace rizer
{
 
    // ------------------------------------------------------------------------------ //
    //   GammaLawThermo - constructor                                                 //
    // ------------------------------------------------------------------------------ //
 
    GammaLawThermo::GammaLawThermo(double gamma,
                                   std::vector<double> molecular_weights,
                                   int electron_index)
        : m_gamma(gamma),
          m_W(std::move(molecular_weights)),
          m_electron_index(electron_index)
    {
        if (!(m_gamma > 1.0))
        {
            throw std::invalid_argument(
                "GammaLawThermo: the heavy-pool gamma must be > 1, got " +
                std::to_string(m_gamma));
        }
        if (m_W.empty())
        {
            throw std::invalid_argument("GammaLawThermo: molecular_weights is empty.");
        }
        for (std::size_t k = 0; k < m_W.size(); k++)
        {
            if (!(m_W[k] > 0.0))
            {
                throw std::invalid_argument(
                    "GammaLawThermo: molecular_weights[" +
                    std::to_string(k) +
                    "] must be > 0, got " +
                    std::to_string(m_W[k]));
            }
        }
        if (m_electron_index < -1 || m_electron_index >= static_cast<int>(m_W.size()))
        {
            throw std::invalid_argument(
                "GammaLawThermo: electron_index must be -1 (no electron species) or "
                "a valid species index, got " +
                std::to_string(m_electron_index));
        }
    }
 
    // ------------------------------------------------------------------------------ //
    //   GammaLawThermo - public methods                                              //
    // ------------------------------------------------------------------------------ //
 
    void GammaLawThermo::pressures(double /*r*/, double /*t*/, double rho, double Tg,
                                   double Te, const double *Y, double &p_heavy,
                                   double &p_electron) const
    {
        // PHYSICS.md (G6), one pool at a time: p_H = n_H k_B T_H and
        // p_e = n_e k_B T_e, with n_pool/rho the pool's moles per kg of mixture.
        p_heavy = rho * heavyMolesPerMass(Y) * GasConstant * Tg;
        p_electron = rho * electronMolesPerMass(Y) * GasConstant * Te;
    }
 
    double GammaLawThermo::soundSpeed(double /*r*/, double /*t*/, double /*rho*/,
                                      double Tg, double Te, const double *Y) const
    {
        // PHYSICS.md appendix A.4, frozen two-pool speed: each pool follows its
        // own adiabat, isentropic bulk moduli add, so
        // c = sqrt((gamma_H p_H + gamma_e p_e)/rho). Written on the moles-per-kg
        // basis the partial pressures are built from, rho cancels, which keeps
        // the speed defined on a state whose density has not been set.
        const double moles_electron = electronMolesPerMass(Y);
        const double moles_heavy = heavyMolesPerMass(Y);
        if (!(moles_heavy + moles_electron > 0.0))
        {
            throw std::invalid_argument(
                "GammaLawThermo::soundSpeed: mass fractions sum to zero moles.");
        }
        return std::sqrt(GasConstant * (m_gamma * moles_heavy * Tg +
                                        ElectronGamma * moles_electron * Te));
    }
 
    double GammaLawThermo::heavyEnergy(double Tg, const double *Y) const
    {
        // Per kg of MIXTURE: sum_{s!=e} Y_s * R*Tg/(W_s*(gamma-1)), which is
        // the heavy moles per kg of mixture times R*Tg/(gamma-1).
        return Tg * GasConstant / (m_gamma - 1.0) * heavyMolesPerMass(Y);
    }
 
    double GammaLawThermo::electronCv() const
    {
        if (m_electron_index < 0)
        {
            throw std::invalid_argument(
                "GammaLawThermo::electronCv: no electron species "
                "(electron_index = -1); the electron energy row is not solved.");
        }
        const std::size_t e = static_cast<std::size_t>(m_electron_index);
        return monatomicElectronCv(m_W[e]);
    }
 
    double GammaLawThermo::densityFromPressure(double pressure, double Tg,
                                               double Te, const double *Y) const
    {
        // Inverting PHYSICS.md (G6). pressures() writes each pool's number
        // density as (moles per kg of mixture) times rho, so
        //
        //   p = p_H + p_e = rho * [ sum_{s!=e} (Y_s/W_s) R T_H
        //                           + (Y_e/W_e) R T_e ]
        //
        // and at fixed temperatures and composition the bracket is a constant:
        // p is exactly linear in rho with no offset, so rho = p / k with k the
        // bracket. Evaluating pressures() at rho = 1 IS that bracket, because
        // every term there carries one factor of rho -- which is why the
        // inverse reuses the forward formula rather than restating (G6) a
        // second time and leaving two copies to keep in sync. Exact, not a
        // linearisation: any probe density returns the same p(rho_0)/rho_0.
        double p_heavy_unit = 0.0;
        double p_electron_unit = 0.0;
        pressures(/*r=*/0.0, /*t=*/0.0, 1.0, Tg, Te, Y, p_heavy_unit,
                  p_electron_unit);
        const double p_unit = p_heavy_unit + p_electron_unit;
        if (!(p_unit > 0.0))
        {
            throw std::invalid_argument(
                "GammaLawThermo::densityFromPressure: the state has no pressure at "
                "unit density; check the temperatures and mass fractions.");
        }
        return pressure / p_unit;
    }
 
    double GammaLawThermo::heavyTemperature(double /*rho*/, const double *Y,
                                            double heavy_internal_energy_target,
                                            double /*Tg_seed*/) const
    {
        // heavyEnergy() is linear in Tg for a calorically perfect gas, so the
        // inversion is exact and the seed is unused. A Cantera-backed accessor
        // needs the seed for its Newton iteration.
        const double moles = heavyMolesPerMass(Y);
        if (!(moles > 0.0))
        {
            throw std::invalid_argument(
                "GammaLawThermo::heavyTemperature: no heavy species present.");
        }
        // The stub's energy is R*Tg/(gamma-1) per heavy mole, strictly positive
        // for a positive temperature, so a non-positive target is an unphysical
        // state and not a cold one. Raised rather than returned: the bare
        // division would hand back a negative Tg, and the next soundSpeed() call
        // takes the square root of a negative bulk modulus, so the NaN surfaces
        // several sub-cycles later against the density row instead. A Cantera-backed
        // accessor has its own bracket and does not need this, since formation
        // enthalpies make a negative target legitimate there.
        if (!(heavy_internal_energy_target > 0.0))
        {
            throw std::runtime_error(
                "GammaLawThermo::heavyTemperature: heavy internal energy is " +
                std::to_string(heavy_internal_energy_target) +
                " J/kg, which no positive temperature can produce. The energy "
                "row has lost more to the kinetic split than it holds; lower "
                "cfl or refine the grid.");
        }
        return heavy_internal_energy_target * (m_gamma - 1.0) / (GasConstant * moles);
    }
 
    // ------------------------------------------------------------------------------ //
    //   GammaLawThermo - private methods                                             //
    // ------------------------------------------------------------------------------ //
 
    double GammaLawThermo::heavyMolesPerMass(const double *Y) const
    {
        double moles = 0.0;
        for (std::size_t k = 0; k < m_W.size(); k++)
        {
            if (static_cast<int>(k) == m_electron_index)
                continue;
            moles += Y[k] / m_W[k];
        }
        return moles;
    }
 
    double GammaLawThermo::electronMolesPerMass(const double *Y) const
    {
        if (m_electron_index < 0)
            return 0.0;
        const std::size_t e = static_cast<std::size_t>(m_electron_index);
        return Y[e] / m_W[e];
    }
 
} // namespace rizer