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 "PropertyTable.h"
#include "PlasmaColumnSolver.h"
 
#include <vector>
 
namespace py = pybind11;
 
namespace {
 
// Build a PropertyTable from two equal-length numpy arrays (T strictly increasing).
// Transform numpy arrays into std::vector<double> for the PropertyTable constructor.
rizer::PropertyTable makeTable(const py::array_t<double>& T,
                               const py::array_t<double>& values)
{
    // Create unchecked accessors to the input arrays, which avoids bounds checking on every access.
    auto t = T.unchecked<1>();
    auto v = values.unchecked<1>();
    // Copy the data into std::vector<double> for the PropertyTable constructor.
    std::vector<double> tv(static_cast<std::size_t>(t.shape(0)));
    std::vector<double> vv(static_cast<std::size_t>(v.shape(0)));
    // Copy the data from the numpy arrays into the std::vectors.
    for (py::ssize_t i = 0; i < t.shape(0); i++) {
        tv[static_cast<std::size_t>(i)] = t(i);
    }
    for (py::ssize_t i = 0; i < v.shape(0); i++) {
        vv[static_cast<std::size_t>(i)] = v(i);
    }
    // Construct and return the PropertyTable from the vectors.
    return rizer::PropertyTable(std::move(tv), std::move(vv));
}
 
} // 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.");
 
    // 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.
    m.def("solve_column",
        [](double R,                    // Outer wall radius [m].
           double electric_field,       // Axial E field [V/m].
           double T_wall,               // Wall temperature [K].
           py::array_t<double> sigma_T, // Temperature values [K] — 1-D array, strictly increasing.
           py::array_t<double> sigma_v, // Corresponding sigma values [S/m] — 1-D array.
           py::array_t<double> kappa_T, // Temperature values [K] — 1-D array, strictly increasing.
           py::array_t<double> kappa_v, // Corresponding kappa values [W/(m·K)] — 1-D array.
           py::object prad_T,           // Temperature values [K] — 1-D array, strictly increasing, or None.
           py::object prad_v,           // Corresponding P_rad values [W/m³] — 1-D array, or None.
           std::size_t n_points,        // Initial uniform grid count.
           double T_center_guess,       // Parabolic-seed centerline temperature [K].
           double rho_cp,               // rho*cp [J/m³/K] for pseudo-transient scaling.
           int loglevel,                // Sim1D verbosity (0 = silent).
           bool refine_grid,            // Enable Cantera adaptive grid refinement.
           double refine_ratio,         // Refiner: max ratio of adjacent cell widths .
           double refine_slope,         // Refiner: max fractional slope change per cell.
           double refine_curve,         // Refiner: max fractional curvature per cell.
           double refine_prune,         // Refiner: remove points below this threshold.
           py::object init_r,           // Optional seed radii [m], strictly increasing, 0..R, or None.
           py::object init_T            // Optional seed temperatures [K] at init_r, or None.
        )
        {
            // 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).");
}