PlasmaChannel1D.cpp

Go to the documentation of this file.

#include "PlasmaChannel1D.h"
#include "PlasmaRates.h"
#include "units.h"
 
#include "cantera/thermo/PlasmaPhase.h"
#include "cantera/kinetics/Kinetics.h"
#include "cantera/kinetics/Reaction.h"
#include "cantera/base/ct_defs.h"
 
#include <algorithm>
#include <cmath>
#include <stdexcept>
 
namespace rizer {
 
// Ideal gas constant on the kmol basis (matches Cantera's molecularWeights(),
// which are in kg/kmol), used by the species-diffusion compressibility term
// below. Mirrors the same constant/convention in ReactorRHS.cpp.
constexpr double GasConstant = units::R_kmol;  // J/kmol/K
 
namespace {
// Register the custom plasma rates (janev / reverse-two-temperature-plasma)
// before loading the mechanism, so YAML files that use them load in C++.
// Idempotent — safe to call on every construction.
Cantera::shared_ptr<Cantera::Solution>
loadMechanism(const std::string& mech, const std::string& phase)
{
    registerPlasmaRates();
    return Cantera::newSolution(mech, phase);
}
} // namespace
 
PlasmaChannel1D::PlasmaChannel1D(const PlasmaChannelConfig& cfg)
    : Cantera::Domain1D(/*nv placeholder*/1, cfg.n_points)
    , m_solution(loadMechanism(cfg.mech, cfg.phase))
    , m_n_species(m_solution->thermo()->nSpecies())
    , m_electron_index(m_solution->thermo()->speciesIndex("e-"))
    , m_outer_radius(cfg.R_max)
    , m_rho(cfg.rho)
    , m_electric_field(cfg.electric_field)
    , m_species_diffusivity(cfg.D_species)
    , m_reacting(cfg.reacting)
    , m_ambient_temperature(cfg.T_amb)
    , m_Y0(cfg.Y0)
    , m_kappa(PropertyTable::fromTableOrScalar(cfg.kappa_T, cfg.kappa_v))
    , m_kappa_e(PropertyTable::fromTableOrScalar(cfg.kappa_e_T, cfg.kappa_e_v))
    , m_initTg(PropertyTable::fromTableOrScalar(cfg.initTg_r, cfg.initTg_v))
    , m_initTe(PropertyTable::fromTableOrScalar(cfg.initTe_r, cfg.initTe_v))
{
    // An explicit grid (if given) defines the mesh, n_points and R_max.
    const bool explicit_grid = cfg.grid.size() >= 2;
    const std::size_t n_points = explicit_grid ? cfg.grid.size() : cfg.n_points;
    if (explicit_grid) {
        m_outer_radius = cfg.grid.back();
    } else if (cfg.R_max <= 0.0 || cfg.n_points < 1) {
        throw std::invalid_argument("PlasmaChannel1D: need R_max>0, n_points>=1 "
                                    "(or supply an explicit grid of >=2 points).");
    }
    if (m_electron_index == Cantera::npos) {
        throw std::invalid_argument("PlasmaChannel1D: mechanism has no 'e-' species.");
    }
    if (m_Y0.size() != m_n_species) {
        throw std::invalid_argument("PlasmaChannel1D: Y0 length != nSpecies.");
    }
    // A composition-resolved collision model (cfg.specs, built from Python's
    // `mtcf` list of MomentumTransferCollisionFrequencyModel) and the legacy
    // sigma/nu_m/nu_E closure are two different ways to get sigma/nu_E into
    // the SAME ReactorRHS::Config.collision slot; mixing them would silently
    // pick one (whichever ReactorRHS::transport() checks first) and ignore
    // the other's inputs, so reject the combination outright instead of
    // guessing which one the caller meant.
    // WHY this matters physically (see ARCHITECTURE.md "Composition-resolved
    // sigma/nu_E"): sigma is far more sensitive to the *ionization state* than
    // to bulk composition. In a fast NRP-style discharge, electron-impact
    // ionization can run orders of magnitude ahead of the channel's initial
    // (often near-neutral) seed composition within a couple of nanoseconds.
    // cfg.specs recomputes sigma/nu_E live from the *evolving* composition
    // every eval() call, tracking that; the legacy sigma(Te)/nu_m(Te) tables
    // are frozen at whatever composition/shape they were tabulated for, and
    // can only re-scale with Te, not with the actual evolving n_e/ionization
    // fraction. Silently combining the two would recombine a frozen sigma(Te)
    // *shape* with the channel's own independently-evolving n_e(t), which is
    // exactly the ~10% error the composition-resolved model was added to fix.
    const bool have_collision = !cfg.specs.empty();
    const bool have_legacy_closure = !cfg.sigma_T.empty() || !cfg.sigma_v.empty()
        || !cfg.nu_m_T.empty() || !cfg.nu_m_v.empty()
        || cfg.nu_m > 0.0 || cfg.nu_E != 0.0;
    if (have_collision && have_legacy_closure) {
        throw std::invalid_argument(
            "PlasmaChannel1D: cfg.specs (composition-resolved collision model) "
            "cannot be combined with sigma_T/sigma_v, nu_m, nu_m_T/nu_m_v, or "
            "nu_E (the legacy closure) -- pick one.");
    }
    // Single-cell 2T physics (chemistry/Joule/exchange) shared with the 0D reactor.
    // The field E is supplied per eval via setElectricField().
    {
        ReactorRHS::Config reactor_config;
        reactor_config.reacting = cfg.reacting;
        if (have_collision) {
            // Same pattern as Plasma0DReactor::Plasma0DReactor: molar_mass and
            // is_electron are derived from the mechanism, not the caller, since
            // they must exactly match this domain's own species indexing.
            std::vector<CollisionModel::SpeciesSpec> specs = cfg.specs;
            specs.resize(m_n_species);
            const auto& molecular_weights = m_solution->thermo()->molecularWeights();
            for (std::size_t k = 0; k < m_n_species; k++) {
                // kg/kmol -> kg/mol
                specs[k].molar_mass = molecular_weights[k] * 1.0e-3;
                specs[k].is_electron = (k == m_electron_index);
            }
            reactor_config.collision =
                CollisionModel(specs, cfg.Te_min, cfg.Te_max, cfg.Te_n, cfg.spitzer);
        } else {
            reactor_config.sigma_T = cfg.sigma_T;  reactor_config.sigma_v = cfg.sigma_v;
            reactor_config.nu_m_T = cfg.nu_m_T;    reactor_config.nu_m_v = cfg.nu_m_v;
            reactor_config.nu_m = cfg.nu_m;        reactor_config.nu_E = cfg.nu_E;
        }
        m_reactor_rhs =
            std::make_unique<ReactorRHS>(m_solution, std::move(reactor_config));
    }
    // Initial radial composition profile: one Tg-style table per species.
    if (!cfg.initY_r.empty() && cfg.initY_r.size() >= 2
        && cfg.initY_v.size() == cfg.initY_r.size() * m_n_species) {
        const std::size_t n_radii = cfg.initY_r.size();
        m_initY.reserve(m_n_species);
        for (std::size_t k = 0; k < m_n_species; k++) {
            std::vector<double> column(n_radii);
            for (std::size_t i = 0; i < n_radii; i++) {
                column[i] = cfg.initY_v[i * m_n_species + k];
            }
            m_initY.push_back(PropertyTable::fromTableOrScalar(cfg.initY_r, column));
        }
    }
    m_has_transport =
        !m_kappa.empty() || !m_kappa_e.empty() || m_species_diffusivity > 0.0;
 
    Cantera::Domain1D::resize(2 + m_n_species, n_points);
    if (explicit_grid) {
        setupGrid(Cantera::span<const double>(cfg.grid.data(), cfg.grid.size()));
    } else if (cfg.n_points >= 2) {
        setupUniformGrid(cfg.n_points, cfg.R_max, 0.0);
    }  // n_points==1: keep the single node at r=0 (resize left m_z={0});
       // setupUniformGrid would divide by (points-1)=0.
 
    setComponentName(0, "Tg");
    setComponentName(1, "Te");
    for (std::size_t k = 0; k < m_n_species; k++) {
        setComponentName(2 + k, m_solution->thermo()->speciesName(k));
    }
    // Newton step-limiting bounds, not a hard clamp (resetBadValues() is the
    // clamp): 200 K keeps both temperatures comfortably above 0 K/the thermo
    // fit floor without constraining genuine cooling; 1e6 K is a generous
    // ceiling no physical run should approach. Species bounds allow a small
    // negative undershoot (-1e-12) since Newton trial states can dip slightly
    // below zero for a trace species without being "bad" -- resetBadValues()
    // clips to [0,1] only after repeated Newton failures, not on every step.
    setBounds(0, 200.0, 1.0e6);
    setBounds(1, 200.0, 1.0e6);
    for (std::size_t k = 0; k < m_n_species; k++) {
        setBounds(2 + k, -1.0e-12, 2.0);
    }
    // Transient absolute tolerance is 100x tighter than steady (1e-12 vs
    // 1e-10), the same ratio ThermalPlasmaColumn1D uses for its own T
    // component -- a deliberate, consistent choice across this codebase's
    // Domain1D subclasses, not an arbitrary pair of numbers.
    setSteadyTolerances(1.0e-4, 1.0e-10);
    setTransientTolerances(1.0e-4, 1.0e-12);
}
 
// Map a Domain1D component index to its name ("Tg", "Te", or a species name),
// the inverse of componentIndex(). Required by Cantera's Domain1D interface
// (used e.g. for solution output, error messages, and by our own callers in
// TransientPlasmaChannelSolver.cpp to enumerate/report components by name).
std::string PlasmaChannel1D::componentName(std::size_t n) const
{
    if (n == 0) return "Tg";
    if (n == 1) return "Te";
    return m_solution->thermo()->speciesName(n - 2);
}
 
// Map a component name back to its Domain1D index (inverse of componentName()).
// The `checkAlias` argument from the Domain1D base signature is unused: this
// domain has no aliases, only the canonical "Tg"/"Te"/species names.
std::size_t PlasmaChannel1D::componentIndex(const std::string& name, bool) const
{
    if (name == "Tg") return 0;
    if (name == "Te") return 1;
    std::size_t k = m_solution->thermo()->speciesIndex(name);
    if (k != Cantera::npos) return 2 + k;
    throw std::invalid_argument("PlasmaChannel1D: no component '" + name + "'.");
}
 
bool PlasmaChannel1D::hasComponent(const std::string& name, bool) const
{
    return name == "Tg" || name == "Te"
           || m_solution->thermo()->speciesIndex(name) != Cantera::npos;
}
 
// Read one component's radial profile out of the domain's own solution buffer
// (m_state, populated by Cantera after getInitialSoln()/a solve). Used by
// TransientPlasmaChannelSolver::solveChannelTransient to seed its flat state
// vector `x` from the freshly-initialized profile (see the comment there on
// why getState() can't be used for that instead).
void PlasmaChannel1D::getValues(const std::string& component,
                                 Cantera::span<double> values) const
{
    if (!m_state) {
        throw std::runtime_error("PlasmaChannel1D::getValues: not installed.");
    }
    std::size_t comp = componentIndex(component);
    const double* soln = m_state->data() + loc();
    for (std::size_t j = 0; j < m_points; j++) {
        values[j] = soln[index(comp, j)];
    }
}
 
// Cantera calls this once per (component, node) to seed the initial guess
// before the first solve. Falls back to a uniform value (m_ambient_temperature
// or the ambient mass fraction m_Y0[k]) wherever no explicit initial-profile
// table was supplied, so a channel with no Tg_profile/Te_profile/Y0_profile
// simply starts uniform.
double PlasmaChannel1D::initialValue(std::size_t n, std::size_t j)
{
    double r = z(j);
    if (n == 0) return m_initTg.empty() ? m_ambient_temperature : m_initTg.eval(r);
    if (n == 1) return m_initTe.empty() ? m_ambient_temperature : m_initTe.eval(r);
    std::size_t k = n - 2;
    return m_initY.empty() ? m_Y0[k] : m_initY[k].eval(r);
}
 
// Clamp a trial state back onto the physically valid manifold after a failed
// Newton/BDF step produces something the thermo layer can't evaluate: floor
// both temperatures at 200 K (Cantera's thermo/kinetics routines are not
// guaranteed valid, and can throw, below their fit range or at/below 0 K) and
// clip every mass fraction into [0, 1] (a valid composition; Newton overshoots
// can otherwise drive Y_k slightly negative or above 1). This does not
// conserve mass Sum(Y_k)=1 -- it is a last-resort rescue to get an evaluable
// state for the next retry, not a physical correction.
void PlasmaChannel1D::resetBadValues(Cantera::span<double> xg)
{
    auto x = xg.subspan(loc(), size());
    for (std::size_t j = 0; j < m_points; j++) {
        x[index(0, j)] = std::max(x[index(0, j)], 200.0);
        x[index(1, j)] = std::max(x[index(1, j)], 200.0);
        for (std::size_t k = 0; k < m_n_species; k++) {
            x[index(2 + k, j)] = std::min(std::max(x[index(2 + k, j)], 0.0), 1.0);
        }
    }
}
 
// Thin pass-through to the shared 0D physics (ReactorRHS), so every caller
// (eval()'s residual, and TransientPlasmaChannelSolver's history recording)
// computes n_e the same way the reaction source terms themselves see it.
double PlasmaChannel1D::electronDensity(double Tg, double Te, const double* Y) const
{
    return m_reactor_rhs->electronDensity(Tg, Te, Y, m_rho);
}
 
// Cantera's Jacobian coloring calls eval() once per perturbed global point
// `jg` (finite-difference Jacobian: only that point's residual needs
// recomputing) plus once with jg==npos (a full residual evaluation, e.g. for
// the initial guess or a converged-solution check). This domain's stencil
// couples each node only to its immediate east/west neighbours (see fvDiv
// below), so perturbing point `jg` can only change the residual at jg-1, jg,
// and jg+1 -- hence the narrow [jmin, jmax] window computed below instead of
// recomputing every node on every perturbation.
void PlasmaChannel1D::eval(std::size_t jg, Cantera::span<const double> xg,
                            Cantera::span<double> rg, Cantera::span<int> maskg,
                            double rdt)
{
    // Points outside this domain's own range (jg is a *global* index across
    // all domains in the Sim1D) can't affect our residual at all; skip early.
    if (jg != Cantera::npos
        && (jg + 1 < firstPoint() || jg > lastPoint() + 1)) {
        return;
    }
    auto x = xg.subspan(loc(), size());
    auto rsd = rg.subspan(loc(), size());
    auto diag = maskg.subspan(loc(), size());
 
    std::size_t jmin, jmax;
    if (jg == Cantera::npos) {
        // Full evaluation: every node's residual needs (re)computing.
        jmin = 0;
        jmax = m_points - 1;
    } else {
        // Narrow window: only the perturbed point and its two neighbours.
        std::size_t local_point = (jg == 0) ? 0 : jg - firstPoint();
        jmin = std::max<std::size_t>(local_point, 1) - 1;
        jmax = std::min(local_point + 1, m_points - 1);
    }
 
    // Finite-volume cylindrical divergence (1/r) d/dr(r * coef * dphi/dr) at
    // node `node` for component `component`, face coefficient from coef(phi).
    // The axis face (node==0) has zero area (symmetry). The outer wall is
    // imposed as a Dirichlet ghost node held at `ghost`, a half-interval
    // beyond the last node -- this keeps every node transient (well scaled),
    // avoiding an ill-conditioned algebraic boundary row.
    // Units below: radius/spacing/volume are always [m]/[m]/[m^2] (geometry is
    // component-independent); `value_*`/`coef_*`/`flux_*`/the return value take
    // whatever unit `component` has -- [K] and [W/m/K] for the Tg/Te conduction
    // calls, [-] and [kg/m/s] for the species-diffusion call below (see the two
    // call sites: coef is kappa/kappa_e or rho*D respectively).
    // Lambda parameters (a plain, non-Doxygen comment: this lambda has no
    // separate declaration for \param to attach to, so a doc-comment here
    // would be misattached to the enclosing eval()):
    //   component  Component index being diffused: 0=Tg, 1=Te, or 2+k for
    //              species k [-]
    //   node       Radial grid node index [-]
    //   ghost      Dirichlet ghost value beyond the wall, in `component`'s
    //              unit ([K] for Tg/Te, [-] for Y_k)
    //   coef       Face coefficient functor value->coef(value): kappa/kappa_e
    //              [W/m/K] for Tg/Te, or rho*D [kg/m/s] for species
    auto fvDiv = [&](std::size_t component, std::size_t node, double ghost,
                      auto coef) -> double {
        double value_here = x[index(component, node)];
        double radius_east, spacing_east, value_east;  // [m], [m], component's unit
        if (node < m_points - 1) {
            value_east = x[index(component, node + 1)];
            radius_east = 0.5 * (z(node) + z(node + 1));
            spacing_east = z(node + 1) - z(node);
        } else {  // last node: ghost at r = z + spacing_west, value = ambient
            // [m]
            double last_spacing =
                (m_points >= 2) ? z(node) - z(node - 1) : m_outer_radius;
            value_east = ghost;
            radius_east = z(node) + 0.5 * last_spacing;
            spacing_east = last_spacing;
        }
        double coef_east = 0.5 * (coef(value_here) + coef(value_east));
        double flux_east =
            radius_east * coef_east * (value_east - value_here) / spacing_east;
        double flux_west = 0.0, radius_west = 0.0;  // [m]
        if (node > 0) {
            double value_west = x[index(component, node - 1)];
            radius_west = 0.5 * (z(node - 1) + z(node));
            double spacing_west = z(node) - z(node - 1);  // [m]
            double coef_west = 0.5 * (coef(value_west) + coef(value_here));
            flux_west =
                radius_west * coef_west * (value_here - value_west) / spacing_west;
        }
        // [m^2]
        double volume =
            0.5 * (radius_east * radius_east - radius_west * radius_west);
        return (flux_east - flux_west) / volume;
    };
 
    // Node-centered radial gradient d(phi)/dr at `node` for `component`,
    // used by the species-diffusion enthalpy-transport term below (a mixed-
    // gradient PRODUCT term, not a flux divergence, so fvDiv() itself does
    // not apply here). Obtained as the arithmetic mean of the same one-sided
    // east/west face differences fvDiv() already forms internally, so the
    // gradient used here stays consistent with the divergence stencil
    // discretizing the rest of the transport terms. At the axis (node==0)
    // this returns 0 exactly -- symmetry makes every component an even
    // function of r there, so its radial derivative vanishes -- rather than
    // a one-sided difference extrapolated across the (zero-area) axis face.
    // At the wall (node==m_points-1) the east-side difference uses the same
    // Dirichlet ghost value `ghost` a half-cell beyond the last node that
    // fvDiv()'s own wall treatment uses.
    //   component  Component index: 0=Tg, 1=Te, or 2+k for species k [-]
    //   node       Radial grid node index [-]
    //   ghost      Dirichlet ghost value beyond the wall, in `component`'s
    //              unit ([K] for Tg/Te, [-] for Y_k)
    auto centeredGradient = [&](std::size_t component, std::size_t node,
                                 double ghost) -> double {
        if (node == 0) return 0.0;  // axis symmetry: d(phi)/dr(r=0) = 0
        double value_here = x[index(component, node)];
        double value_west = x[index(component, node - 1)];
        double spacing_west = z(node) - z(node - 1);  // [m]
        double gradient_west = (value_here - value_west) / spacing_west;
 
        double gradient_east;
        if (node < m_points - 1) {
            double value_east = x[index(component, node + 1)];
            double spacing_east = z(node + 1) - z(node);  // [m]
            gradient_east = (value_east - value_here) / spacing_east;
        } else {  // wall: ghost a half-cell beyond the last node
            double last_spacing =
                (m_points >= 2) ? z(node) - z(node - 1) : m_outer_radius;
            gradient_east = (ghost - value_here) / last_spacing;
        }
        return 0.5 * (gradient_west + gradient_east);
    };
 
    // Molecular weights W_k [kg/kmol], node-independent; fetched once and
    // reused every node by the species-diffusion compressibility term below.
    const std::vector<double>& molecular_weights = m_reactor_rhs->molecularWeights();
 
    std::vector<double> dYdt(m_n_species);  // dY_k/dt [1/s], one entry per species
    std::vector<double> raw_Y_divergence(m_n_species);  // see the species-diffusion
                                                         // block below [1/m^2]
    for (std::size_t j = jmin; j <= jmax; j++) {
        double Tg = x[index(0, j)];  // gas temperature [K]
        double Te = x[index(1, j)];  // electron temperature [K]
        const double* Y = &x[index(2, j)];  // species mass fractions [-]
 
        // lhsTg/lhsTe are the effective heat capacities per unit volume
        // [J/m^3/K] that ReactorRHS::rates already divided the chemistry/
        // Joule/exchange source terms by, so dTg/dTe below are already true
        // dT/dt [K/s] for those terms. The conduction term computed via
        // fvDiv() is a bare flux divergence, not yet divided by capacitance,
        // so it must be divided by the SAME lhsTg/lhsTe here to add
        // consistently onto dTg/dTe (see the class-level header comment for
        // the governing PDEs, rho*cv_h*dTg/dt = div(...) + ...).
        // [K/s], [K/s], [J/m^3/K], [J/m^3/K]
        double dTg = 0.0, dTe = 0.0, lhsTg = 1.0, lhsTe = 1.0;
        std::fill(dYdt.begin(), dYdt.end(), 0.0);
        m_reactor_rhs->rates(Tg, Te, Y, m_rho, m_electric_field, dTg, dTe,
                              dYdt.data(), lhsTg, lhsTe);
 
        if (m_has_transport) {
            if (!m_kappa.empty()) {
                dTg += fvDiv(0, j, m_ambient_temperature,
                             [&](double T) { return m_kappa.eval(T); }) / lhsTg;
            }
            if (!m_kappa_e.empty()) {
                dTe += fvDiv(1, j, m_ambient_temperature,
                             [&](double T) { return m_kappa_e.eval(T); }) / lhsTe;
            }
            if (m_species_diffusivity > 0.0) {
                // Raw (unit-coefficient) divergence div(grad(Y_k)) =
                // (1/r) d/dr(r dY_k/dr) at this node, for every species k --
                // computed once per species and reused below for the
                // ordinary Fickian diffusion term (Term A) AND the species-
                // diffusion compressibility term (Term B). Term A's actual
                // face coefficient rho*D is a spatially uniform constant
                // (does not depend on Y_k), so fvDiv(coef=rho*D) equals
                // rho*D*raw_Y_divergence[k] exactly -- no need to call
                // fvDiv() twice per species with two different constant
                // coefficients.
                for (std::size_t k = 0; k < m_n_species; k++) {
                    raw_Y_divergence[k] =
                        fvDiv(2 + k, j, m_Y0[k], [](double) { return 1.0; });
                }
 
                // --- Term A: ordinary Fickian species diffusion,
                // div(rho*D*grad(Y_k)), converted to dY_k/dt by /rho ==
                // D*raw_Y_divergence[k] (see above; rho assumed spatially
                // uniform -- this model has no continuity equation).
                // Matches the mass-fraction species equation in the
                // class-level header comment.
                for (std::size_t k = 0; k < m_n_species; k++) {
                    dYdt[k] += m_species_diffusivity * raw_Y_divergence[k];
                }
 
                // --- Term B: species-diffusion compressibility. Expanding
                // DP/Dt for the two-temperature ideal-gas EOS under Fick's
                // law contributes an extra energy-equation source (see
                // ARCHITECTURE.md, "species-diffusion compressibility"):
                //   electron: +rho*D*(R*Te/We) * (1/r) d/dr(r dYe/dr)
                //   gas: +rho*D*R*Tg * sum_{k!=e} (1/Wk) * (1/r) d/dr(r dYk/dr)
                // (summed over heavy species only -- they share one Tg).
                const double We = molecular_weights[m_electron_index];
                dTe += m_rho * m_species_diffusivity * (GasConstant * Te / We)
                       * raw_Y_divergence[m_electron_index] / lhsTe;
 
                double gas_compressibility_sum = 0.0;
                for (std::size_t k = 0; k < m_n_species; k++) {
                    if (k == m_electron_index) continue;
                    gas_compressibility_sum +=
                        raw_Y_divergence[k] / molecular_weights[k];
                }
                dTg += m_rho * m_species_diffusivity * GasConstant * Tg
                       * gas_compressibility_sum / lhsTg;
 
                // --- Term C: enthalpy diffusion (species heat transport).
                // The full multicomponent energy equation's heat flux
                // includes +sum_k h_k*j_k (species carry their own sensible
                // enthalpy as they diffuse); with Fick's law j_k=-rho*D*
                // grad(Y_k), its divergence -div(sum_k h_k*j_k) reduces to a
                // mixed-gradient PRODUCT term rather than a divergence (see
                // ARCHITECTURE.md, "enthalpy diffusion"):
                //   electron: +rho*D*cp_e * (dYe/dr)(dTe/dr)
                //   gas: +rho*D*(dTg/dr) * sum_{k!=e} cp_k * (dYk/dr)
                const std::vector<double>& species_heat_capacity =
                    m_reactor_rhs->speciesCp();
                const double Tg_gradient =
                    centeredGradient(0, j, m_ambient_temperature);
                const double Te_gradient =
                    centeredGradient(1, j, m_ambient_temperature);
                const double Ye_gradient = centeredGradient(
                    2 + m_electron_index, j, m_Y0[m_electron_index]);
                dTe += m_rho * m_species_diffusivity
                       * species_heat_capacity[m_electron_index] * Ye_gradient
                       * Te_gradient / lhsTe;
 
                double gas_enthalpy_diffusion_sum = 0.0;
                for (std::size_t k = 0; k < m_n_species; k++) {
                    if (k == m_electron_index) continue;
                    double Yk_gradient = centeredGradient(2 + k, j, m_Y0[k]);
                    gas_enthalpy_diffusion_sum +=
                        species_heat_capacity[k] * Yk_gradient;
                }
                dTg += m_rho * m_species_diffusivity * Tg_gradient
                       * gas_enthalpy_diffusion_sum / lhsTg;
            }
        }
 
        // Steady residual = source+transport rate; the -rdt*(phi-phi_prev) term
        // is Cantera's pseudo-transient damping (rdt = 1/dt during a time step,
        // 0 during a steady solve or when computing a bare method-of-lines RHS
        // -- see ChannelRHS::eval in TransientPlasmaChannelSolver.cpp). diag=1
        // marks every component here as a true (non-algebraic) unknown, so
        // Cantera includes the transient term for all of them; there is no
        // Dirichlet/algebraic row in this domain (the wall is a ghost-cell
        // flux, not a fixed unknown).
        rsd[index(0, j)] = dTg - rdt * (Tg - prevSoln(0, j));
        rsd[index(1, j)] = dTe - rdt * (Te - prevSoln(1, j));
        diag[index(0, j)] = 1;
        diag[index(1, j)] = 1;
        for (std::size_t k = 0; k < m_n_species; k++) {
            rsd[index(2 + k, j)] = dYdt[k] - rdt * (Y[k] - prevSoln(2 + k, j));
            diag[index(2 + k, j)] = 1;
        }
    }
}
 
} // namespace rizer