PlasmaChannel1D.cpp

Go to the documentation of this file.

#include "PlasmaChannel1D.h"
 
#include "models/nrp/RadialDiffusion.h"
#include "kinetics/PlasmaRates.h"
#include "numeric/DiffusiveFlux.h"
#include "numeric/Format.h"
#include "thermo/PlasmaState.h"
#include "units.h"
 
#include "cantera/thermo/PlasmaPhase.h"
#include "cantera/base/ct_defs.h"
 
#ifdef _OPENMP
#include <omp.h>
#endif
 
#include <algorithm>
#include <cmath>
#include <numeric>
#include <sstream>
#include <stdexcept>
 
namespace rizer
{
 
    // ------------------------------------------------------------------------------ //
    //   Constructor                                                                  //
    // ------------------------------------------------------------------------------ //
 
    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_electric_field(cfg.electric_field),
          m_species_diffusivity(cfg.D_species),
          m_reacting(cfg.reacting),
          m_hydrodynamics(cfg.hydrodynamics),
          m_ambient_temperature(cfg.T_amb),
          m_Y_amb(cfg.Y_amb),
          m_kappa_H(PropertyTable::fromTableOrScalar(
              cfg.kappa_H_T, cfg.kappa_H_v, "kappa_H_T/kappa_H_v")),
          m_kappa_e(PropertyTable::fromTableOrScalar(
              cfg.kappa_e_T, cfg.kappa_e_v, "kappa_e_T/kappa_e_v")),
          m_initTg(PropertyTable::fromTableOrScalar(
              cfg.initTg_r, cfg.initTg_v, "initTg_r/initTg_v")),
          m_initTe(PropertyTable::fromTableOrScalar(
              cfg.initTe_r, cfg.initTe_v, "initTe_r/initTe_v"))
    {
        // Eight steps, in order:
        //   1. resolve the mesh, then the inputs every mode needs
        //   2. pick the transport closure and build the shared PlasmaCell
        //   3. the initial composition profile, one table per species
        //   4. size the domain, install the grid, require electrons at every node
        //   5. the radial density profile, sampled once the mesh exists
        //   6. the hydro-only inputs, accepted or refused as a set
        //   7. the wall ghost density, from everything validated above
        //   8. solver tolerances
 
        // ---- 1. mesh, then the inputs every mode needs ----------------------
        // 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)
        {
            // faceRadius(0) is pinned at the axis, so a first node off it would
            // make the first cell the whole disc inside that node.
            if (cfg.grid.front() != 0.0)
            {
                throw std::invalid_argument(
                    "PlasmaChannel1D: an explicit grid must start on the axis, "
                    "grid[0] = 0, but grid[0] = " +
                    std::to_string(cfg.grid.front()) +
                    " m. The first face is pinned at r = 0, so the first cell "
                    "would span the whole disc inside that node.");
            }
            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.");
        }
        // Required outright, never derived from the initial composition: the
        // two are independent, and the one the caller forgot is the one a
        // fallback would hide (see PlasmaChannelConfig::Y_amb).
        if (m_Y_amb.empty())
        {
            throw std::invalid_argument(
                "PlasmaChannel1D: Y_amb is required. It is the composition "
                "OUTSIDE the domain, which the wall ghost and the hydro inflow "
                "read; it is not derived from the initial composition. Pass the "
                "ambient composition explicitly, repeating the initial one only "
                "if that is genuinely what the boundary holds.");
        }
        if (m_Y_amb.size() != m_n_species)
        {
            throw std::invalid_argument(
                "PlasmaChannel1D: Y_amb has " +
                std::to_string(m_Y_amb.size()) +
                " entries, but the mechanism has " +
                std::to_string(m_n_species) +
                " species.");
        }
        // Checked whatever the mode: the species-diffusion wall ghost reads
        // Y_amb as well as the hydro boundary, and a ghost that does not sum
        // to 1 injects or removes mass through the wall face.
        const double ambient_sum =
            std::accumulate(m_Y_amb.begin(), m_Y_amb.end(), 0.0);
        if (!(std::abs(ambient_sum - 1.0) <= 1.0e-10))
        {
            throw std::invalid_argument(
                "PlasmaChannel1D: Y_amb, the composition of the outer boundary, "
                "sums to " +
                std::to_string(ambient_sum) + ", not 1 to within 1e-10.");
        }
 
        // ---- 2. transport closure, then the shared cell ---------------------
        // 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 PlasmaCellConfig.collision slot; mixing them would silently
        // pick one (whichever PlasmaCell::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 PHYSICS.md section 4.1, the
        // 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().
        {
            PlasmaCellConfig reactor_config;
            reactor_config.reacting = cfg.reacting;
            reactor_config.species_h0k = cfg.species_h0k;
            if (have_collision)
            {
                // Same pattern as PlasmaReactor0D::PlasmaReactor0D: molar_mass and
                // is_electron are derived from the mechanism, not the caller, since
                // they must exactly match this domain's own species indexing.
                // Refused rather than resized: padding gives the missing
                // species no collision model and truncating drops the caller's
                // tail, so either way a species silently leaves every
                // collision sum.
                if (cfg.specs.size() != m_n_species)
                {
                    throw std::invalid_argument(
                        "PlasmaChannel1D: got collision specs for " +
                        std::to_string(cfg.specs.size()) +
                        " species, but the mechanism has " +
                        std::to_string(m_n_species) +
                        " species. Give one spec per mechanism species, in its "
                        "order.");
                }
                std::vector<CollisionModel::SpeciesSpec> specs = cfg.specs;
                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);
            }
            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;
            }
            // Kept so setEvaluatorPoolSize() can build further evaluators with the
            // identical closure rather than re-deriving it from the config.
            m_reactor_config = reactor_config;
            m_mech = cfg.mech;
            m_phase = cfg.phase;
            m_cell =
                std::make_unique<PlasmaCell>(m_solution, std::move(reactor_config));
        }
 
        // ---- 3. initial composition profile ---------------------------------
        // Initial radial composition profile: one Tg-style table per species.
        // Required, like initTg/initTe -- a uniform composition is this table
        // with two radii, which is what the Python layer builds from `Y0`.
        if (cfg.initY_r.size() < 2 ||
            cfg.initY_v.size() != cfg.initY_r.size() * m_n_species)
        {
            throw std::invalid_argument(
                "PlasmaChannel1D: initY_r needs at least 2 radii and initY_v "
                "exactly initY_r.size() * nSpecies entries, row-major; got " +
                std::to_string(cfg.initY_r.size()) + " and " +
                std::to_string(cfg.initY_v.size()) + " for " +
                std::to_string(m_n_species) + " 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, "initY_r/initY_v"));
            }
        }
        m_has_transport =
            !m_kappa_H.empty() || !m_kappa_e.empty() || m_species_diffusivity > 0.0;
 
        // ---- 4. size the domain, install the grid, require electrons -------
        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.
 
        // Every initial profile has to span the grid: outside its radii a table
        // holds its end value, so a profile ending short of the wall would start
        // every node past it at the last radius's value.
        const double outer_node = z(n_points - 1);
        auto requireSpansGrid = [&](const PropertyTable &table, const std::string &name)
        {
            if (!table.empty() && !(table.covers(0.0) && table.covers(outer_node)))
            {
                throw std::invalid_argument(
                    "PlasmaChannel1D: the initial profile " + name + " covers [" +
                    significant(table.first()) + ", " + significant(table.last()) +
                    "] m but the grid spans [0, " + significant(outer_node) + "] m.");
            }
        };
        requireSpansGrid(m_initTg, "initTg_r/initTg_v");
        requireSpansGrid(m_initTe, "initTe_r/initTe_v");
        for (const PropertyTable &table : m_initY)
        {
            requireSpansGrid(table, "initY_r/initY_v");
        }
 
        // Every node has to start with electrons. Te is the temperature OF the
        // electrons, so in a cell holding none it is undefined, and electrons
        // that chemistry later creates there would start at whatever value the
        // cell carried. Checked at the nodes, once they exist, because the
        // composition is a radial table and a zero can sit between its radii.
        // The ambient is exempt: it only feeds the wall ghost, and the cells
        // next to it keep electrons through the face fluxes.
        for (std::size_t j = 0; j < n_points; j++)
        {
            const double Y_electron = m_initY[m_electron_index].eval(z(j));
            if (!(Y_electron > 0.0))
            {
                std::ostringstream message;
                message << "PlasmaChannel1D: the initial electron mass fraction is "
                        << Y_electron << " at node " << j << " (r=" << z(j)
                        << " m). The electron temperature is undefined where "
                           "there are no electrons; seed a small electron "
                           "fraction at every node, with a matching positive "
                           "ion to keep the charge balance.";
                throw std::invalid_argument(message.str());
            }
        }
 
        // A mismatched pair would make fromTableOrScalar() hand back an empty
        // table, leaving the scalar rho at every node -- a silently different
        // physical problem, and the one NUMERICS.md section "Validation ladder" warns reads as a
        // solver bug rather than as the setup error it is.
        const bool want_density_profile =
            !cfg.init_rho_r.empty() || !cfg.init_rho_v.empty();
        if (want_density_profile && (cfg.init_rho_r.size() < 2 ||
                                     cfg.init_rho_r.size() != cfg.init_rho_v.size()))
        {
            throw std::invalid_argument(
                "PlasmaChannel1D: init_rho_r and init_rho_v must be the same "
                "length and have at least 2 points. For a density that is the "
                "same everywhere, leave both empty and set the scalar `rho`.");
        }
 
        // ---- 5. radial density profile --------------------------------------
        // Sampled after the mesh exists, since the profile is a function of radius.
        const PropertyTable init_rho =
            PropertyTable::fromTableOrScalar(cfg.init_rho_r, cfg.init_rho_v,
                                             "init_rho_r/init_rho_v");
        requireSpansGrid(init_rho, "init_rho_r/init_rho_v");
        m_rho.assign(n_points, cfg.rho);
        if (!init_rho.empty())
        {
            for (std::size_t j = 0; j < n_points; j++)
            {
                m_rho[j] = init_rho.eval(z(j));
            }
        }
        for (std::size_t j = 0; j < n_points; j++)
        {
            if (!(m_rho[j] > 0.0))
            {
                std::ostringstream message;
                message
                    << "PlasmaChannel1D: density at node "
                    << j << " (r=" << z(j) << " m) is " << m_rho[j]
                    << " kg/m^3. The species equation divides by it, so a zero or "
                       "negative density is not a state this domain can evaluate.";
                throw std::invalid_argument(message.str());
            }
        }
 
        // ---- 6. hydro-only inputs, accepted or refused as a set -------------
        // Parsed whether or not the sub-step runs: an unknown name is a typo in
        // either mode, and a mode the caller chose that is never looked at is the
        // same silent drop as one that is misspelled.
        m_outer_bc_mode = parseOuterBc(cfg.outer_bc_mode);
 
        // Nothing outside the hydro sub-step reads these, so carrying them would
        // report a configuration the run never applied -- an initial velocity most
        // of all, which would come back in the recorded history as a flow this
        // mode does not have.
        if (!m_hydrodynamics)
        {
            std::string unused;
            auto name = [&unused](const std::string &field)
            {
                unused += unused.empty() ? field : ", " + field;
            };
            if (!cfg.init_ur_r.empty() || !cfg.init_ur_v.empty())
            {
                name("init_ur_r/init_ur_v");
            }
            if (cfg.p_amb >= 0.0)
                name("p_amb");
            if (cfg.outer_bc_mode != "ambient_ghost")
                name("outer_bc_mode");
            if (!unused.empty())
            {
                throw std::invalid_argument(
                    "PlasmaChannel1D: " +
                    unused +
                    " is read only by the hydrodynamic sub-step, which this "
                    "config leaves off. Set hydrodynamics=true, or drop it.");
            }
        }
 
        // The hydro sub-step treats a zero p_amb as a vacuum, a state a
        // defaulted config could reach without the caller ever choosing it.
        if (m_hydrodynamics)
        {
            if (!(cfg.p_amb > 0.0))
            {
                throw std::invalid_argument(
                    "PlasmaChannel1D: hydrodynamics needs p_amb > 0 Pa, the "
                    "ambient pressure its outer boundary imposes.");
            }
            if (!(cfg.cfl > 0.0) || cfg.cfl > 1.0)
            {
                throw std::invalid_argument(
                    "PlasmaChannel1D: hydrodynamics needs 0 < cfl <= 1.");
            }
            if (n_points < 2)
            {
                throw std::invalid_argument(
                    "PlasmaChannel1D: hydrodynamics needs at least 2 radial "
                    "nodes; a single node has no face to carry a flux.");
            }
        }
 
        // ---- 7. wall ghost density ------------------------------------------
        // Below every check above, because it consumes what they validate: Y0 is
        // handed to setMassFractions_NoNorm, which normalizes nothing, so a sum
        // of zero would reach the equation of state as an infinite mixture weight
        // and come back as an infinite density or a Cantera error naming the
        // phase instead of the field the caller set.
        //
        // The wall ghost is already pinned to the ambient composition and the
        // ambient temperature, so taking its density from the outermost node would
        // build one face coefficient out of two different states. With an ambient
        // pressure in hand the equation of state closes it properly; with
        // hydrodynamics off there is no ambient pressure to close it with, and a
        // zero-gradient extrapolation is the only thing the config supports.
        m_wall_ghost_density = m_rho.back();
        if (m_hydrodynamics)
        {
            auto thermo = m_solution->thermo();
            thermo->setMassFractions_NoNorm(
                Cantera::span<const double>(m_Y_amb.data(), m_Y_amb.size()));
            // Both pools at the ambient temperature, set explicitly. setState_TP
            // leaves Te untouched, and PlasmaPhase closes the density from the
            // MEAN of the two, so without this line the ambient ghost would depend
            // on whatever electron temperature the phase happened to carry -- the
            // same (p, T_amb, T_amb, Y_amb) inversion the hydro sub-step's ambient
            // ghost performs through ThermoAccessor::densityFromPressure().
            thermo->setElectronTemperature(cfg.T_amb);
            thermo->setState_TP(cfg.T_amb, cfg.p_amb);
            m_wall_ghost_density = thermo->density();
        }
 
        const PropertyTable init_ur =
            PropertyTable::fromTableOrScalar(cfg.init_ur_r, cfg.init_ur_v,
                                             "init_ur_r/init_ur_v");
        requireSpansGrid(init_ur, "init_ur_r/init_ur_v");
        m_u_r.assign(n_points, 0.0);
        if (!init_ur.empty())
        {
            for (std::size_t j = 0; j < n_points; j++)
            {
                m_u_r[j] = init_ur.eval(z(j));
            }
        }
 
        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));
        }
 
        // ---- 8. solver tolerances -------------------------------------------
        // 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(/*rtol=*/1.0e-4, /*atol=*/1.0e-10);
        setTransientTolerances(/*rtol=*/1.0e-4, /*atol=*/1.0e-12);
    }
 
    // ------------------------------------------------------------------------------ //
    //   Public methods                                                               //
    // ------------------------------------------------------------------------------ //
 
    // 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
    // RadialDiffusion.h), 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)
    {
        // Five steps, in order:
        //   1. is this global point ours at all, and which slice of it
        //   2. the node window the perturbation can reach
        //   3. the transport coefficients and ghosts, assembled once
        //   4. per node: source terms, then the radial coupling on top
        //   5. per node: residual = rate - rdt*(phi - phi_prev), mask all ones
 
        // ---- 1. is this point ours ------------------------------------------
        // 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());
 
        // ---- 2. the node window to recompute --------------------------------
        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);
        }
 
        // ---- 3. node-independent coefficients, assembled once ---------------
        // Assembled once and reused for every node of the window: the transport
        // coefficients and ghosts do not vary with j.
        const RadialDiffusionCoefficients diffusion = diffusionCoefficients();
 
        std::vector<double> dYdt(m_n_species); // dY_k/dt [1/s], one entry per species
        // div(rho*D*grad(Y_k)) [kg/m^3/s], the species-diffusion block's Fickian
        // flux divergence, reused by Terms A and B below
        std::vector<double> species_flux_divergence(m_n_species);
 
        // ---- 4. per node: source terms, then the radial coupling ------------
        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 PlasmaCell::sourceTerms 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
            // faceDivergence() 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_cell->sourceTerms(Tg, Te, Y, m_rho[j], m_electric_field,
                                dTg, dTe, dYdt.data(), lhsTg, lhsTe);
 
            if (m_has_transport)
            {
                addRadialDiffusion(*this, x, diffusion, j, lhsTg, lhsTe, dTg, dTe,
                                   dYdt.data(), species_flux_divergence.data());
            }
 
            // ---- 5. per node: residual and transient mask -------------------
            // 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 PlasmaChannelSolver.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;
            }
        }
    }
 
    // 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
    // PlasmaChannelSolver.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
    // PlasmaChannelSolver::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. The temperatures fall back to a uniform
    // m_ambient_temperature where no profile table was supplied; the
    // composition has no fallback, since initY is required (a uniform
    // composition is that table with two radii).
    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[k].eval(r);
    }
 
    // The single writer of the two fields the hydro sub-step owns, so the
    // positivity the species equation needs is checked in one place.
    void PlasmaChannel1D::setHydroState(const std::vector<double> &rho,
                                        const std::vector<double> &u_r)
    {
        if (rho.size() != m_points || u_r.size() != m_points)
        {
            throw std::invalid_argument(
                "PlasmaChannel1D::setHydroState: rho and u_r must each have one "
                "entry per radial node.");
        }
        for (std::size_t j = 0; j < m_points; j++)
        {
            if (!(rho[j] > 0.0))
            {
                std::ostringstream message;
                message << "PlasmaChannel1D::setHydroState: density at node " << j
                        << " (r=" << z(j) << " m) is " << rho[j]
                        << " kg/m^3. The hydro sub-step has driven it to zero or "
                           "below; lower cfl or refine the grid.";
                throw std::invalid_argument(message.str());
            }
        }
        m_rho = rho;
        m_u_r = u_r;
    }
 
    // Thin pass-through to the shared 0D physics (PlasmaCell), so the recorded
    // history reports n_e and p the way the reaction source terms themselves see
    // them. The channel carries no pressure of its own: the state is
    // (rho, u_r, Tg, Te, Y_s) and p follows from the equation of state, so both
    // of these are diagnostic reads through the same state-setting the residual
    // uses -- one state-setting, which is why they are returned together.
    void PlasmaChannel1D::nodeDiagnostics(double Tg, double Te, const double *Y,
                                          std::size_t j, double &pressure,
                                          double &electron_density) const
    {
        m_cell->pressureAndElectronDensity(Tg, Te, Y, m_rho[j], pressure,
                                           electron_density);
    }
 
    void PlasmaChannel1D::nodeRates(std::size_t j, double Tg, double Te,
                                    const double *Y, double electric_field,
                                    double *ydot) const
    {
        if (m_has_transport)
        {
            throw std::runtime_error(
                "PlasmaChannel1D::nodeRates: this domain has conduction or "
                "species diffusion enabled, so a node's rates depend on its "
                "neighbours and cannot be evaluated one node at a time. Use "
                "eval() instead.");
        }
        // lhs_Tg/lhs_Te are the capacitances rates() has already divided its own
        // sources by, so dTg/dTe come back as true dT/dt. eval() only keeps them
        // because its conduction term is a bare flux divergence that still needs
        // dividing by the same capacitance; with no transport there is nothing
        // left to divide.
        double lhs_Tg = 1.0, lhs_Te = 1.0; // [J/m^3/K]
 
        // Which evaluator this call may touch. Outside a parallel region, and in
        // a build without OpenMP, this is always slot 0.
#ifdef _OPENMP
        const std::size_t slot = static_cast<std::size_t>(omp_get_thread_num());
        // Refused rather than wrapped into range. A slot shared by two threads is
        // two of them rewriting one PlasmaCell's scratch and one Cantera phase
        // state, which returns plausible wrong numbers instead of crashing --
        // exactly what the pool exists to prevent. The caller sizes the pool to
        // the team it is about to start, so reaching here means that invariant
        // broke and the answer cannot be trusted.
        if (slot >= evaluatorPoolSize())
        {
            throw std::runtime_error(
                "PlasmaChannel1D::nodeRates: thread " +
                std::to_string(slot) +
                " has no evaluator of its own; the pool holds " +
                std::to_string(evaluatorPoolSize()) +
                ". Call setEvaluatorPoolSize() with the number of threads that "
                "will enter this function before starting the parallel region.");
        }
#else
        const std::size_t slot = 0;
#endif
        const PlasmaCell &reactor =
            (slot == 0) ? *m_cell : *m_extra_rhs[slot - 1];
        reactor.sourceTerms(Tg, Te, Y, m_rho[j], electric_field,
                            ydot[0], ydot[1], ydot + 2, lhs_Tg, lhs_Te);
    }
 
    double PlasmaChannel1D::nodeElectricalConductivity(double Tg, double Te,
                                                       const double *Y,
                                                       std::size_t j) const
    {
        return m_cell->electricalConductivity(Tg, Te, Y, m_rho[j]);
    }
 
    void PlasmaChannel1D::wallEnergyOutflow(Cantera::span<const double> x,
                                            double &conduction_out,
                                            double &diffusion_out,
                                            std::vector<double> &species_mass_out) const
    {
        conduction_out = 0.0;
        diffusion_out = 0.0;
        species_mass_out.assign(m_n_species, 0.0);
        if (!m_has_transport)
        {
            return;
        }
        const std::size_t wall_node = m_points - 1;
        std::vector<double> species_enthalpy(m_n_species); // [J/kg]
        m_cell->speciesEnthalpies(x[index(0, wall_node)], x[index(1, wall_node)],
                                  &x[index(2, wall_node)], m_rho[wall_node],
                                  species_enthalpy.data());
        const WallEnergyOutflow outflow = rizer::wallEnergyOutflow(
            *this, x, diffusionCoefficients(), species_enthalpy.data());
        conduction_out = outflow.conduction;
        diffusion_out = outflow.diffusion;
        if (!outflow.species_mass.empty())
        {
            species_mass_out = outflow.species_mass;
        }
    }
 
    double PlasmaChannel1D::faceRadius(std::size_t j) const
    {
        if (j == 0)
        {
            return 0.0;
        }
        if (j < m_points)
        {
            return 0.5 * (z(j - 1) + z(j));
        }
        // The wall face, mirrored half a spacing past the last node. A single-point
        // grid has no spacing to mirror, so the domain radius stands in for it,
        // through faceSpacing() like every other face.
        return z(m_points - 1) + 0.5 * faceSpacing(m_points);
    }
    double PlasmaChannel1D::faceSpacing(std::size_t j) const
    {
        if (j < m_points)
        {
            return z(j) - z(j - 1);
        }
        return (m_points >= 2) ? z(m_points - 1) - z(m_points - 2) : m_outer_radius;
    }
 
    void PlasmaChannel1D::setEvaluatorPoolSize(std::size_t slots)
    {
        if (slots < 1)
        {
            throw std::runtime_error(
                "PlasmaChannel1D::setEvaluatorPoolSize: the pool must hold at "
                "least the evaluator built at construction.");
        }
        m_extra_rhs.resize(slots - 1);
        for (auto &evaluator : m_extra_rhs)
        {
            if (!evaluator)
            {
                evaluator = std::make_unique<PlasmaCell>(
                    loadMechanism(m_mech, m_phase), m_reactor_config);
            }
        }
    }
 
    // ------------------------------------------------------------------------------ //
    //   Private methods                                                              //
    // ------------------------------------------------------------------------------ //
 
    RadialDiffusionCoefficients PlasmaChannel1D::diffusionCoefficients() const
    {
        RadialDiffusionCoefficients diffusion;
        diffusion.kappa = &m_kappa_H;
        diffusion.kappa_e = &m_kappa_e;
        diffusion.species_diffusivity = m_species_diffusivity;
        diffusion.ambient_temperature = m_ambient_temperature;
        diffusion.wall_ghost_density = m_wall_ghost_density;
        diffusion.Y_amb = &m_Y_amb;
        diffusion.molecular_weights = &m_cell->molecularWeights();
        diffusion.species_cp = &m_cell->speciesCp();
        diffusion.electron_index = m_electron_index;
        diffusion.n_species = m_n_species;
        return diffusion;
    }
 
} // namespace rizer