bindings.cpp

Go to the documentation of this file.

// pybind11 bindings for the rizer Cantera plasma extension (_plasma1d).
//
// The extension holds rizer's C++ plasma models and this file is the only
// place they meet Python. It covers three unrelated problems, which is why the
// exposed API below reads as three groups rather than one:
//
//   - a steady LTE arc column (Elenbaas-Heller), solved on Cantera's Sim1D
//  with its damped-Newton driver;
//   - a transient two-temperature NRP plasma channel, 0D and 1D radial, with
//  finite-rate chemistry, Joule heating and electron-heavy exchange,
//  integrated by CVODE/BDF;
//   - the explicit hydrodynamic half-step that wraps that channel under the
//  Strang split of NUMERICS.md section 4.2, plus the thermodynamic
//  accessors it needs.
//
// 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 (abscissa, value) numpy arrays and are
//  converted to PropertyTable objects inside the binding lambda.
//   - Results leave as plain Python dicts of numpy arrays and Python scalars.
//
// Exposed API:
//   Library
//  cantera_version()           -> str   version of the linked libcantera
//  register_plasma_rates()            register rizer's native reaction rates
//  mean_cross_section(...)     -> array Maxwellian-averaged Qbar(Te)
//   Steady LTE arc column
//  solve_column(...)           -> dict  converged radial column profile
//   Transient NRP channel
//  Reactor0D                          persistent 0D two-temperature reactor
//  solve_1D_unsteady_discharge(...) -> dict  1D radial channel history
//   Hydrodynamic sub-step and its thermodynamics
//  solve_hydro_step_only(...)  -> dict  one explicit Strang half-step
//  hydro_thermo_probe(...)     -> dict  ThermoAccessor at one state
 
#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/PlasmaReactor0D.h"
#include "models/steady_lte_arc/PlasmaColumnSolver.h"
#include "kinetics/PlasmaRates.h"
#include "models/nrp/PlasmaChannelSolver.h"
#include "models/nrp/BoundaryConditions.h"
#include "models/nrp/HydroStep.h"
#include "models/nrp/PlasmaChannel1D.h"
#include "models/nrp/RadialDiffusion.h"
#include "thermo/CanteraThermo.h"
#include "thermo/GammaLawThermo.h"
#include "thermo/PrescribedPressureThermo.h"
#include "thermo/ThermoAccessor.h"
 
#include <cmath>
#include <filesystem>
#include <map>
#include <memory>
#include <optional>
#include <sstream>
#include <stdexcept>
#include <string>
#include <vector>
 
namespace py = pybind11;
 
namespace
{
 
    // ------------------------------------------------------------------------------ //
    //   File-local helpers                                                           //
    // ------------------------------------------------------------------------------ //

    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());
    }

    std::vector<double> optVec(const py::object &obj)
    {
        if (obj.is_none())
            return {};
        return toVec(obj.cast<py::array_t<double>>());
    }

    void requireLength(std::size_t size, std::size_t expected, const char *name)
    {
        if (size != expected)
        {
            throw std::invalid_argument(
                std::string(name) + " has length " + std::to_string(size) +
                ", expected " + std::to_string(expected) + ".");
        }
    }

    void requireTeGrid(int Te_n)
    {
        if (Te_n < 2)
        {
            throw std::invalid_argument(
                "Te_n must be at least 2 to span [Te_min, Te_max], got " +
                std::to_string(Te_n) + ".");
        }
    }

    void requireMatrixShape(const py::array &a, std::size_t rows,
                            std::optional<std::size_t> cols, const char *name)
    {
        const bool shape_ok =
            a.ndim() == 2 && static_cast<std::size_t>(a.shape(0)) == rows &&
            (!cols || static_cast<std::size_t>(a.shape(1)) == *cols);
        if (!shape_ok)
        {
            std::string got = "(";
            for (py::ssize_t d = 0; d < a.ndim(); d++)
            {
                got += (d ? ", " : "") + std::to_string(a.shape(d));
            }
            got += a.ndim() == 1 ? ",)" : ")";
            throw std::invalid_argument(
                std::string(name) + " has shape " + got + ", expected (" +
                std::to_string(rows) + ", " +
                (cols ? std::to_string(*cols) : std::string("n_species")) + ").");
        }
    }

    std::vector<int> toChargeNumbers(const py::array_t<double> &a)
    {
        const std::vector<double> values = toVec(a);
        std::vector<int> charges(values.size());
        for (std::size_t k = 0; k < values.size(); k++)
        {
            if (!std::isfinite(values[k]) || values[k] != std::floor(values[k]))
            {
                throw std::invalid_argument(
                    "ion_Z[" + std::to_string(k) + "] is " + std::to_string(values[k]) +
                    "; a charge number must be a whole number.");
            }
            charges[k] = static_cast<int>(values[k]);
        }
        return charges;
    }

    std::vector<rizer::CollisionModel::SpeciesSpec>
    collisionSpecs(const py::list &xsec_energy_J, const py::list &xsec_sigma_m2,
                   const py::array_t<double> &radius_m, const py::array_t<double> &ion_Z)
    {
        const std::size_t n_species = static_cast<std::size_t>(xsec_energy_J.size());
        requireLength(static_cast<std::size_t>(xsec_sigma_m2.size()), n_species,
                      "xsec_sigma_m2");
        requireLength(static_cast<std::size_t>(radius_m.size()), n_species, "radius_m");
        requireLength(static_cast<std::size_t>(ion_Z.size()), n_species, "ion_Z");
        const std::vector<double> radius = toVec(radius_m);
        const std::vector<int> charges = toChargeNumbers(ion_Z);
 
        std::vector<rizer::CollisionModel::SpeciesSpec> specs(n_species);
        for (std::size_t k = 0; k < n_species; k++)
        {
            specs[k].radius = radius[k];
            specs[k].Z = charges[k];
            const py::object energy = xsec_energy_J[k];
            const py::object sigma = xsec_sigma_m2[k];
            if (energy.is_none() != sigma.is_none())
            {
                throw std::invalid_argument(
                    "species " + std::to_string(k) + " has " +
                    (energy.is_none() ? "xsec_sigma_m2 but no xsec_energy_J"
                                      : "xsec_energy_J but no xsec_sigma_m2") +
                    "; a tabulated cross section needs both, or neither for a "
                    "hard sphere or no model.");
            }
            if (!energy.is_none())
            {
                specs[k].energy_J = toVec(energy.cast<py::array_t<double>>());
                specs[k].sigma_m2 = toVec(sigma.cast<py::array_t<double>>());
            }
        }
        return specs;
    }

    rizer::PropertyTable makeTable(const py::array_t<double> &T,
                                   const py::array_t<double> &values,
                                   const char *name = "property table")
    {
        std::vector<double> abscissa = toVec(T);
        std::vector<double> ordinate = toVec(values);
        if (abscissa.size() != ordinate.size())
        {
            throw std::invalid_argument(
                std::string(name) +
                ": the abscissa has " +
                std::to_string(abscissa.size()) +
                " point(s) against " +
                std::to_string(ordinate.size()) +
                " value(s); they tabulate one curve and must be the same length.");
        }
        for (std::size_t i = 1; i < abscissa.size(); i++)
        {
            if (!(abscissa[i] > abscissa[i - 1]))
            {
                throw std::invalid_argument(
                    std::string(name) +
                    ": the abscissa must be strictly increasing, but entry " +
                    std::to_string(i) +
                    " (" + std::to_string(abscissa[i]) +
                    ") does not exceed entry " +
                    std::to_string(i - 1) +
                    " (" + std::to_string(abscissa[i - 1]) + ").");
            }
        }
        return rizer::PropertyTable(std::move(abscissa), std::move(ordinate));
    }
 
} // namespace
 
namespace
{
 
    // ------------------------------------------------------------------------------ //
    //   Cached Cantera thermo accessor                                               //
    // ------------------------------------------------------------------------------ //

    rizer::CanteraThermo &cachedCanteraThermo(const std::string &mech,
                                              const std::string &phase)
    {
        struct Entry
        {
            std::filesystem::file_time_type modified;
            std::unique_ptr<rizer::CanteraThermo> thermo;
        };
        static std::map<std::pair<std::string, std::string>, Entry> cache;
        const std::string path = Cantera::findInputFile(mech);
        const auto modified = std::filesystem::last_write_time(path);
        Entry &entry = cache[std::make_pair(path, phase)];
        if (!entry.thermo || entry.modified != modified)
        {
            entry.thermo = std::make_unique<rizer::CanteraThermo>(path, phase);
            entry.modified = modified;
        }
        return *entry.thermo;
    }
 
} // namespace
 
// ---------------------------------------------------------------------------------- //
//   Module definition                                                                //
// ---------------------------------------------------------------------------------- //
 
PYBIND11_MODULE(_plasma1d, m)
{
    m.doc() =
        "rizer Cantera plasma extension: a steady LTE arc column, a "
        "transient two-temperature NRP plasma channel (0D and 1D radial) "
        "with its explicit hydrodynamic Strang half-step, and the "
        "thermodynamic accessors they share. Flat numerics only across "
        "the boundary: no Cantera C++ object is ever exposed, so this "
        "module and the stock `cantera` package coexist in one process.";
 
    m.def("cantera_version", []()
          { return Cantera::version(); }, "Version string of the linked libcantera.");
 
    // Whether this build compiled the per-node reaction loop with OpenMP.
    // Exported because its ABSENCE is otherwise invisible: `n_threads` is
    // still accepted and still returns the right answer, one thread at a
    // time, so a CMake change that stops attaching OpenMP to plasma1d_core
    // costs the parallel speed-up with nothing failing. This flag is read
    // from bindings.cpp's own translation unit, which inherits the compile
    // option through plasma1d_core, so it tracks the core's setting.
#ifdef _OPENMP
    m.attr("has_openmp") = true;
#else
    m.attr("has_openmp") = false;
#endif
 
    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.");
 
    m.def("mean_cross_section",
          // 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, held at the
          // nearest edge value outside the tabulated range rather than zeroed.
          [](
              // Parameters
              // ----------
              py::array_t<double> energy_J,
              //  Tabulated collision energy grid [J], strictly increasing.
              py::array_t<double> sigma_m2,
              //  Tabulated cross section Q(energy_J) [m^2], same length as energy_J.
              py::array_t<double> Te,
              //  Electron temperatures to evaluate Qbar(Te) at [K].
              double x_max,
              //  Upper limit of the reduced-energy (E/kTe) integration [-].
              //  [Optional] By default, 20.0
              int N
              //  Number of quadrature points for the integration [-].
              //  [Optional] By default, 1000
              //
              // Returns
              // -------
              // NumPy array
              //  Maxwellian momentum-transfer kernel Qbar(Te) [m^2], one entry per Te.
          )
          {
            std::vector<double> e = toVec(energy_J);
            std::vector<double> s = toVec(sigma_m2);
            std::vector<double> t = toVec(Te);
            // Checked here rather than left to the table: without it a
            // mismatched pair would interpolate to 0.0 at every energy and
            // return a correctly shaped array of zeros instead of raising.
            if (e.size() != s.size()) {
                throw std::invalid_argument(
                    "mean_cross_section: energy_J has "
                    + std::to_string(e.size()) + " entries but sigma_m2 has "
                    + std::to_string(s.size()) + "; they tabulate the same "
                    "cross section and must be the same length.");
            }
            if (e.size() < 2) {
                throw std::invalid_argument(
                    "mean_cross_section: energy_J/sigma_m2 need at least 2 "
                    "points to interpolate, got " + std::to_string(e.size())
                    + ".");
            }
            std::vector<double> q =
                rizer::CollisionModel::meanCrossSectionTable(e, s, t, x_max, N);
            return py::array_t<double>(q.size(), q.data()); },
          // --- required
          py::arg("energy_J"), py::arg("sigma_m2"), py::arg("Te"),
          // --- optional
          py::arg("x_max") = 20.0, py::arg("N") = 1000);
 
    // 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::PlasmaReactor0D>(m, "Reactor0D")
        .def(py::init(
                 // Construct a Reactor0D from a mechanism and a collision model.
                 //
                 // The per-species cross-section model is passed once here as ragged
                 // lists and baked into a Qbar(Te) table, rather than integrated live
                 // on every right-hand-side evaluation.
                 [](
                     // Parameters
                     // ----------
                     //
                     // Mechanism
                     const std::string &mech,
                     //  Cantera mechanism (YAML) path
                     const std::string &phase,
                     //  Phase name within the mechanism
                     bool reacting,
                     //  Include finite-rate chemistry source terms [-]; false freezes
                     //  the composition
                     //
                     // Collision model, per species
                     py::list xsec_energy_J,
                     //  Tabulated cross-section energy grids [J] (list of 1-D arrays,
                     //  or None per species), length nsp
                     py::list xsec_sigma_m2,
                     //  Tabulated cross sections at xsec_energy_J [m^2] (list of 1-D
                     //  arrays, or None per species), length nsp
                     py::array_t<double> radius_m,
                     //  Hard-sphere effective radii [m], NaN where unused, length nsp
                     py::array_t<double> ion_Z,
                     //  Ion charge numbers [-], 0 for neutrals, length nsp
                     double Te_min,
                     //  Lower bound of the log-spaced Qbar(Te) grid [K]
                     double Te_max,
                     //  Upper bound of the log-spaced Qbar(Te) grid [K]
                     int Te_n,
                     //  Points in the Qbar(Te) grid [-]
                     //
                     // Geometry and thermochemistry
                     double gap,
                     //  Inter-electrode gap [m], used by the polytropic volume equation
                     //  [Optional] By default, 0.0
                     py::array_t<double> species_h0k
                     //  Standard enthalpies of formation at 0 K [J/kmol], length nsp;
                     //  required when `reacting` and the mechanism carries an
                     //  electron-impact reaction (see rizer.kinetics.electron_reactions)
                     //  [Optional] By default, empty
                 )
                 {
                     rizer::PlasmaReactor0DConfig cfg;
                     cfg.mech = mech;
                     cfg.phase = phase;
                     cfg.reacting = reacting;
                     cfg.Te_min = Te_min;
                     cfg.Te_max = Te_max;
                     requireTeGrid(Te_n);
                     cfg.Te_n = static_cast<std::size_t>(Te_n);
                     cfg.gap = gap;
                     cfg.species_h0k = toVec(species_h0k);
                     cfg.specs = collisionSpecs(xsec_energy_J, xsec_sigma_m2, radius_m,
                                                ion_Z);
                     return std::make_unique<rizer::PlasmaReactor0D>(cfg);
                 }),
             // --- required: mechanism
             py::arg("mech"),
             py::arg("phase"),
             // --- optional: mechanism
             py::arg("reacting") = true,
             // --- required: collision model
             py::arg("xsec_energy_J"),
             py::arg("xsec_sigma_m2"),
             py::arg("radius_m"),
             py::arg("ion_Z"),
             // --- optional: collision model
             py::arg("Te_min") = 200.0,
             py::arg("Te_max") = 1.0e5,
             py::arg("Te_n") = 400,
             // --- optional: geometry and thermochemistry
             py::arg("gap") = 0.0,
             py::arg("species_h0k") = py::array_t<double>(0))
 
        .def("electrical_conductivity",
             // Electrical conductivity at one state.
             [](
                 // Parameters
                 // ----------
                 rizer::PlasmaReactor0D &self,
                 //  The bound Reactor0D instance
                 double Tg,
                 //  Heavy-species temperature [K]
                 double Te,
                 //  Electron temperature [K]
                 py::array_t<double> Y,
                 //  Species mass fractions [-], length nsp
                 double rho
                 //  Mass density [kg/m^3]
                 //
                 // Returns
                 // -------
                 // float
                 //  Electrical conductivity [S/m]
             )
             {
                auto y = toVec(Y);
                requireLength(y.size(), self.nSpecies(), "Y");
                return self.electricalConductivity(Tg, Te, y.data(), rho); },
             // --- required
             py::arg("Tg"), py::arg("Te"), py::arg("Y"), py::arg("rho"))
 
        .def("electron_density",
             // Electron number density [1/m^3]
             [](
                 // Parameters
                 // ----------
                 rizer::PlasmaReactor0D &self,
                 //  The bound Reactor0D instance
                 double Tg,
                 //  Gas (heavy-species) temperature [K]
                 double Te,
                 //  Electron temperature [K]
                 py::array_t<double> Y,
                 //  Species mass fractions [-], length nsp
                 double rho
                 //  Mass density [kg/m^3]
                 //
                 // Returns
                 // -------
                 // float
                 //  Electron number density [1/m^3]
             )
             {
                auto y = toVec(Y);
                requireLength(y.size(), self.nSpecies(), "Y");
                return self.electronDensity(Tg, Te, y.data(), rho); },
             // --- required
             py::arg("Tg"), py::arg("Te"), py::arg("Y"), py::arg("rho"))
 
        .def("rhs",
             // 0D RHS dydt = [dTg, dTe, dV, dY...] for state [Tg, Te, V, Y...];
             [](
                 // Parameters
                 // ----------
                 rizer::PlasmaReactor0D &self,
                 //  The bound Reactor0D instance
                 double Tg,
                 //  Gas (heavy-species) temperature [K]
                 double Te,
                 //  Electron temperature [K]
                 double V,
                 //  Volume [m^3]
                 py::array_t<double> Y,
                 //  Species mass fractions [-], length nsp
                 double E,
                 //  Applied electric field [V/m]
                 double mass,
                 //  Fixed total mass [kg] (rho = mass/V)
                 double p_ext,
                 //  External pressure [Pa] (drives the polytropic volume equation)
                 double polytropic_index
                 //  Polytropic index [-] (inf = constant volume)
                 //
                 // Returns
                 // -------
                 // NumPy array
                 //  dydt = [dTg, dTe, dV, dY...] (length 3 + nsp)
 
             )
             {
                auto y = toVec(Y);
                requireLength(y.size(), self.nSpecies(), "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()); },
             // --- required
             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"))
 
        .def_property_readonly("species", [](rizer::PlasmaReactor0D &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_H/P_chemical_e/P_Joule/
        // nu_eH/nu_ee/nu_eI/nu_eH_mass_weighted attributes, natively).
        .def_property_readonly("P_elastic", [](rizer::PlasmaReactor0D &self)
                               { return self.diagnostics().P_elastic; })
        .def_property_readonly("P_inelastic", [](rizer::PlasmaReactor0D &self)
                               { return self.diagnostics().P_inelastic; })
        .def_property_readonly("P_chemical_H", [](rizer::PlasmaReactor0D &self)
                               { return self.diagnostics().P_chemical_H; })
        .def_property_readonly("P_chemical_e", [](rizer::PlasmaReactor0D &self)
                               { return self.diagnostics().P_chemical_e; })
        .def_property_readonly("P_Joule", [](rizer::PlasmaReactor0D &self)
                               { return self.diagnostics().P_Joule; })
        .def_property_readonly("nu_eH", [](rizer::PlasmaReactor0D &self)
                               { return self.diagnostics().nu_eH; })
        .def_property_readonly("nu_ee", [](rizer::PlasmaReactor0D &self)
                               { return self.diagnostics().nu_ee; })
        .def_property_readonly("nu_eI", [](rizer::PlasmaReactor0D &self)
                               { return self.diagnostics().nu_eI; })
        .def_property_readonly("nu_eH_mass_weighted", [](rizer::PlasmaReactor0D &self)
                               { return self.diagnostics().nu_eH_mass_weighted; });
 
    m.def("solve_column",
          // Solve the steady radial Elenbaas-Heller LTE arc column.
          //
          // Assembles [Empty1D | ThermalPlasmaColumn1D | Empty1D] in a Cantera
          // Sim1D, solves with its damped-Newton driver plus a pseudo-transient
          // fallback, and optionally refines the radial grid adaptively.
          [](
              // Parameters
              // ----------
              //
              // Geometry and operating point
              double R,
              //  Outer wall radius [m]
              double electric_field,
              //  Axial electric field [V/m]
              double T_wall,
              //  Wall temperature [K]
              //
              // Material properties
              py::array_t<double> sigma_T,
              //  Temperatures for sigma_v [K], strictly increasing
              py::array_t<double> sigma_v,
              //  Electrical conductivity at sigma_T [S/m]
              py::array_t<double> kappa_T,
              //  Temperatures for kappa_v [K], strictly increasing
              py::array_t<double> kappa_v,
              //  Thermal conductivity at kappa_T [W/m/K]
              py::object prad_T,
              //  Temperatures for prad_v [K], strictly increasing
              //  [Optional] By default, None (radiation off)
              py::object prad_v,
              //  Radiated power density at prad_T [W/m^3]
              //  [Optional] By default, None (radiation off)
              //
              // Initial guess
              std::size_t n_points,
              //  Initial uniform grid count [-]
              //  [Optional] By default, 41
              double T_center_guess,
              //  Parabolic-seed centerline temperature [K]
              //  [Optional] By default, 10000.0
              py::object init_r,
              //  Seed radii [m], strictly increasing over 0..R, overriding the
              //  parabolic seed
              //  [Optional] By default, None
              py::object init_T,
              //  Seed temperatures at init_r [K]
              //  [Optional] By default, None
              //
              // Solver controls
              double rho_cp,
              //  rho*cp for the pseudo-transient scaling [J/m^3/K]
              //  [Optional] By default, 1.0e3
              int loglevel,
              //  Sim1D verbosity [-]; 0 is silent
              //  [Optional] By default, 0
              //
              // Grid refinement
              bool refine_grid,
              //  Enable Cantera's adaptive grid refinement [-]
              //  [Optional] By default, true
              double refine_ratio,
              //  Refiner: max ratio of adjacent cell widths [-]
              //  [Optional] By default, 4.0
              double refine_slope,
              //  Refiner: max fractional slope change per cell [-]
              //  [Optional] By default, 0.05
              double refine_curve,
              //  Refiner: max fractional curvature per cell [-]
              //  [Optional] By default, 0.10
              double refine_prune
              //  Refiner: remove points below this threshold [-]
              //  [Optional] By default, 0.0
              //
              // Returns
              // -------
              // dict
              //  The converged profile: r [m], T [K], sigma [S/m], kappa
              //  [W/m/K], the total current [A], the electric field [V/m] and
              //  the final n_points.
          )
          {
            // Convert the input numpy arrays into PropertyTable objects.
            rizer::PropertyTable sigma =
                makeTable(sigma_T, sigma_v, "solve_column: sigma_T/sigma_v");
            rizer::PropertyTable kappa =
                makeTable(kappa_T, kappa_v, "solve_column: 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>>(),
                                "solve_column: prad_T/prad_v");
 
            // 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.
            // solveColumn refuses a pair it cannot use, so both go through even
            // when only one is given.
            opts.init_r = optVec(init_r);
            opts.init_T = optVec(init_T);
 
            // 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; },
          // --- required: geometry and operating point (outer radius, field, wall temperature)
          py::arg("R"), py::arg("electric_field"), py::arg("T_wall"),
          // --- required: plasma properties (electrical and thermal conductivity)
          py::arg("sigma_T"), py::arg("sigma_v"), py::arg("kappa_T"), py::arg("kappa_v"),
          // --- optional: plasma properties (radiation)
          py::arg("prad_T") = py::none(), py::arg("prad_v") = py::none(),
          // --- optional: initial guess (parabolic seed or explicit profile)
          py::arg("n_points") = 41, py::arg("T_center_guess") = 10000.0, py::arg("init_r") = py::none(), py::arg("init_T") = py::none(),
          // --- optional: solver controls (pseudo-transient scaling and verbosity)
          py::arg("rho_cp") = 1.0e3, py::arg("loglevel") = 0,
          // --- optional: grid refinement (Cantera's adaptive grid refinement)
          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);
 
    m.def("solve_1D_unsteady_discharge",
          // Integrate the transient two-temperature NRP plasma channel.
          //
          // One radial domain of `n_points` nodes carrying Tg, Te and the K
          // species mass fractions, advanced by CVODE/BDF over `n_steps` steps
          // of `dt`. With `hydrodynamics` on, each step is the Strang sequence
          // of NUMERICS.md section 4.2 -- half an explicit hydro step, the
          // reaction step, half an explicit hydro step -- and rho and u_r
          // become solved fields rather than fixed inputs.
          //
          // Setting n_points = 1 gives the 0D reactor as the single-node limit.
          [](
              // Parameters
              // ----------
              //
              // Mechanism
              const std::string &mech,
              //  Cantera mechanism (YAML) path
              const std::string &phase,
              //  Phase name within the mechanism
              //
              // Grid
              std::size_t n_points,
              //  Number of radial nodes [-]; 1 is the 0D limit
              double R_max,
              //  Outer radius of the domain [m]
              py::object grid,
              //  Explicit node radii [m], overriding the uniform grid built from
              //  n_points/R_max
              //  [Optional] By default, None (uniform)
              //
              // Initial state, uniform
              double rho,
              //  Mass density [kg/m^3], held fixed unless `hydrodynamics`
              double Tg0,
              //  Initial heavy-pool temperature [K]
              double Te0,
              //  Initial electron temperature [K]
              py::array_t<double> Y0,
              //  Initial mass fractions [-], length K, summing to 1
              //
              // Initial state, radial profiles (each a radius/value pair)
              py::object initTg_r,
              //  Radii for initTg_v [m]
              //  [Optional] By default, None (uniform Tg0)
              py::object initTg_v,
              //  Heavy-pool temperature at initTg_r [K]
              //  [Optional] By default, None
              py::object initTe_r,
              //  Radii for initTe_v [m]
              //  [Optional] By default, None (uniform Te0)
              py::object initTe_v,
              //  Electron temperature at initTe_r [K]
              //  [Optional] By default, None
              py::object initY_r,
              //  Radii for initY_v [m]
              //  [Optional] By default, None (uniform Y0)
              py::object initY_v,
              //  Mass fractions at initY_r [-], shape [n_radii, K]
              //  [Optional] By default, None
              py::object init_rho_r,
              //  Radii for init_rho_v [m]
              //  [Optional] By default, None (uniform rho)
              py::object init_rho_v,
              //  Mass density at init_rho_r [kg/m^3]
              //  [Optional] By default, None
              py::object init_ur_r,
              //  Radii for init_ur_v [m]
              //  [Optional] By default, None (gas at rest)
              py::object init_ur_v,
              //  Radial velocity at init_ur_r [m/s]; needs `hydrodynamics`
              //  [Optional] By default, None
              //
              // Chemistry
              bool reacting,
              //  Include finite-rate chemistry source terms [-]
              //  [Optional] By default, true
              py::object species_h0k,
              //  Per-species 0 K formation enthalpies [J/kmol], length K, for the
              //  inelastic threshold energies
              //  [Optional] By default, None (taken from the mechanism)
              //
              // Discharge
              double electric_field,
              //  Axial electric field [V/m], constant in time
              //  [Optional] By default, 0.0
              py::object Efield_t,
              //  Times for Efield_v [s], overriding the constant field
              //  [Optional] By default, None
              py::object Efield_v,
              //  Axial electric field at Efield_t [V/m]
              //  [Optional] By default, None
              double gap,
              //  Electrode gap [m], needed to turn `plasma_voltage` into a field
              //  [Optional] By default, 0.0
              py::object plasma_voltage,
              //  Callable V_p(t [s], R_p [Ohm]) -> volts across the gap,
              //  evaluated at both ends of every global step with R_p frozen
              //  over it; incompatible with n_threads > 1
              //  [Optional] By default, None
              //
              // Transport and collisions, bulk closures
              py::object sigma_T,
              //  Temperatures for sigma_v [K]
              //  [Optional] By default, None
              py::object sigma_v,
              //  Electrical conductivity at sigma_T [S/m]
              //  [Optional] By default, None
              double nu_E,
              //  Fixed electron-heavy energy-exchange frequency [1/s]
              //  [Optional] By default, 0.0
              double nu_m,
              //  Fixed electron momentum-transfer frequency [1/s], from which
              //  sigma follows by the Drude relation
              //  [Optional] By default, 0.0
              py::object nu_m_T,
              //  Temperatures for nu_m_v [K]
              //  [Optional] By default, None
              py::object nu_m_v,
              //  Momentum-transfer frequency at nu_m_T [1/s]
              //  [Optional] By default, None
              py::object kappa_H_T,
              //  Temperatures for kappa_H_v [K]
              //  [Optional] By default, None
              py::object kappa_H_v,
              //  Heavy-pool thermal conductivity at kappa_H_T [W/m/K]
              //  [Optional] By default, None (conduction off)
              py::object kappa_e_T,
              //  Temperatures for kappa_e_v [K]
              //  [Optional] By default, None
              py::object kappa_e_v,
              //  Electron thermal conductivity at kappa_e_T [W/m/K]
              //  [Optional] By default, None (conduction off)
              double D_species,
              //  Species diffusivity [m^2/s]; 0 disables species diffusion
              //  [Optional] By default, 0.0
              //
              // Transport and collisions, composition-resolved model
              py::object xsec_energy_J,
              //  Per-species cross-section energy grids [J] (list of 1-D arrays
              //  or None per species), length K
              //  [Optional] By default, None
              py::object xsec_sigma_m2,
              //  Per-species cross sections at xsec_energy_J [m^2], length K
              //  [Optional] By default, None
              py::object radius_m,
              //  Per-species hard-sphere effective radii [m], length K
              //  [Optional] By default, None
              py::object ion_Z,
              //  Per-species ion charge numbers [-], 0 for neutrals, length K
              //  [Optional] By default, None
              double Te_min,
              //  Lower edge of the precomputed Qbar(Te) table [K]
              //  [Optional] By default, 300.0
              double Te_max,
              //  Upper edge of the precomputed Qbar(Te) table [K]
              //  [Optional] By default, 1.0e5
              int Te_n,
              //  Points in the precomputed Qbar(Te) table [-]
              //  [Optional] By default, 400
              //
              // Ambient and boundary
              double T_amb,
              //  Ambient temperature at the wall [K]
              //  [Optional] By default, 300.0
              double p_amb,
              //  Ambient pressure [Pa]; required when `hydrodynamics`
              //  [Optional] By default, -1.0 (unset)
              py::object Y_amb,
              //  Ambient mass fractions [-], length K, summing to 1; required by
              //  outer_bc_mode="ambient_ghost"
              //  [Optional] By default, None
              const std::string &outer_bc_mode,
              //  Outer boundary treatment, "ambient_ghost" or "pressure_outlet"
              //  [Optional] By default, "ambient_ghost"
              //
              // Hydrodynamics
              bool hydrodynamics,
              //  Wrap each step in the explicit hydro Strang half-steps [-]
              //  [Optional] By default, false
              double cfl,
              //  Courant number for the hydro sub-cycling [-], in (0, 1]
              //  [Optional] By default, 0.9
              double hydro_substep,
              //  Pinned hydro sub-step [s], dividing dt/2 exactly; a test control
              //  for convergence studies, with no CFL check
              //  [Optional] By default, 0.0 (the CFL bound chooses)
              //
              // Time stepping
              double dt,
              //  Global (Strang) step [s]
              //  [Optional] By default, 1.0e-9
              std::size_t n_steps,
              //  Number of global steps [-]
              //  [Optional] By default, 1000
              std::size_t record_every,
              //  Record one frame every this many steps [-]
              //  [Optional] By default, 1
              //
              // Integrator controls
              double rtol,
              //  CVODE relative tolerance, shared by every component [-]
              //  [Optional] By default, 1.0e-7
              double atol_T,
              //  CVODE absolute tolerance on both temperatures [K]
              //  [Optional] By default, 1.0e-3
              double atol_Y,
              //  CVODE absolute tolerance on each mass fraction [-]
              //  [Optional] By default, 1.0e-9
              bool independent_nodes,
              //  Give each node its own CVODE instead of one banded system;
              //  legal only with every transport term off
              //  [Optional] By default, false
              std::size_t n_threads,
              //  OpenMP threads over the per-node loop [-]; needs
              //  `independent_nodes`
              //  [Optional] By default, 1
              int loglevel
              //  Solver verbosity [-]
              //  [Optional] By default, 0
              //
              // Returns
              // -------
              // dict
              //  Recorded history: t [s], r [m], the Tg/Te/ne/Y/rho/u_r/p fields
              //  shaped [nframes, npts] (Y as [nframes, npts, nsp]), the
              //  per-frame solver counters, and the species names.
          )
          {
            // Every optional property/profile table below (sigma, kappa,
            // kappa_e, nu_m, grid, Et) arrives as `None` when disabled;
            // ChannelOptions/PlasmaChannelConfig use empty vectors as that
            // sentinel instead, which optVec() converts.
            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.species_h0k = optVec(species_h0k);
            cfg.sigma_v = optVec(sigma_v);
            cfg.sigma_T = optVec(sigma_T);
            cfg.kappa_H_v = optVec(kappa_H_v);
            cfg.kappa_H_T = optVec(kappa_H_T);
            cfg.kappa_e_v = optVec(kappa_e_v);
            cfg.kappa_e_T = optVec(kappa_e_T);
            // Explicit radial mesh (empty -> uniform).
            cfg.grid = optVec(grid);
            // Initial profiles: the explicit radial table when given, else a
            // two-radius table carrying the uniform scalar. The domain takes
            // only tables, so the uniform case is built here rather than
            // doubling every field in PlasmaChannelConfig. The uniform table
            // spans the grid the domain will use: an explicit grid's own extent,
            // which need not be R_max.
            const double outer_radius = cfg.grid.size() >= 2 ? cfg.grid.back() : R_max;
            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, outer_radius}
                                            : 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, outer_radius}
                                            : optVec(initTe_r);
            if (initY_v.is_none()) {
                // Row-major [2 * K]: the same composition at both radii.
                const std::vector<double> Y0_v = toVec(Y0);
                cfg.initY_r = {0.0, outer_radius};
                cfg.initY_v = Y0_v;
                cfg.initY_v.insert(cfg.initY_v.end(), Y0_v.begin(), Y0_v.end());
            } else {
                // initY_v is [nr, nsp]; toVec's C-order copy is the row-major
                // layout PlasmaChannelConfig::initY_v expects.
                cfg.initY_r = optVec(initY_r);
                const auto initY_array = initY_v.cast<py::array_t<double>>();
                requireMatrixShape(initY_array, cfg.initY_r.size(),
                                   static_cast<std::size_t>(Y0.size()), "initY_v");
                cfg.initY_v = toVec(initY_array);
            }
            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);
            // Field-vs-time table (empty -> constant).
            options.Et_t = optVec(Efield_t);
            options.Et_v = optVec(Efield_v);
            cfg.init_rho_r = optVec(init_rho_r);
            cfg.init_rho_v = optVec(init_rho_v);
            cfg.hydrodynamics = hydrodynamics;
            cfg.init_ur_r = optVec(init_ur_r);
            cfg.init_ur_v = optVec(init_ur_v);
            cfg.outer_bc_mode = outer_bc_mode;
            cfg.p_amb = p_amb;
            cfg.cfl = cfl;
            cfg.hydro_substep = hydro_substep;
            // 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_1D_unsteady_discharge: 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.");
                }
                cfg.specs = collisionSpecs(xsec_energy_J.cast<py::list>(),
                                           xsec_sigma_m2.cast<py::list>(),
                                           radius_m.cast<py::array_t<double>>(),
                                           ion_Z.cast<py::array_t<double>>());
                cfg.Te_min = Te_min;
                cfg.Te_max = Te_max;
                requireTeGrid(Te_n);
                cfg.Te_n = static_cast<std::size_t>(Te_n);
            }
 
            options.dt = dt;
            options.n_steps = n_steps;
            options.record_every = record_every;
            options.loglevel = loglevel;
            options.rtol = rtol;
            options.atol_T = atol_T;
            options.atol_Y = atol_Y;
            options.independent_nodes = independent_nodes;
            options.n_threads = n_threads;
            if (Y_amb.is_none()) {
                throw std::invalid_argument(
                    "solve_1D_unsteady_discharge: Y_amb is required. It is the "
                    "composition OUTSIDE the domain -- the wall ghost and the "
                    "hydro inflow read it -- and it is deliberately not derived "
                    "from Y0, which is the composition INSIDE at t=0. Pass "
                    "Y_amb=Y0 only when the boundary genuinely holds the "
                    "initial composition.");
            }
            cfg.Y_amb = toVec(Y_amb.cast<py::array_t<double>>());
            options.gap = gap;
            if (!plasma_voltage.is_none()) {
                // The solve below releases the GIL, so every call back into
                // Python has to take it again. Captured by value: the solve
                // outlives this lambda's enclosing scope only through
                // ChannelOptions, which owns the std::function.
                py::function callback = plasma_voltage;
                options.plasma_voltage =
                    [callback](double t, double R_p) -> double {
                        py::gil_scoped_acquire gil;
                        return callback(t, R_p).cast<double>();
                    };
            }
 
            // 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["rho"] = py::array_t<double>(
                {history.nframes, history.npts}, history.rho.data());
            out["u_r"] = py::array_t<double>(
                {history.nframes, history.npts}, history.u_r.data());
            out["p"] = py::array_t<double>(
                {history.nframes, history.npts}, history.p.data());
            out["boundary_mass_out"] = history.boundary_mass_out;
            out["boundary_species_mass_out"] = py::array_t<double>(
                history.boundary_species_mass_out.size(),
                history.boundary_species_mass_out.data());
            out["boundary_total_energy_out"] = history.boundary_total_energy_out;
            out["boundary_electron_energy_out"] =
                history.boundary_electron_energy_out;
            out["pressure_work_electron"] = history.pressure_work_electron;
            out["boundary_conduction_energy_out"] =
                history.boundary_conduction_energy_out;
            out["boundary_diffusion_energy_out"] =
                history.boundary_diffusion_energy_out;
            // Per-frame CVODE counters.
            for (const auto& [key, values] : {
                     std::pair{"solver_steps", &history.solver_steps},
                     std::pair{"solver_rhs_evals", &history.solver_rhs_evals},
                     std::pair{"solver_jac_evals", &history.solver_jac_evals},
                     std::pair{"solver_err_test_fails",
                               &history.solver_err_test_fails},
                     std::pair{"solver_nonlin_conv_fails",
                               &history.solver_nonlin_conv_fails},
                     std::pair{"solver_last_order", &history.solver_last_order},
                     std::pair{"solver_lin_rhs_evals",
                               &history.solver_lin_rhs_evals},
                     std::pair{"solver_lin_setups", &history.solver_lin_setups},
                 }) {
                out[key] = py::array_t<double>(values->size(), values->data());
            }
            // Circuit history ("plasma_voltage" only; empty arrays without it).
            for (const auto& [key, values] : {
                     std::pair{"Rp", &history.Rp},
                     std::pair{"Rp_state", &history.Rp_state},
                     std::pair{"Vp", &history.Vp},
                     std::pair{"I", &history.I},
                     std::pair{"E", &history.E},
                 }) {
                out[key] = py::array_t<double>(values->size(), values->data());
            }
            out["species"] = history.species;
            out["nframes"] = history.nframes;
            out["npts"] = history.npts;
            out["nsp"] = history.nsp;
            return out; },
          // --- required: mechanism
          py::arg("mech"), py::arg("phase"),
          // --- required: grid
          py::arg("n_points"), py::arg("R_max"),
          // --- optional: grid
          py::arg("grid") = py::none(),
          // --- required: initial state
          py::arg("rho"), py::arg("Tg0"), py::arg("Te0"), py::arg("Y0"),
          // --- optional: initial radial profiles
          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("initY_r") = py::none(), py::arg("initY_v") = py::none(), py::arg("init_rho_r") = py::none(), py::arg("init_rho_v") = py::none(), py::arg("init_ur_r") = py::none(), py::arg("init_ur_v") = py::none(),
          // --- optional: chemistry
          py::arg("reacting") = true, py::arg("species_h0k") = py::none(),
          // --- optional: discharge
          py::arg("electric_field") = 0.0, py::arg("Efield_t") = py::none(), py::arg("Efield_v") = py::none(), py::arg("gap") = 0.0, py::arg("plasma_voltage") = py::none(),
          // --- optional: transport and collisions, bulk closures
          py::arg("sigma_T") = py::none(), py::arg("sigma_v") = py::none(), py::arg("nu_E") = 0.0, py::arg("nu_m") = 0.0, py::arg("nu_m_T") = py::none(), py::arg("nu_m_v") = py::none(), py::arg("kappa_H_T") = py::none(), py::arg("kappa_H_v") = py::none(), py::arg("kappa_e_T") = py::none(), py::arg("kappa_e_v") = py::none(), py::arg("D_species") = 0.0,
          // --- optional: transport and collisions, composition-resolved model
          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") = 200.0, py::arg("Te_max") = 1.0e5, py::arg("Te_n") = 400,
          // --- optional: ambient and boundary
          py::arg("T_amb") = 300.0, py::arg("p_amb") = -1.0, py::arg("Y_amb") = py::none(), py::arg("outer_bc_mode") = "ambient_ghost",
          // --- optional: hydrodynamics
          py::arg("hydrodynamics") = false, py::arg("cfl") = 0.9, py::arg("hydro_substep") = 0.0,
          // --- optional: time stepping
          py::arg("dt") = 1.0e-9, py::arg("n_steps") = 1000, py::arg("record_every") = 1,
          // --- optional: integrator controls
          py::arg("rtol") = 1.0e-7, py::arg("atol_T") = 1.0e-3, py::arg("atol_Y") = 1.0e-9, py::arg("independent_nodes") = false, py::arg("n_threads") = 1, py::arg("loglevel") = 0);
 
    // ---- The Cantera-backed thermo accessor, one state at a time -------
    m.def("hydro_thermo_probe",
          // Evaluate every ThermoAccessor method at one state.
          //
          // No time stepping and no fluxes, so the equation of state can be
          // asserted as the state function it is -- which is what NUMERICS.md
          // sections "Two-temperature sound speed" and "Energy-inversion round trip", cases 8 and 9 ask for.
          [](
              // Parameters
              // ----------
              const std::string &mech,
              //  Cantera mechanism (YAML) path
              const std::string &phase,
              //  PlasmaPhase name within the mechanism
              double rho,
              //  Mass density [kg/m^3]
              double Tg,
              //  Heavy-pool temperature [K]
              double Te,
              //  Electron temperature [K]
              py::array_t<double> Y,
              //  Mass fractions [-], length nsp, summing to 1
              py::object Tg_seed
              //  Seed for the heavy-temperature Newton solve [K]
              //  [Optional] By default, None (seeded at Tg)
              //
              // Returns
              // -------
              // dict
              //  p_heavy, p_electron [Pa]; sound_speed [m/s]; heavy_energy
              //  [J/kg of mixture]; electron_cv [J/kg of electrons/K];
              //  density_from_pressure [kg/m^3], which round-trips rho; and
              //  heavy_temperature [K], which round-trips Tg from Tg_seed.
          )
          {
            const std::vector<double> Yv = toVec(Y);
            rizer::CanteraThermo& thermo = cachedCanteraThermo(mech, phase);
            if (Yv.size() != thermo.nSpecies()) {
                throw std::invalid_argument(
                    "hydro_thermo_probe: Y has length "
                    + std::to_string(Yv.size()) + ", expected "
                    + std::to_string(thermo.nSpecies()) + ".");
            }
            double p_heavy = 0.0;
            double p_electron = 0.0;
            thermo.pressures(0.0, 0.0, rho, Tg, Te, Yv.data(), p_heavy, p_electron);
            const double heavy_energy = thermo.heavyEnergy(Tg, Yv.data());
            const double seed =
                Tg_seed.is_none() ? Tg : Tg_seed.cast<double>();
 
            py::dict out;
            out["p_heavy"] = p_heavy;
            out["p_electron"] = p_electron;
            out["sound_speed"] =
                thermo.soundSpeed(0.0, 0.0, rho, Tg, Te, Yv.data());
            out["heavy_energy"] = heavy_energy;
            out["electron_cv"] = thermo.electronCv();
            out["density_from_pressure"] = thermo.densityFromPressure(
                p_heavy + p_electron, Tg, Te, Yv.data());
            out["heavy_temperature"] =
                thermo.heavyTemperature(rho, Yv.data(), heavy_energy, seed);
            return out; },
          // --- required
          py::arg("mech"), py::arg("phase"), py::arg("rho"), py::arg("Tg"), py::arg("Te"), py::arg("Y"),
          // --- optional
          py::arg("Tg_seed") = py::none());
 
    // ---- One explicit hydrodynamic Strang half-step, plain arrays -------
    m.def("radial_diffusion_probe",
          // One node-by-node evaluation of the parabolic operator, nothing else.
          //
          // RadialDiffusion is not a sub-step: it runs inside CVODE's residual
          // and adds into rates the caller already holds, so there is no
          // "solve" to expose. This evaluates it on a state handed in, with
          // both capacitances set to 1, which makes every returned rate the
          // raw divergence -- W/m^3 for the two pools, kg/m^3/s per species.
          // That is what a manufactured solution can be compared against.
          //
          // A mechanism is still loaded, because the operator reads the
          // domain's grid, species count and per-node density. It runs no
          // chemistry: the probe never calls eval().
          [](
              // Parameters
              // ----------
              const std::string &mech,
              //     Cantera mechanism (YAML) path
              const std::string &phase,
              //     Phase name within the mechanism
              py::array_t<double> grid,
              //     Node radii [m], length N, ascending, grid[0] = 0
              py::array_t<double> rho,
              //     Mass density per node [kg/m^3], length N
              py::array_t<double> Tg,
              //     Heavy-pool temperature per node [K], length N
              py::array_t<double> Te,
              //     Electron temperature per node [K], length N
              py::array_t<double> Y,
              //     Mass fractions [-], shape [N, K], each row summing to 1
              double T_amb,
              //     Wall ghost temperature [K], the Dirichlet value both
              //     conduction calls sit against
              py::array_t<double> Y_amb,
              //     Wall ghost mass fractions [-], length K
              double wall_ghost_density,
              //     Mass density of the wall ghost half-cell [kg/m^3]
              py::object kappa_H_T,
              //     Heavy-pool conductivity abscissa [K], or None for no
              //     heavy conduction
              //     [Optional] By default, None
              py::object kappa_H_v,
              //     Heavy-pool conductivity [W/m/K] at kappa_H_T
              //     [Optional] By default, None
              py::object kappa_e_T,
              //     Electron conductivity abscissa [K], or None
              //     [Optional] By default, None
              py::object kappa_e_v,
              //     Electron conductivity [W/m/K] at kappa_e_T
              //     [Optional] By default, None
              double D_species
              //     Scalar species mass diffusivity [m^2/s]; 0 switches the
              //     species block and both energy terms it drives off
              //     [Optional] By default, 0.0
              //
              // Returns
              // -------
              // dict
              //     dTg, dTe [W/m^3] and dY [kg/m^3/s, shape (N, K)]: the
              //     divergences this node-by-node operator contributes.
          )
          {
            const auto g = toVec(grid);
            const auto rho_v = toVec(rho);
            const auto Tg_v = toVec(Tg);
            const auto Te_v = toVec(Te);
            const auto Y_v = toVec(Y);
            const auto Y_amb_v = toVec(Y_amb);
            const std::size_t n_points = g.size();
            if (n_points < 2) {
                throw std::invalid_argument(
                    "radial_diffusion_probe: the grid needs at least 2 nodes.");
            }
            requireLength(rho_v.size(), n_points, "rho");
            requireLength(Tg_v.size(), n_points, "Tg");
            requireLength(Te_v.size(), n_points, "Te");
            requireMatrixShape(Y, n_points, Y_amb_v.size(), "Y");
 
            rizer::PlasmaChannelConfig cfg;
            cfg.mech = mech;
            cfg.phase = phase;
            cfg.grid = g;
            cfg.n_points = n_points;
            cfg.R_max = g.back();
            cfg.rho = rho_v.front();
            // The initial profiles only have to make the constructor's state
            // valid; the probe overwrites x with the caller's field below. The
            // composition comes from that field's end rows, not from Y_amb: the
            // ambient may hold no electrons, and the constructor refuses a node
            // without them.
            cfg.initTg_r = {g.front(), g.back()};
            cfg.initTg_v = {Tg_v.front(), Tg_v.front()};
            cfg.initTe_r = {g.front(), g.back()};
            cfg.initTe_v = {Te_v.front(), Te_v.front()};
            const std::size_t n_species_in = Y_amb_v.size();
            if (Y_v.size() != n_points * n_species_in) {
                throw std::invalid_argument(
                    "radial_diffusion_probe: Y has " + std::to_string(Y_v.size()) +
                    " entries, expected " + std::to_string(n_points * n_species_in) +
                    " (one row of Y_amb's length per grid node).");
            }
            cfg.initY_r = {g.front(), g.back()};
            cfg.initY_v.assign(Y_v.begin(), Y_v.begin() + n_species_in);
            cfg.initY_v.insert(cfg.initY_v.end(), Y_v.end() - n_species_in, Y_v.end());
            cfg.Y_amb = Y_amb_v;
            cfg.T_amb = T_amb;
            cfg.reacting = false;
            cfg.D_species = D_species;
            cfg.kappa_H_T = optVec(kappa_H_T);
            cfg.kappa_H_v = optVec(kappa_H_v);
            cfg.kappa_e_T = optVec(kappa_e_T);
            cfg.kappa_e_v = optVec(kappa_e_v);
 
            rizer::PlasmaChannel1D channel(cfg);
            const std::size_t n_species = channel.nSpecies();
            if (Y_v.size() != n_points * n_species) {
                throw std::invalid_argument(
                    "radial_diffusion_probe: Y has " + std::to_string(Y_v.size()) +
                    " entries, expected " + std::to_string(n_points * n_species) + ".");
            }
            channel.setHydroState(rho_v, std::vector<double>(n_points, 0.0));
 
            // The domain's own component layout, filled by hand: this probe
            // never goes through eval(), so nothing else populates it.
            const std::size_t n_vars = 2 + n_species;
            std::vector<double> x(n_vars * n_points, 0.0);
            for (std::size_t j = 0; j < n_points; j++) {
                x[channel.index(0, j)] = Tg_v[j];
                x[channel.index(1, j)] = Te_v[j];
                for (std::size_t k = 0; k < n_species; k++) {
                    x[channel.index(2 + k, j)] = Y_v[j * n_species + k];
                }
            }
 
            rizer::RadialDiffusionCoefficients coefficients;
            coefficients.kappa = &channel.heavyConductivity();
            coefficients.kappa_e = &channel.electronConductivity();
            coefficients.species_diffusivity = D_species;
            coefficients.ambient_temperature = T_amb;
            coefficients.wall_ghost_density = wall_ghost_density;
            coefficients.Y_amb = &Y_amb_v;
            coefficients.molecular_weights = &channel.molecularWeights();
            coefficients.species_cp = &channel.speciesCp();
            coefficients.electron_index = channel.electronIndex();
            coefficients.n_species = n_species;
 
            std::vector<double> dTg(n_points, 0.0), dTe(n_points, 0.0);
            std::vector<double> dY(n_points * n_species, 0.0);
            std::vector<double> scratch(n_species, 0.0);
            for (std::size_t j = 0; j < n_points; j++) {
                // The enthalpy-diffusion term reads speciesCp(), a per-call
                // cache that eval() refreshes through sourceTerms(). Nothing
                // here calls sourceTerms, so this probe refreshes it itself;
                // without it the term would silently contribute zero.
                channel.updateSpeciesCp(Tg_v[j], Te_v[j],
                                        &Y_v[j * n_species], rho_v[j]);
                rizer::addRadialDiffusion(
                    channel, Cantera::span<const double>(x.data(), x.size()),
                    coefficients, j, /*lhsTg=*/1.0, /*lhsTe=*/1.0, dTg[j], dTe[j],
                    &dY[j * n_species], scratch.data());
            }
 
            py::dict out;
            out["dTg"] = py::array_t<double>(dTg.size(), dTg.data());
            out["dTe"] = py::array_t<double>(dTe.size(), dTe.data());
            out["dY"] = py::array_t<double>(
                {static_cast<py::ssize_t>(n_points),
                 static_cast<py::ssize_t>(n_species)},
                dY.data());
            return out; },
          // --- required
          py::arg("mech"),
          py::arg("phase"),
          py::arg("grid"),
          py::arg("rho"),
          py::arg("Tg"),
          py::arg("Te"),
          py::arg("Y"),
          py::arg("T_amb"),
          py::arg("Y_amb"),
          py::arg("wall_ghost_density"),
          // --- optional
          py::arg("kappa_H_T") = py::none(),
          py::arg("kappa_H_v") = py::none(),
          py::arg("kappa_e_T") = py::none(),
          py::arg("kappa_e_v") = py::none(),
          py::arg("D_species") = 0.0);
 
    m.def("solve_hydro_step_only",
          // Advance one explicit hydrodynamic Strang half-step on plain arrays.
          //
          // Solves the reduced hyperbolic system of PHYSICS.md section 5.4 --
          // continuity, radial momentum, species advection and the two pool
          // energies -- with the Rusanov flux and forward Euler of NUMERICS.md
          // section 5, sub-cycled to the acoustic CFL limit. No chemistry, no
          // Joule heating, no exchange, no conduction, no species diffusion.
          //
          // Deliberately separate from solve_1D_unsteady_discharge: it takes no
          // Cantera::Solution unless asked for one, so the pure-array unit tests
          // can cover the flux scheme, the boundary treatment and the axis
          // factor with nothing but numpy arrays and a gamma-law equation of
          // state.
          [](
              // Parameters
              // ----------
              //
              // Grid, from NUMERICS.md section 2
              py::array_t<double> r_node,
              //  Node radii [m], length N, ascending
              py::array_t<double> r_face,
              //  Face radii [m], length N+1, non-decreasing. r_face[0] = 0 makes
              //  the axis face vanish, which IS the symmetry condition;
              //  r_face = 1 everywhere is the planar entry point
              py::array_t<double> V_cell,
              //  Cell "volumes", cross-section per unit axial length [m^2],
              //  length N
              bool cylindrical_source,
              //  Apply the cylindrical geometric momentum source p*dr/V [-].
              //  MUST be false in planar mode: the face radii alone do not
              //  disable it
              //
              // State in, advanced in place
              py::array_t<double> rho,
              //  Mass density per node [kg/m^3], length N
              py::array_t<double> u_r,
              //  Radial velocity per node [m/s], length N
              py::array_t<double> Tg,
              //  Heavy-pool temperature per node [K], length N
              py::array_t<double> Te,
              //  Electron temperature per node [K], length N
              py::array_t<double> Y,
              //  Mass fractions [-], shape [N, K], each row summing to 1
              //
              // Equation of state
              const std::string &eos,
              //  Which ThermoAccessor closes the system. Three values:
              //
              //  "gamma" -- a calorically perfect gas-law accessor built from
              //    `gamma` and `molecular_weights`. Needs no mechanism and no
              //    Cantera object, so it is what the pure-array unit tests run
              //    on: the pressures, the sound speed and the energy inversion
              //    are all closed form, which makes the reference exact rather
              //    than iterative.
              //
              //  "cantera" -- the real two-temperature equation of state, read
              //    from `mech`/`phase` through Cantera's PlasmaPhase. Species
              //    thermo comes from the mechanism's NASA polynomials, so the
              //    heavy-pool inversion is a Newton solve rather than a
              //    division, and `electron_index` must be a real species. This
              //    is the production path, and the one the split solver uses.
              //
              //  "prescribed_p" -- TEST-ONLY. Pressure comes from the
              //    `prescribed_p` callable p(r, t) while every other property
              //    falls back to the gamma-law accessor, so a validation case
              //    can impose a closed-form pressure field and isolate the
              //    momentum row against it (NUMERICS.md section "Momentum alone", cases
              //    1-2). It needs `ghost_state`, because pressures() and
              //    densityFromPressure() are then no longer mutually inverse
              //    and neither outer mode could close a consistent ghost.
              //  [Optional] By default, "gamma"
              double gamma,
              //  Ratio of the HEAVY pool's specific heats,
              //  gamma_H = c_p,H / c_v,H, for eos="gamma" [-], > 1. The
              //  electron pool carries its own, fixed at 5/3
              //  [Optional] By default, 1.4
              py::array_t<double> molecular_weights,
              //  Per-species molecular weights for eos="gamma" [kg/kmol],
              //  length K
              //  [Optional] By default, empty
              int electron_index,
              //  Index of the electron species in Y [-], or -1 for none.
              //  Required by eos="cantera"
              //  [Optional] By default, -1
              py::object prescribed_p,
              //  Callable p(r, t) -> Pa for eos="prescribed_p"
              //  [Optional] By default, None
              const std::string &mech,
              //  Cantera mechanism (YAML) path for eos="cantera"
              //  [Optional] By default, "" (unset)
              const std::string &phase,
              //  Phase name within the mechanism for eos="cantera"
              //  [Optional] By default, "" (unset)
              //
              // Outer boundary, from NUMERICS.md section 7
              const std::string &outer_bc_mode,
              //  "ambient_ghost" writes the far-field state into the ghost and
              //  anchors the domain; "pressure_outlet" imposes p = p_amb and
              //  extrapolates the rest, and reflects
              //  [Optional] By default, "ambient_ghost"
              double p_amb,
              //  Ambient pressure closing the ghost density [Pa]
              //  [Optional] By default, 0.0
              double T_amb,
              //  Ambient temperature [K]
              //  [Optional] By default, 0.0
              py::array_t<double> Y_amb,
              //  Ambient mass fractions [-], length K, summing to 1; needed by
              //  "ambient_ghost"
              //  [Optional] By default, empty
              py::object ghost_state,
              //  TEST-ONLY callable
              //  (r_ghost, t, is_inner_face, rho, u_r, Tg, Te, Y) -> ghost state,
              //  overriding outer_bc_mode. `r_ghost` is the ghost CELL's node
              //  radius, one node spacing beyond the end node, so a closed-form
              //  field evaluated there fills the ghost where the ghost actually
              //  is. Deliberately absent from every
              //  production config: it needs an answer the solver is supposed to
              //  be computing (NUMERICS.md section 7.6)
              //  [Optional] By default, None
              //
              // Stepping
              double dt_total,
              //  Interval to advance over [s]; 0 probes the state and returns
              //  dt_cfl without advancing it
              std::size_t n_subcycles,
              //  Force this many equal sub-steps instead of CFL sub-cycling
              //  [-]; TEST-ONLY, for cases needing a known step size
              //  [Optional] By default, 0 (CFL)
              double cfl,
              //  Courant number for the sub-cycling [-], in (0, 1]
              //  [Optional] By default, 0.9
              double t0
              //  Time the interval starts at [s], passed to the accessors and
              //  the ghost callback
              //  [Optional] By default, 0.0
              //
              // Returns
              // -------
              // dict
              //  The advanced state (rho, u_r, Tg, Te, Y), the sub-cycle
              //  diagnostics (n_subcycles_used, dt_cfl, rho_min, p_min) and the
              //  accumulated boundary fluxes that
              //  close NUMERICS.md section 10's identities.
          )
          {
 
            rizer::HydroStepConfig cfg;
            cfg.r_node = toVec(r_node);
            cfg.r_face = toVec(r_face);
            cfg.V_cell = toVec(V_cell);
            cfg.cylindrical_source = cylindrical_source;
            cfg.p_amb = p_amb;
            cfg.T_amb = T_amb;
            cfg.Y_amb = toVec(Y_amb);
            cfg.dt_total = dt_total;
            cfg.n_subcycles = n_subcycles;
            cfg.cfl = cfl;
            cfg.t0 = t0;
            cfg.electron_index = electron_index;
 
            cfg.outer_bc_mode = rizer::parseOuterBc(outer_bc_mode);
 
            rizer::HydroStepState state;
            state.rho = toVec(rho);
            state.u_r = toVec(u_r);
            state.Tg = toVec(Tg);
            state.Te = toVec(Te);
            state.Y = toVec(Y);
            const std::size_t n_points = state.rho.size();
            if (n_points == 0) {
                throw std::invalid_argument(
                    "solve_hydro_step_only: rho is empty.");
            }
            requireMatrixShape(Y, n_points, std::nullopt, "Y");
            state.n_species = static_cast<std::size_t>(Y.shape(1));
 
            // Every accessor indexes Y by its OWN species list -- the
            // mechanism's for the Cantera one, molecular_weights for the
            // gamma-law one -- so a Y row shorter than that list is read past
            // its end, and for the last node past the array.
            const std::size_t n_weights =
                static_cast<std::size_t>(molecular_weights.size());
            if (eos != "cantera" && n_weights != state.n_species) {
                throw std::invalid_argument(
                    "solve_hydro_step_only: molecular_weights has "
                    + std::to_string(n_weights) + " entries but Y carries "
                    + std::to_string(state.n_species)
                    + " species per node.");
            }
 
            // Only the accessors this call constructs are owned here; the
            // Cantera one outlives the call in the cache above.
            std::unique_ptr<rizer::ThermoAccessor> owned_thermo;
            const rizer::ThermoAccessor* thermo = nullptr;
            if (eos == "cantera") {
                if (mech.empty() || phase.empty()) {
                    throw std::invalid_argument(
                        "solve_hydro_step_only: eos=\"cantera\" needs mech and "
                        "phase; the equation of state comes from the mechanism, "
                        "not from gamma and molecular_weights.");
                }
                rizer::CanteraThermo& cantera = cachedCanteraThermo(mech, phase);
                if (cantera.nSpecies() != state.n_species) {
                    throw std::invalid_argument(
                        "solve_hydro_step_only: Y carries "
                        + std::to_string(state.n_species)
                        + " species per node but phase \"" + phase
                        + "\" has " + std::to_string(cantera.nSpecies())
                        + ".");
                }
                if (electron_index < 0) {
                    // The two settings disagree about whether an electron pool
                    // exists: CanteraThermo keeps e- out of heavyEnergy() while
                    // electron_index = -1 skips the electron energy row, so the
                    // pool's internal energy would sit in neither and Te would
                    // be carried forward untouched. -1 is this binding's
                    // default, so the combination is reachable by omission.
                    throw std::invalid_argument(
                        "solve_hydro_step_only: eos=\"cantera\" needs a real "
                        "electron_index; -1 skips the electron energy row while "
                        "the Cantera accessor keeps the electron out of the "
                        "heavy pool, which drops its energy from both rows.");
                }
                if (static_cast<std::size_t>(electron_index) != cantera.electronIndex()) {
                    // The energy rows have to agree on which species is the
                    // electron: heavyEnergy() leaves out the phase's own
                    // electron, so any other index would run the electron row on
                    // one species while the real electron's energy sits in
                    // neither row.
                    throw std::invalid_argument(
                        "solve_hydro_step_only: electron_index = "
                        + std::to_string(electron_index)
                        + " but phase \"" + phase + "\" has its electron at index "
                        + std::to_string(cantera.electronIndex()) + ".");
                }
                thermo = &cantera;
            } else if (eos == "gamma") {
                owned_thermo = std::make_unique<rizer::GammaLawThermo>(
                    gamma, toVec(molecular_weights), electron_index);
                thermo = owned_thermo.get();
            } else if (eos == "prescribed_p") {
                if (prescribed_p.is_none()) {
                    throw std::invalid_argument(
                        "solve_hydro_step_only: eos=\"prescribed_p\" needs a "
                        "prescribed_p callable p(r, t) -> Pa.");
                }
                if (ghost_state.is_none()) {
                    // The prescribed accessor answers pressures() from p(r, t)
                    // but densityFromPressure() from its gamma-law backing, so
                    // the two are not mutually inverse. Either outer mode would
                    // then close the ghost density from p_amb while its
                    // pressure came back as p(r_face, t), discarding p_amb and
                    // leaving the ghost internally inconsistent.
                    throw std::invalid_argument(
                        "solve_hydro_step_only: eos=\"prescribed_p\" needs a "
                        "ghost_state callback. outer_bc_mode closes the ghost "
                        "density from p_amb through the gamma-law backing, "
                        "which does not invert the prescribed p(r, t).");
                }
                auto fn = prescribed_p.cast<py::function>();
                owned_thermo = std::make_unique<rizer::PrescribedPressureThermo>(
                    [fn](double r, double t) {
                        return fn(r, t).cast<double>();
                    },
                    rizer::GammaLawThermo(gamma, toVec(molecular_weights),
                                          electron_index));
                thermo = owned_thermo.get();
            } else {
                throw std::invalid_argument(
                    "solve_hydro_step_only: eos must be \"gamma\", \"cantera\" "
                    "or \"prescribed_p\", got \"" + eos + "\".");
            }
            cfg.thermo = thermo;
 
            if (!ghost_state.is_none()) {
                auto fn = ghost_state.cast<py::function>();
                const std::size_t nsp = state.n_species;
                cfg.ghost_state = [fn, nsp](double r_ghost, double t,
                                            bool is_inner_face,
                                            const rizer::NodeState& interior) {
                    py::object out = fn(
                        r_ghost, t, is_inner_face, interior.rho,
                        interior.u_r, interior.Tg, interior.Te,
                        py::array_t<double>(interior.Y.size(),
                                            interior.Y.data()));
                    auto tup = out.cast<py::tuple>();
                    if (tup.size() != 5) {
                        throw std::invalid_argument(
                            "solve_hydro_step_only: ghost_state must return "
                            "(rho, u_r, Tg, Te, Y).");
                    }
                    rizer::NodeState ghost;
                    ghost.rho = tup[0].cast<double>();
                    ghost.u_r = tup[1].cast<double>();
                    ghost.Tg = tup[2].cast<double>();
                    ghost.Te = tup[3].cast<double>();
                    ghost.Y = toVec(tup[4].cast<py::array_t<double>>());
                    if (ghost.Y.size() != nsp) {
                        throw std::invalid_argument(
                            "solve_hydro_step_only: ghost_state returned a Y "
                            "of the wrong length.");
                    }
                    return ghost;
                };
            }
 
            const rizer::HydroStepResult result = rizer::hydroStep(cfg, state);
 
            py::dict out;
            out["rho"] = py::array_t<double>(state.rho.size(), state.rho.data());
            out["u_r"] = py::array_t<double>(state.u_r.size(), state.u_r.data());
            out["Tg"] = py::array_t<double>(state.Tg.size(), state.Tg.data());
            out["Te"] = py::array_t<double>(state.Te.size(), state.Te.data());
            out["Y"] = py::array_t<double>(
                std::vector<py::ssize_t>{static_cast<py::ssize_t>(n_points),
                                         static_cast<py::ssize_t>(state.n_species)},
                state.Y.data());
            out["n_subcycles_used"] = result.n_subcycles_used;
            out["dt_cfl"] = result.dt_cfl;
            out["rho_min"] = result.rho_min;
            out["p_min"] = result.p_min;
            out["boundary_mass_out"] = result.boundary_mass_out;
            out["boundary_species_mass_out"] = py::array_t<double>(
                result.boundary_species_mass_out.size(),
                result.boundary_species_mass_out.data());
            out["boundary_total_energy_out"] = result.boundary_total_energy_out;
            out["boundary_electron_energy_out"] =
                result.boundary_electron_energy_out;
            out["pressure_work_electron"] = result.pressure_work_electron;
            return out; },
          // --- required: grid
          py::arg("r_node"), py::arg("r_face"), py::arg("V_cell"),
          // --- optional: grid
          py::arg("cylindrical_source") = true,
          // --- required: state in, advanced in place
          py::arg("rho"), py::arg("u_r"), py::arg("Tg"), py::arg("Te"), py::arg("Y"),
          // --- optional: equation of state
          py::arg("eos") = "gamma", py::arg("gamma") = 1.4, py::arg("molecular_weights") = py::array_t<double>(0), py::arg("electron_index") = -1, py::arg("prescribed_p") = py::none(), py::arg("mech") = "", py::arg("phase") = "",
          // --- optional: outer boundary
          py::arg("outer_bc_mode") = "ambient_ghost", py::arg("p_amb") = 0.0, py::arg("T_amb") = 0.0, py::arg("Y_amb") = py::array_t<double>(0), py::arg("ghost_state") = py::none(),
          // --- required: stepping
          py::arg("dt_total"),
          // --- optional: stepping
          py::arg("n_subcycles") = 0, py::arg("cfl") = 0.9, py::arg("t0") = 0.0);
}