bindings.cpp

Go to the documentation of this file.

// pybind11 bindings for the rizer Cantera 1D plasma extension (_plasma1d).
//
// This extension is a pure C++ implementation of the Elenbaas-Heller column
// problem, using Cantera's 1D simulation framework (Sim1D) to solve the
// steady-state radial energy balance with a user-specified axial electric field.
// This file contains the pybind11 bindings, which are a thin wrapper around the
// C++ implementation in PlasmaColumnSolver.h/cpp.
//
// Design rule: "flat numerics only" across the Python boundary.
//   - No Cantera C++ type (Solution, Domain1D, Sim1D, ...) is ever exposed to
//     Python.  This sidesteps all ABI-sharing concerns with the stock `cantera`
//     Python package; both can be loaded in the same process without conflicts.
//   - Properties enter as parallel (T, value) numpy arrays; they are converted
//     to PropertyTable objects inside the binding lambda.
//   - The converged profile leaves as a plain Python dict of numpy arrays and
//     Python scalars — no custom C++ types cross the boundary.
//
// Exposed API:
//   cantera_version() -> str          Version of the linked libcantera
//   solve_column(...) -> dict         Solve and return (r, T, sigma, kappa,
//                                     current, electric_field, n_points)
 
#ifndef NOMINMAX
#define NOMINMAX
#endif
#include <pybind11/pybind11.h>
#include <pybind11/stl.h>
#include <pybind11/numpy.h>
 
#include "cantera/base/global.h"
 
#include "numeric/PropertyTable.h"
#include "transport/CollisionModel.h"
#include "models/nrp/Plasma0DReactor.h"
#include "models/steady_lte_arc/PlasmaColumnSolver.h"
#include "kinetics/PlasmaRates.h"
#include "models/nrp/TransientPlasmaChannelSolver.h"
 
#include <vector>
 
namespace py = pybind11;
 
namespace {

std::vector<double> toVec(const py::array_t<double>& a)
{
    // Re-wrapping as c_style|forcecast forces a contiguous double array (copying
    // only if `a` isn't already one), so the raw-pointer range below is always
    // valid -- unlike raw pointer arithmetic, this does not assume `a` itself is
    // contiguous (e.g. a strided view).
    const py::array_t<double, py::array::c_style | py::array::forcecast> c(a);
    return std::vector<double>(c.data(), c.data() + c.size());
}

rizer::PropertyTable makeTable(const py::array_t<double>& T,
                                const py::array_t<double>& values)
{
    return rizer::PropertyTable(toVec(T), toVec(values));
}
 
} // namespace
 
PYBIND11_MODULE(_plasma1d, m)
{
    m.doc() = "rizer Cantera 1D plasma extension (ThermalPlasmaColumn1D)";
 
    m.def("cantera_version", []() { return Cantera::version(); },
          "Version string of the linked libcantera.");
 
    m.def("register_plasma_rates", &rizer::registerPlasmaRates,
          "(Re-)register the native C++ plasma reaction rates "
          "(janev-dissociative-recombination-C2Hy, Druyvesteyn, "
          "reverse-two-temperature-plasma) with Cantera's ReactionRateFactory. "
          "Idempotent; replaces the Python ExtensibleRate registrations of the "
          "same names, so mechanisms load and evaluate natively afterwards.");
 
    // Maxwellian momentum-transfer kernel Qbar(Te) (Mitchner II-6.30), exposed for
    // validation against the Python TabulatedSpeciesCrossSection integral. Vectorized
    // over Te. energy_J/sigma_m2 are the tabulated cross-section (zero outside range).
    // 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 PYBIND11_MODULE function):
    //   energy_J  Tabulated collision energy grid [J], strictly increasing
    //   sigma_m2  Tabulated cross section Q(energy_J) [m^2], same length as
    //             energy_J
    //   Te        Electron temperatures to evaluate Qbar(Te) at [K]
    //   x_max     Upper limit of the reduced-energy (E/kTe) integration [-]
    //   N         Number of quadrature points for the integration [-]
    m.def("mean_cross_section",
        [](py::array_t<double> energy_J, py::array_t<double> sigma_m2,
           py::array_t<double> Te, double x_max, int N) {
            std::vector<double> e = toVec(energy_J);
            std::vector<double> s = toVec(sigma_m2);
            std::vector<double> t = toVec(Te);
            std::vector<double> q =
                rizer::CollisionModel::meanCrossSectionTable(e, s, t, x_max, N);
            return py::array_t<double>(q.size(), q.data());
        },
        py::arg("energy_J"), py::arg("sigma_m2"), py::arg("Te"),
        py::arg("x_max") = 20.0, py::arg("N") = 1000,
        "Maxwellian-averaged momentum-transfer cross section Qbar(Te) [m^2] "
        "(Mitchner II-6.30), interpolating Q from (energy_J [J], sigma_m2 [m^2]).");
 
    // Native C++ two-temperature 0D plasma reactor: a persistent evaluator that the
    // Python driver (scipy + NRP circuit) calls per RHS. Flat numerics only; the
    // per-species cross-section model is passed once at construction as ragged lists.
    py::class_<rizer::Plasma0DReactor>(m, "Reactor0D")
        // Construct a Reactor0D from a Cantera mechanism and an optional
        // per-species collision model. Lambda parameters (plain comment --
        // see the note at mean_cross_section above for why this isn't `//!`):
        //   mech           Cantera mechanism (YAML) path
        //   phase          Cantera mechanism (YAML) phase
        //   reacting       true: include finite-rate chemistry source terms;
        //                  false: frozen composition
        //   xsec_energy_J  Per-species tabulated cross-section energy grids
        //                  [J] (list of 1-D arrays, or None per species),
        //                  length nsp
        //   xsec_sigma_m2  Per-species tabulated cross-section values [m^2]
        //                  (list of 1-D arrays, or None per species),
        //                  length nsp
        //   radius_m       Per-species hard-sphere effective radius [m]
        //                  (NaN if not used), length nsp
        //   ion_Z          Per-species ion charge number [-] (0 for
        //                  neutrals), length nsp
        //   Te_min         Collision-model Te grid lower bound [K]
        //   Te_max         Collision-model Te grid upper bound [K]
        //   Te_n           Collision-model Te grid point count [-]
        //                  (log-spaced)
        //   spitzer        true: apply the Spitzer e-e correction to the
        //                  collision model
        //   gap            Inter-electrode gap [m] (optional; used in the
        //                  polytropic volume equation)
        //   species_h0k    Per-species standard enthalpy of formation at 0 K
        //                  [J/kmol], length nsp; required when `reacting` and
        //                  the mechanism has an electron-impact reaction (see
        //                  rizer.kinetics.electron_reactions)
        .def(py::init([](const std::string& mech, const std::string& phase,
                         bool reacting, py::list xsec_energy_J,
                         py::list xsec_sigma_m2, py::array_t<double> radius_m,
                         py::array_t<int> ion_Z, double Te_min, double Te_max,
                         int Te_n, bool spitzer, double gap,
                         py::array_t<double> species_h0k) {
                rizer::Plasma0DReactor::Config cfg;
                cfg.mech = mech; cfg.phase = phase; cfg.reacting = reacting;
                cfg.Te_min = Te_min; cfg.Te_max = Te_max;
                cfg.Te_n = static_cast<std::size_t>(std::max(2, Te_n));
                cfg.spitzer = spitzer; cfg.gap = gap;
                cfg.species_h0k = toVec(species_h0k);
                const std::size_t nsp = static_cast<std::size_t>(xsec_energy_J.size());
                cfg.specs.resize(nsp);
                auto r = radius_m.unchecked<1>();
                auto z = ion_Z.unchecked<1>();
                for (std::size_t k = 0; k < nsp; k++) {
                    cfg.specs[k].radius = r(static_cast<py::ssize_t>(k));
                    cfg.specs[k].Z = z(static_cast<py::ssize_t>(k));
                    py::object eo = xsec_energy_J[k];
                    py::object so = xsec_sigma_m2[k];
                    if (!eo.is_none() && !so.is_none()) {
                        cfg.specs[k].energy_J = toVec(eo.cast<py::array_t<double>>());
                        cfg.specs[k].sigma_m2 = toVec(so.cast<py::array_t<double>>());
                    }
                }
                return std::make_unique<rizer::Plasma0DReactor>(cfg);
            }),
            py::arg("mech"), py::arg("phase"), py::arg("reacting") = true,
            py::arg("xsec_energy_J"), py::arg("xsec_sigma_m2"),
            py::arg("radius_m"), py::arg("ion_Z"),
            py::arg("Te_min") = 300.0, py::arg("Te_max") = 1.0e5,
            py::arg("Te_n") = 400, py::arg("spitzer") = true, py::arg("gap") = 0.0,
            py::arg("species_h0k") = py::array_t<double>(0))
        // Electrical conductivity [S/m] at (Tg, Te, Y, rho). Lambda parameters
        // (plain comment, see the note at mean_cross_section above):
        //   self  The bound Reactor0D instance
        //   Tg    Gas (heavy-species) temperature [K]
        //   Te    Electron temperature [K]
        //   Y     Species mass fractions [-], length nsp
        //   rho   Mass density [kg/m^3]
        .def("conductivity",
            [](rizer::Plasma0DReactor& self, double Tg, double Te,
               py::array_t<double> Y, double rho) {
                auto y = toVec(Y);
                return self.conductivity(Tg, Te, y.data(), rho);
            },
            py::arg("Tg"), py::arg("Te"), py::arg("Y"), py::arg("rho"),
            "Electrical conductivity [S/m] at (Tg, Te, Y, rho).")
        // Electron number density [1/m^3] at (Tg, Te, Y, rho). Lambda
        // parameters (plain comment, see the note at mean_cross_section above):
        //   self  The bound Reactor0D instance
        //   Tg    Gas (heavy-species) temperature [K]
        //   Te    Electron temperature [K]
        //   Y     Species mass fractions [-], length nsp
        //   rho   Mass density [kg/m^3]
        .def("electron_density",
            [](rizer::Plasma0DReactor& self, double Tg, double Te,
               py::array_t<double> Y, double rho) {
                auto y = toVec(Y);
                return self.electronDensity(Tg, Te, y.data(), rho);
            },
            py::arg("Tg"), py::arg("Te"), py::arg("Y"), py::arg("rho"),
            "Electron number density [1/m^3] at (Tg, Te, Y, rho).")
        // 0D RHS dydt = [dTg, dTe, dV, dY...] for state [Tg, Te, V, Y...];
        // rho = mass/V; polytropic_index=inf -> constant volume. Lambda
        // parameters (plain comment, see the note at mean_cross_section above):
        //   self              The bound Reactor0D instance
        //   Tg                Gas (heavy-species) temperature [K]
        //   Te                Electron temperature [K]
        //   V                 Volume [m^3]
        //   Y                 Species mass fractions [-], length nsp
        //   E                 Applied electric field [V/m]
        //   mass              Fixed total mass [kg] (rho = mass/V)
        //   p_ext             External pressure [Pa] (drives the polytropic
        //                     volume equation)
        //   polytropic_index  Polytropic index [-] (inf = constant volume)
        .def("rhs",
            [](rizer::Plasma0DReactor& self, double Tg, double Te, double V,
               py::array_t<double> Y, double E, double mass, double p_ext,
               double polytropic_index) {
                auto y = toVec(Y);
                std::vector<double> dydt(3 + self.nSpecies(), 0.0);
                self.rhs(Tg, Te, V, y.data(), E, mass, p_ext, polytropic_index,
                         dydt.data());
                return py::array_t<double>(dydt.size(), dydt.data());
            },
            py::arg("Tg"), py::arg("Te"), py::arg("V"), py::arg("Y"), py::arg("E"),
            py::arg("mass"), py::arg("p_ext"), py::arg("polytropic_index"),
            "0D RHS dydt = [dTg, dTe, dV, dY...] for state [Tg, Te, V, Y...]; "
            "rho = mass/V; polytropic_index=inf -> constant volume.")
        .def_property_readonly("species",
            [](rizer::Plasma0DReactor& self) { return self.speciesNames(); })
        // Power / Maxwellian-condition breakdown from the most recent rhs()
        // call (rizer.models.nrp.isomass_2T_volume_reactor_cpp.Isomass2TVolumeReactor's
        // stored P_elastic/P_inelastic/P_chemical/P_chemical_e/P_Joule/
        // nu_eH/nu_ee/nu_eI/nu_eH_mass_weighted attributes, natively).
        .def_property_readonly("P_elastic",
            [](rizer::Plasma0DReactor& self) { return self.diagnostics().P_elastic; })
        .def_property_readonly("P_inelastic",
            [](rizer::Plasma0DReactor& self) { return self.diagnostics().P_inelastic; })
        .def_property_readonly("P_chemical",
            [](rizer::Plasma0DReactor& self) { return self.diagnostics().P_chemical; })
        .def_property_readonly("P_chemical_e",
            [](rizer::Plasma0DReactor& self) { return self.diagnostics().P_chemical_e; })
        .def_property_readonly("P_Joule",
            [](rizer::Plasma0DReactor& self) { return self.diagnostics().P_Joule; })
        .def_property_readonly("nu_eH",
            [](rizer::Plasma0DReactor& self) { return self.diagnostics().nu_eH; })
        .def_property_readonly("nu_ee",
            [](rizer::Plasma0DReactor& self) { return self.diagnostics().nu_ee; })
        .def_property_readonly("nu_eI",
            [](rizer::Plasma0DReactor& self) { return self.diagnostics().nu_eI; })
        .def_property_readonly("nu_eH_mass_weighted",
            [](rizer::Plasma0DReactor& self) { return self.diagnostics().nu_eH_mass_weighted; });
 
    // Define the solve_column() binding.
    // The lambda converts the Python inputs into C++ types, calls solveColumn(),
    // and converts the result back into a Python dict of numpy arrays and scalars.
    // Lambda parameters (plain comment, see the note at mean_cross_section
    // above -- this lambda has no separate declaration for \param to attach
    // to):
    //   R               Outer wall radius [m]
    //   electric_field  Axial E field [V/m]
    //   T_wall          Wall temperature [K]
    //   sigma_T         Temperature values [K], 1-D array, strictly increasing
    //   sigma_v         Corresponding sigma values [S/m], 1-D array
    //   kappa_T         Temperature values [K], 1-D array, strictly increasing
    //   kappa_v         Corresponding kappa values [W/(m*K)], 1-D array
    //   prad_T          Temperature values [K], 1-D array, strictly increasing,
    //                   or None
    //   prad_v          Corresponding P_rad values [W/m^3], 1-D array, or None
    //   n_points        Initial uniform grid count
    //   T_center_guess  Parabolic-seed centerline temperature [K]
    //   rho_cp          rho*cp [J/m^3/K] for pseudo-transient scaling
    //   loglevel        Sim1D verbosity (0 = silent)
    //   refine_grid     Enable Cantera adaptive grid refinement
    //   refine_ratio    Refiner: max ratio of adjacent cell widths
    //   refine_slope    Refiner: max fractional slope change per cell
    //   refine_curve    Refiner: max fractional curvature per cell
    //   refine_prune    Refiner: remove points below this threshold
    //   init_r          Optional seed radii [m], strictly increasing, 0..R, or
    //                   None
    //   init_T          Optional seed temperatures [K] at init_r, or None
    m.def("solve_column",
        [](double R, double electric_field, double T_wall, py::array_t<double> sigma_T,
           py::array_t<double> sigma_v, py::array_t<double> kappa_T,
           py::array_t<double> kappa_v, py::object prad_T, py::object prad_v,
           std::size_t n_points, double T_center_guess, double rho_cp, int loglevel,
           bool refine_grid, double refine_ratio, double refine_slope,
           double refine_curve, double refine_prune, py::object init_r,
           py::object init_T)
        {
            // Convert the input numpy arrays into PropertyTable objects.
            rizer::PropertyTable sigma = makeTable(sigma_T, sigma_v);
            rizer::PropertyTable kappa = makeTable(kappa_T, kappa_v);
            // For the optional radiation property, check if the inputs are None.
            // If either is None, create a constant PropertyTable with value 0.0.
            // Otherwise, create a PropertyTable from the provided arrays.
            rizer::PropertyTable prad =
                (prad_T.is_none() || prad_v.is_none())
                    ? rizer::PropertyTable::constant(0.0)
                    : makeTable(prad_T.cast<py::array_t<double>>(),
                                prad_v.cast<py::array_t<double>>());
            // Set up the solver options from the provided parameters.
            rizer::ColumnOptions opts;
            opts.n_points = n_points;
            opts.T_center_guess = T_center_guess;
            opts.rho_cp = rho_cp;
            opts.loglevel = loglevel;
            opts.refine_grid = refine_grid;
            opts.refine_ratio = refine_ratio;
            opts.refine_slope = refine_slope;
            opts.refine_curve = refine_curve;
            opts.refine_prune = refine_prune;
 
            // If the user provided an initial seed profile,
            // convert the numpy arrays into std::vectors.
            if (!init_r.is_none() && !init_T.is_none()) {
                auto ir = init_r.cast<py::array_t<double>>();
                auto it = init_T.cast<py::array_t<double>>();
                auto ira = ir.unchecked<1>();
                auto ita = it.unchecked<1>();
                opts.init_r.resize(static_cast<std::size_t>(ira.shape(0)));
                opts.init_T.resize(static_cast<std::size_t>(ita.shape(0)));
                for (py::ssize_t i = 0; i < ira.shape(0); i++) {
                    opts.init_r[static_cast<std::size_t>(i)] = ira(i);
                }
                for (py::ssize_t i = 0; i < ita.shape(0); i++) {
                    opts.init_T[static_cast<std::size_t>(i)] = ita(i);
                }
            }
 
            // Call the C++ solver function, releasing the GIL since it is pure C++.
            // The result is a ColumnResult struct containing the solution.
            rizer::ColumnResult res;
            {
                py::gil_scoped_release release;  // solve is pure C++
                res = rizer::solveColumn(R, electric_field, T_wall,
                                         sigma, kappa, prad, opts);
            }
 
            // Convert the result into a Python dict of numpy arrays and scalars.
            py::dict out;
            out["r"] = py::array_t<double>(res.r.size(), res.r.data());
            out["T"] = py::array_t<double>(res.T.size(), res.T.data());
            out["sigma"] = py::array_t<double>(res.sigma.size(), res.sigma.data());
            out["kappa"] = py::array_t<double>(res.kappa.size(), res.kappa.data());
            out["current"] = res.current;
            out["electric_field"] = res.electric_field;
            out["n_points"] = res.n_points;
            return out;
        },
        py::arg("R"), py::arg("electric_field"), py::arg("T_wall"),
        py::arg("sigma_T"), py::arg("sigma_v"),
        py::arg("kappa_T"), py::arg("kappa_v"),
        py::arg("prad_T") = py::none(), py::arg("prad_v") = py::none(),
        py::arg("n_points") = 41, py::arg("T_center_guess") = 10000.0,
        py::arg("rho_cp") = 1.0e3, py::arg("loglevel") = 0,
        py::arg("refine_grid") = true, py::arg("refine_ratio") = 4.0,
        py::arg("refine_slope") = 0.05, py::arg("refine_curve") = 0.10,
        py::arg("refine_prune") = 0.0,
        py::arg("init_r") = py::none(), py::arg("init_T") = py::none(),
        "Solve the steady radial Elenbaas-Heller column.\n\n"
        "Assembles [Empty1D | ThermalPlasmaColumn1D | Empty1D] in a Cantera\n"
        "Sim1D, solves with Newton + pseudo-transient fallback, and optionally\n"
        "refines the radial grid adaptively.\n\n"
        "Parameters\n"
        "----------\n"
        "R             : outer wall radius [m]\n"
        "electric_field: axial E field [V/m]\n"
        "T_wall        : wall temperature [K]\n"
        "sigma_T/v     : tabulated sigma(T) [S/m] — 1-D arrays of equal length\n"
        "kappa_T/v     : tabulated kappa(T) [W/m/K]\n"
        "prad_T/v      : tabulated P_rad(T) [W/m^3], or None for no radiation\n"
        "n_points      : initial uniform grid count\n"
        "T_center_guess: parabolic-seed centerline temperature [K]\n"
        "rho_cp        : rho*cp [J/m^3/K] for pseudo-transient scaling\n"
        "loglevel      : Sim1D verbosity (0 = silent)\n"
        "refine_grid   : enable adaptive grid refinement\n"
        "refine_ratio/slope/curve/prune: Refiner parameters\n"
        "init_r/T      : optional seed profile T(r); overrides the parabola\n\n"
        "Returns\n"
        "-------\n"
        "dict with keys r, T, sigma, kappa (numpy arrays) and\n"
        "current [A], electric_field [V/m], n_points (scalars).");
 
    // Advance the transient, two-temperature, finite-rate reacting 1-D radial
    // plasma channel (PlasmaChannel1D) from t=0 to t = dt*n_steps. Lambda
    // parameters (plain comment, see the note at mean_cross_section above --
    // this lambda has no separate declaration for \param to attach to):
    //   mech            Cantera mechanism (YAML) path
    //   phase           Cantera mechanism (YAML) phase
    //   n_points        Radial grid points [-] (ignored if `grid` is given)
    //   R_max           Outer (wall) radius [m]
    //   rho             Fixed mass density [kg/m^3] (constant-volume model)
    //   Tg0             Initial/ambient gas temperature [K] (uniform default
    //                   if initTg_r/v not given)
    //   Te0             Initial/ambient electron temperature [K] (uniform
    //                   default if initTe_r/v not given)
    //   Y0              Initial + ambient mass fractions [-], length nsp
    //   sigma_T         Legacy sigma(Te) table: temperature grid [K], or None
    //   sigma_v         Legacy sigma(Te) table: conductivity values [S/m], or
    //                   None
    //   nu_E            Legacy constant electron-heavy elastic exchange
    //                   frequency [1/s]
    //   electric_field  Constant applied axial field E [V/m] (superseded by
    //                   Efield_t/v)
    //   kappa_T         Gas thermal conductivity table: temperature grid [K],
    //                   or None
    //   kappa_v         Gas thermal conductivity table: kappa(Tg) values
    //                   [W/m/K], or None
    //   kappa_e_T       Electron thermal conductivity table: temperature grid
    //                   [K], or None
    //   kappa_e_v       Electron thermal conductivity table: kappa_e(Te)
    //                   values [W/m/K], or None
    //   D_species       Constant scalar species mass diffusivity [m^2/s]
    //                   (0 = frozen transport)
    //   reacting        true: include finite-rate chemistry source terms;
    //                   false: frozen composition
    //   T_amb           Wall/ambient temperature [K] (Dirichlet BC)
    //   initTg_r        Optional initial Tg(r) profile: radii [m], or None
    //   initTg_v        Optional initial Tg(r) profile: temperatures [K], or
    //                   None
    //   initTe_r        Optional initial Te(r) profile: radii [m], or None
    //   initTe_v        Optional initial Te(r) profile: temperatures [K], or
    //                   None
    //   grid            Optional explicit radial mesh [m] (overrides
    //                   n_points/R_max), or None
    //   nu_m            Legacy constant electron momentum-transfer frequency
    //                   [1/s] (>0 activates sigma = n_e e^2/(m_e nu_m))
    //   initY_r         Optional initial composition profile: radii [m],
    //                   length nr, or None
    //   initY_v         Optional initial composition profile: mass fractions
    //                   [-], row-major shape [nr, nsp], or None
    //   Efield_t        Optional field-vs-time schedule: time grid [s]
    //                   (overrides electric_field), or None
    //   Efield_v        Optional field-vs-time schedule: field values [V/m],
    //                   or None
    //   nu_m_T          Optional Te-dependent nu_m(Te) table: temperature
    //                   grid [K], or None (empty -> scalar nu_m)
    //   nu_m_v          Optional Te-dependent nu_m(Te) table: frequency
    //                   values [1/s], or None
    //   xsec_energy_J   Composition-resolved collision model: per-species
    //                   cross-section energy grids [J] (list of 1-D arrays
    //                   or None per species), or None to disable
    //   xsec_sigma_m2   Composition-resolved collision model: per-species
    //                   cross-section values [m^2] (list of 1-D arrays or
    //                   None per species), or None to disable
    //   radius_m        Composition-resolved collision model: per-species
    //                   hard-sphere effective radius [m] (NaN if unused), or
    //                   None to disable
    //   ion_Z           Composition-resolved collision model: per-species
    //                   ion charge number [-] (0 for neutrals), or None to
    //                   disable
    //   Te_min          Collision-model Te grid lower bound [K] (used only
    //                   if the collision model is enabled)
    //   Te_max          Collision-model Te grid upper bound [K] (used only
    //                   if the collision model is enabled)
    //   Te_n            Collision-model Te grid point count [-] (log-spaced;
    //                   used only if the collision model is enabled)
    //   spitzer         true: apply the Spitzer e-e correction to the
    //                   collision model
    //   species_h0k     Per-species standard enthalpy of formation at 0 K
    //                   [J/kmol], length nsp, or None; required when
    //                   `reacting` and the mechanism has an electron-impact
    //                   reaction (see rizer.kinetics.electron_reactions)
    //   dt              Output cadence / initial step size [s]
    //   n_steps         Number of output steps [-] (t_total = dt*n_steps)
    //   record_every    Record a frame every N output steps [-]
    //   loglevel        Diagnostic verbosity (0 = silent)
    //   integrator      "newton" (adaptive Backward-Euler) or "bdf" (stiff
    //                   CVODE)
    m.def("solve_channel_transient",
        [](const std::string& mech, const std::string& phase, std::size_t n_points,
           double R_max, double rho, double Tg0, double Te0, py::array_t<double> Y0,
           py::object sigma_T, py::object sigma_v, double nu_E, double electric_field,
           py::object kappa_T, py::object kappa_v, py::object kappa_e_T,
           py::object kappa_e_v, double D_species, bool reacting, double T_amb,
           py::object initTg_r, py::object initTg_v, py::object initTe_r,
           py::object initTe_v, py::object grid, double nu_m, py::object initY_r,
           py::object initY_v, py::object Efield_t, py::object Efield_v,
           py::object nu_m_T, py::object nu_m_v, py::object xsec_energy_J,
           py::object xsec_sigma_m2, py::object radius_m, py::object ion_Z,
           double Te_min, double Te_max, int Te_n, bool spitzer,
           py::object species_h0k, double dt,
           std::size_t n_steps, std::size_t record_every, int loglevel,
           const std::string& integrator)
        {
            // Every optional property/profile table below (sigma, kappa, kappa_e,
            // nu_m, grid, Et) is passed from Python as `None` when disabled;
            // ChannelOptions/PlasmaChannelConfig instead use empty vectors as
            // their "disabled" sentinel (see makeTable() in PlasmaChannel1D.cpp),
            // so this converts the Python convention to the C++ one in one place.
            auto optVec = [](py::object obj) -> std::vector<double> {
                if (obj.is_none()) return {};
                return toVec(obj.cast<py::array_t<double>>());
            };
 
            rizer::ChannelOptions options;
            auto& cfg = options.cfg;
            cfg.mech = mech;
            cfg.phase = phase;
            cfg.n_points = n_points;
            cfg.R_max = R_max;
            cfg.rho = rho;
            cfg.electric_field = electric_field;
            cfg.nu_E = nu_E;
            cfg.D_species = D_species;
            cfg.reacting = reacting;
            cfg.T_amb = T_amb;
            cfg.Y0 = toVec(Y0);
            cfg.species_h0k = optVec(species_h0k);
            cfg.sigma_v = optVec(sigma_v);
            cfg.sigma_T = optVec(sigma_T);
            cfg.kappa_v = optVec(kappa_v);
            cfg.kappa_T = optVec(kappa_T);
            cfg.kappa_e_v = optVec(kappa_e_v);
            cfg.kappa_e_T = optVec(kappa_e_T);
            // Initial profiles: explicit Tg(r)/Te(r) if given, else uniform Tg0/Te0.
            cfg.initTg_v = initTg_v.is_none() ? std::vector<double>{Tg0, Tg0}
                                            : optVec(initTg_v);
            cfg.initTg_r = initTg_r.is_none() ? std::vector<double>{0.0, R_max}
                                            : optVec(initTg_r);
            cfg.initTe_v = initTe_v.is_none() ? std::vector<double>{Te0, Te0}
                                            : optVec(initTe_v);
            cfg.initTe_r = initTe_r.is_none() ? std::vector<double>{0.0, R_max}
                                            : optVec(initTe_r);
            cfg.nu_m = nu_m;
            // Te-dependent nu_m(Te) table (empty -> scalar nu_m).
            cfg.nu_m_T = optVec(nu_m_T);
            cfg.nu_m_v = optVec(nu_m_v);
            // Explicit radial mesh (empty -> uniform).
            cfg.grid = optVec(grid);
            // Field-vs-time table (empty -> constant).
            cfg.Et_t = optVec(Efield_t);
            cfg.Et_v = optVec(Efield_v);
            // Initial composition profile: initY_r (nr) + initY_v (2D [nr, nsp]).
            if (!initY_r.is_none() && !initY_v.is_none()) {
                cfg.initY_r = toVec(initY_r.cast<py::array_t<double>>());
                auto initYArray = initY_v.cast<py::array_t<double>>();
                auto initYView = initYArray.unchecked<2>();
                cfg.initY_v.resize(
                    static_cast<std::size_t>(initYView.shape(0) * initYView.shape(1)));
                for (py::ssize_t i = 0; i < initYView.shape(0); i++) {
                    for (py::ssize_t k = 0; k < initYView.shape(1); k++) {
                        cfg.initY_v[
                            static_cast<std::size_t>(i * initYView.shape(1) + k)] =
                            initYView(i, k);
                    }
                }
            }
            // Composition-resolved collision model (empty -> legacy sigma/nu_m/nu_E
            // closure above). xsec_energy_J/xsec_sigma_m2/radius_m/ion_Z mirror the
            // Reactor0D binding exactly, one entry per species: all four must be
            // given together (or all omitted) since a per-species spec needs its
            // radius/Z slot even when it carries no tabulated cross-section.
            const bool any_collision_arg = !xsec_energy_J.is_none()
                || !xsec_sigma_m2.is_none() || !radius_m.is_none() || !ion_Z.is_none();
            if (any_collision_arg) {
                if (xsec_energy_J.is_none() || xsec_sigma_m2.is_none()
                    || radius_m.is_none() || ion_Z.is_none()) {
                    throw std::invalid_argument(
                        "solve_channel_transient: xsec_energy_J, xsec_sigma_m2, "
                        "radius_m, and ion_Z must all be given together (or all "
                        "omitted) to enable the composition-resolved collision "
                        "model.");
                }
                py::list e_list = xsec_energy_J.cast<py::list>();
                py::list s_list = xsec_sigma_m2.cast<py::list>();
                auto radius_arr = radius_m.cast<py::array_t<double>>();
                auto ion_z_arr = ion_Z.cast<py::array_t<int>>();
                const std::size_t nsp_specs = static_cast<std::size_t>(e_list.size());
                cfg.specs.resize(nsp_specs);
                auto r = radius_arr.unchecked<1>();
                auto z = ion_z_arr.unchecked<1>();
                for (std::size_t k = 0; k < nsp_specs; k++) {
                    cfg.specs[k].radius = r(static_cast<py::ssize_t>(k));
                    cfg.specs[k].Z = z(static_cast<py::ssize_t>(k));
                    py::object eo = e_list[k];
                    py::object so = s_list[k];
                    if (!eo.is_none() && !so.is_none()) {
                        cfg.specs[k].energy_J = toVec(eo.cast<py::array_t<double>>());
                        cfg.specs[k].sigma_m2 = toVec(so.cast<py::array_t<double>>());
                    }
                }
                cfg.Te_min = Te_min;
                cfg.Te_max = Te_max;
                cfg.Te_n = static_cast<std::size_t>(Te_n);
                cfg.spitzer = spitzer;
            }
 
            options.dt = dt;
            options.n_steps = n_steps;
            options.record_every = record_every;
            options.loglevel = loglevel;
            options.integrator = integrator;
 
            // A transient BDF/Newton solve over many steps can run for seconds
            // to minutes; release the GIL for its duration so other Python
            // threads (and, e.g., a Ctrl-C from the interpreter) aren't blocked.
            rizer::ChannelHistory history;
            {
                py::gil_scoped_release release;
                history = rizer::solveChannelTransient(options);
            }
 
            py::dict out;
            out["t"] = py::array_t<double>(history.t.size(), history.t.data());
            out["r"] = py::array_t<double>(history.r.size(), history.r.data());
            // Shaped [nframes, npts] and [nframes, npts, nsp].
            out["Tg"] = py::array_t<double>(
                {history.nframes, history.npts}, history.Tg.data());
            out["Te"] = py::array_t<double>(
                {history.nframes, history.npts}, history.Te.data());
            out["ne"] = py::array_t<double>(
                {history.nframes, history.npts}, history.ne.data());
            out["Y"] = py::array_t<double>(
                {history.nframes, history.npts, history.nsp}, history.Y.data());
            out["species"] = history.species;
            out["nframes"] = history.nframes;
            out["npts"] = history.npts;
            out["nsp"] = history.nsp;
            return out;
        },
        py::arg("mech"), py::arg("phase"), py::arg("n_points"),
        py::arg("R_max"), py::arg("rho"), py::arg("Tg0"), py::arg("Te0"),
        py::arg("Y0"), py::arg("sigma_T") = py::none(),
        py::arg("sigma_v") = py::none(), py::arg("nu_E") = 0.0,
        py::arg("electric_field") = 0.0,
        py::arg("kappa_T") = py::none(), py::arg("kappa_v") = py::none(),
        py::arg("kappa_e_T") = py::none(), py::arg("kappa_e_v") = py::none(),
        py::arg("D_species") = 0.0, py::arg("reacting") = true,
        py::arg("T_amb") = 300.0,
        py::arg("initTg_r") = py::none(), py::arg("initTg_v") = py::none(),
        py::arg("initTe_r") = py::none(), py::arg("initTe_v") = py::none(),
        py::arg("grid") = py::none(), py::arg("nu_m") = 0.0,
        py::arg("initY_r") = py::none(), py::arg("initY_v") = py::none(),
        py::arg("Efield_t") = py::none(), py::arg("Efield_v") = py::none(),
        py::arg("nu_m_T") = py::none(), py::arg("nu_m_v") = py::none(),
        py::arg("xsec_energy_J") = py::none(), py::arg("xsec_sigma_m2") = py::none(),
        py::arg("radius_m") = py::none(), py::arg("ion_Z") = py::none(),
        py::arg("Te_min") = 300.0, py::arg("Te_max") = 1.0e5,
        py::arg("Te_n") = 400, py::arg("spitzer") = true,
        py::arg("species_h0k") = py::none(),
        py::arg("dt") = 1.0e-9,
        py::arg("n_steps") = 1000, py::arg("record_every") = 1,
        py::arg("loglevel") = 0, py::arg("integrator") = "newton",
        "Advance the transient, two-temperature, finite-rate reacting 1-D radial\n"
        "plasma channel (PlasmaChannel1D) from t=0 to t=dt*n_steps.\n\n"
        "Parameters\n"
        "----------\n"
        "mech, phase   : Cantera mechanism (YAML) path and phase name\n"
        "n_points      : radial grid points [-] (ignored if `grid` is given)\n"
        "R_max         : outer (wall) radius [m]\n"
        "rho           : fixed mass density [kg/m^3] (constant-volume model)\n"
        "Tg0, Te0      : initial gas/electron temperature [K] (uniform default\n"
        "                if initTg_r/v or initTe_r/v are not given)\n"
        "Y0            : initial + ambient mass fractions [-], length nsp\n"
        "sigma_T/v, nu_E, nu_m, nu_m_T/v : legacy Joule/exchange closure -- sigma\n"
        "                from a sigma(Te) table, or n_e e^2/(m_e nu_m) from the\n"
        "                local n_e if nu_m>0 or a nu_m(Te) table is given, plus a\n"
        "                constant elastic frequency nu_E. Mutually exclusive with\n"
        "                xsec_energy_J/xsec_sigma_m2/radius_m/ion_Z below -- an\n"
        "                error is raised if both are supplied.\n"
        "electric_field: constant applied field E [V/m] (superseded by Efield_t/v)\n"
        "kappa_T/v     : tabulated gas thermal conductivity kappa(Tg) [K],[W/m/K]\n"
        "kappa_e_T/v   : tabulated electron thermal conductivity kappa_e(Te) "
        "[K],[W/m/K]\n"
        "D_species     : constant scalar species mass diffusivity [m^2/s] "
        "(0 = frozen)\n"
        "reacting      : include finite-rate chemistry source terms\n"
        "T_amb         : wall/ambient temperature [K] (Dirichlet BC)\n"
        "initTg_r/v    : initial radial profile Tg(r) [m],[K]\n"
        "initTe_r/v    : initial radial profile Te(r) [m],[K]\n"
        "grid          : explicit radial mesh [m] (overrides n_points/R_max)\n"
        "initY_r/v     : initial radial composition profile: radii [m] (length nr)\n"
        "                and row-major mass fractions [-] (shape [nr, nsp])\n"
        "Efield_t/v    : field-vs-time table [s],[V/m] (overrides electric_field)\n"
        "xsec_energy_J, xsec_sigma_m2, radius_m, ion_Z :\n"
        "                composition-resolved collision model (preferred over the\n"
        "                legacy closure above): per-species lists/arrays, in\n"
        "                species-index order (length nsp), exactly matching the\n"
        "                Reactor0D binding. xsec_energy_J/xsec_sigma_m2 are lists\n"
        "                of 1D arrays (or None per species) giving a tabulated\n"
        "                cross-section Q(E) [J],[m^2]; radius_m is a hard-sphere\n"
        "                effective radius [m] (NaN if not used); ion_Z is the ion\n"
        "                charge number (0 for neutrals). All four must be given\n"
        "                together (or all omitted).\n"
        "Te_min/Te_max/Te_n : collision-model Te grid [K],[K],[-] (log-spaced)\n"
        "spitzer       : apply the Spitzer e-e correction to the collision model\n"
        "species_h0k   : per-species standard enthalpy of formation at 0 K "
        "[J/kmol],\n"
        "                length nsp, or None; required when `reacting` and the\n"
        "                mechanism has an electron-impact reaction (see\n"
        "                rizer.kinetics.electron_reactions)\n"
        "dt            : output cadence / initial step [s]\n"
        "n_steps       : number of output steps [-] (dt*n_steps = total time [s])\n"
        "record_every  : record a frame every N output steps [-]\n"
        "loglevel      : diagnostic verbosity (0 = silent)\n"
        "integrator    : \"newton\" (adaptive Backward-Euler) or \"bdf\" "
        "(stiff CVODE)\n\n"
        "Returns\n"
        "-------\n"
        "dict with keys:\n"
        "  t       : recorded times [s], shape [nframes]\n"
        "  r       : radial grid [m], shape [npts]\n"
        "  Tg      : gas temperature [K], shape [nframes, npts]\n"
        "  Te      : electron temperature [K], shape [nframes, npts]\n"
        "  ne      : electron number density [1/m^3], shape [nframes, npts]\n"
        "  Y       : species mass fractions [-], shape [nframes, npts, nsp]\n"
        "  species : species names, length nsp (order matches the Y axis)\n"
        "  nframes, npts, nsp : array-shape scalars [-]");
}