PlasmaChannelSolver.cpp

Go to the documentation of this file.

#include "PlasmaChannelSolver.h"
#include "PlasmaChannel1D.h"
#include "CircuitCoupling.h"
#include "HydroStep.h"
#include "thermo/CanteraThermo.h"
#include "thermo/ThermoAccessor.h"
#include "numeric/DiffusiveFlux.h"
#include "numeric/Format.h"
 
#include "cantera/oneD/Sim1D.h"
#include "cantera/oneD/Boundary1D.h"
#include "cantera/numerics/Integrator.h"
#include "cantera/numerics/FuncEval.h"
#include "cantera/base/ctexceptions.h"
 
#include <algorithm>
#include <cmath>
#include <cstdio>
#include <limits>
#include <sstream>
#include <functional>
#include <memory>
#include <stdexcept>
 
namespace rizer
{
 
    namespace
    {
 
        // -------------------------------------------------------------------------- //
        //   File-local helpers                                                       //
        // -------------------------------------------------------------------------- //
 
        // Method-of-lines RHS for the stiff (CVODE) backend:
        // @f[
        //     \dot{y} = f(y)
        // @f]
        // is the channel's spatial operator (transport + chemistry +
        // Joule/elastic). Because PlasmaChannel1D::eval returns each node's
        // @f$d\phi/dt@f$ (source/capacitance plus transport/capacitance) minus the
        // transient term @f$r_{dt}(\phi-\phi_{prev})@f$, calling it with
        // @f$r_{dt}=0@f$ yields exactly @f$\dot{y}@f$. The wall is a ghost-cell flux
        // (no algebraic constraint), so the whole system is a clean ODE.
        class ChannelRHS : public Cantera::FuncEval
        {
        public:
            ChannelRHS(PlasmaChannel1D *channel, Cantera::OneDim *sim,
                       const std::vector<double> &state, std::size_t n_vars,
                       std::size_t channel_offset, std::size_t n_points,
                       std::function<double(double)> electricFieldAt)
                : m_channel(channel), m_sim(sim), m_state(&state), m_n_vars(n_vars),
                  m_channel_offset(channel_offset), m_n_points(n_points),
                  m_field(std::move(electricFieldAt)) {}
 
            std::size_t neq() const override { return m_state->size(); }

            void getState(Cantera::span<double> y) override
            {
                std::copy(m_state->begin(), m_state->end(), y.begin());
            }

            void eval(double t,
                      Cantera::span<const double> y,
                      Cantera::span<double> ydot,
                      Cantera::span<const double> /*params*/) override
            {
                m_channel->setElectricField(m_field(t));
                // @f$r_{dt}=0@f$ -> the OneDim residual is exactly the
                // method-of-lines RHS @f$d\phi/dt@f$. Route through OneDim::eval
                // (the proven path that manages the mask workspace), not the
                // domain's eval directly.
                m_sim->eval(/*j=*/Cantera::npos, y, ydot, /*rdt=*/0.0, /*count=*/0);
            }

            void getConstraints(Cantera::span<double> constraints) override
            {
                std::fill(constraints.begin(), constraints.end(), 0.0);
                for (std::size_t j = 0; j < m_n_points; j++)
                {
                    for (std::size_t comp = 2; comp < m_n_vars; comp++)
                    {
                        // 1.0 => y_i >= 0
                        constraints[m_channel_offset + m_n_vars * j + comp] = 1.0;
                    }
                }
            }
 
        private:
            // Non-owning; owned by solveChannelTransient's `channel`
            PlasmaChannel1D *m_channel;
            // Non-owning; owned by solveChannelTransient's `sim`
            Cantera::OneDim *m_sim;
            // Non-owning; solveChannelTransient's live global state,
            // units per-component (Tg/Te [K], Y_k [-])
            const std::vector<double> *m_state;
            // Components per node, 2+K [-]
            std::size_t m_n_vars;
            // This domain's offset into the global state [-]
            std::size_t m_channel_offset;
            // Number of radial grid points [-]
            std::size_t m_n_points;
            // E(t) [V/m], t in [s]
            std::function<double(double)> m_field;
        };

        class NodeRHS : public Cantera::FuncEval
        {
        public:
            NodeRHS(PlasmaChannel1D *channel, const std::vector<double> &state,
                    std::size_t n_vars, std::size_t offset, std::size_t node,
                    std::function<double(double)> electricFieldAt)
                : m_channel(channel), m_state(&state), m_n_vars(n_vars),
                  m_offset(offset), m_node(node), m_field(std::move(electricFieldAt)) {}
 
            std::size_t neq() const override { return m_n_vars; }
 
            void getState(Cantera::span<double> y) override
            {
                std::copy(m_state->begin() + static_cast<std::ptrdiff_t>(m_offset),
                          m_state->begin() + static_cast<std::ptrdiff_t>(m_offset + m_n_vars),
                          y.begin());
            }
 
            void eval(double t, Cantera::span<const double> y,
                      Cantera::span<double> ydot,
                      Cantera::span<const double> /*params*/) override
            {
                m_channel->nodeRates(m_node, y[0], y[1], &y[2], m_field(t),
                                     ydot.data());
            }

            void getConstraints(Cantera::span<double> constraints) override
            {
                std::fill(constraints.begin(), constraints.end(), 0.0);
                for (std::size_t comp = 2; comp < m_n_vars; comp++)
                {
                    constraints[comp] = 1.0;
                }
            }
 
        private:
            PlasmaChannel1D *m_channel;         // Non-owning
            const std::vector<double> *m_state; // Non-owning
            std::size_t m_n_vars;
            std::size_t m_offset;
            std::size_t m_node;
            std::function<double(double)> m_field;
        };

        std::vector<double> faceRadii(const PlasmaChannel1D &channel)
        {
            const std::size_t n_points = channel.nPoints();
            std::vector<double> r_face(n_points + 1);
            for (std::size_t j = 0; j <= n_points; j++)
            {
                r_face[j] = channel.faceRadius(j);
            }
            return r_face;
        }

        std::vector<double> cellVolumes(const std::vector<double> &r_face)
        {
            std::vector<double> V_cell(r_face.size() - 1);
            for (std::size_t j = 0; j < V_cell.size(); j++)
            {
                V_cell[j] = numeric::annularVolume(r_face[j + 1], r_face[j]);
            }
            return V_cell;
        }

        HydroStepConfig buildHydroConfig(const PlasmaChannel1D &channel,
                                         const ChannelOptions &opts,
                                         const ThermoAccessor &thermo)
        {
            const std::size_t n_points = channel.nPoints();
            HydroStepConfig cfg;
            cfg.r_node.resize(n_points);
            for (std::size_t j = 0; j < n_points; j++)
                cfg.r_node[j] = channel.z(j);
            cfg.r_face = faceRadii(channel);
            cfg.V_cell = cellVolumes(cfg.r_face);
 
            cfg.cylindrical_source = true;
            cfg.outer_bc_mode = channel.outerBcMode();
            cfg.p_amb = opts.cfg.p_amb;
            cfg.T_amb = opts.cfg.T_amb;
            cfg.Y_amb = channel.ambientMassFractions();
            cfg.cfl = opts.cfg.cfl;
            if (opts.cfg.hydro_substep > 0.0)
            {
                const double count = 0.5 * opts.dt / opts.cfg.hydro_substep;
                const double whole = std::round(count);
                if (whole < 1.0 || std::abs(count - whole) > 1.0e-9 * count)
                {
                    throw std::invalid_argument(
                        "solveChannelTransient: the half-step dt/2 = " +
                        significant(0.5 * opts.dt) +
                        " s is not an integer multiple of hydro_substep = " +
                        significant(opts.cfg.hydro_substep) + " s; it comes to " +
                        significant(count) + " sub-steps.");
                }
                cfg.n_subcycles = static_cast<std::size_t>(whole);
                cfg.refuse_above_cfl = true;
            }
            cfg.electron_index = static_cast<int>(channel.electronIndex());
            cfg.thermo = &thermo;
            return cfg;
        }
 
        // One hydro half-step, reading and writing the two representations of the
        // same state: rho/u_r live on the domain, Tg/Te/Y live in the flat solution
        // vector the integrator owns. Both have to be current on the way in and on
        // the way out, which is the whole content of this function.
        void hydroHalfStep(PlasmaChannel1D &channel, HydroStepConfig &cfg,
                           std::vector<double> &x, std::size_t channel_offset,
                           std::size_t n_vars, double t0, double dt,
                           ChannelHistory &history)
        {
            const std::size_t n_points = channel.nPoints();
            const std::size_t n_species = channel.nSpecies();
 
            HydroStepState state;
            state.n_species = n_species;
            state.rho = channel.density();
            state.u_r = channel.velocity();
            state.Tg.resize(n_points);
            state.Te.resize(n_points);
            state.Y.resize(n_points * n_species);
            for (std::size_t j = 0; j < n_points; j++)
            {
                const std::size_t base = channel_offset + n_vars * j;
                state.Tg[j] = x[base + 0];
                state.Te[j] = x[base + 1];
                for (std::size_t k = 0; k < n_species; k++)
                {
                    state.Y[j * n_species + k] = x[base + 2 + k];
                }
            }
 
            // The only two fields that are not the same at every half-step.
            cfg.dt_total = dt;
            cfg.t0 = t0;
 
            const HydroStepResult result = hydroStep(cfg, state);
            history.boundary_mass_out += result.boundary_mass_out;
            for (std::size_t k = 0; k < n_species; k++)
            {
                history.boundary_species_mass_out[k] +=
                    result.boundary_species_mass_out[k];
            }
            history.boundary_total_energy_out += result.boundary_total_energy_out;
            history.boundary_electron_energy_out += result.boundary_electron_energy_out;
            history.pressure_work_electron += result.pressure_work_electron;
 
            channel.setHydroState(state.rho, state.u_r);
            for (std::size_t j = 0; j < n_points; j++)
            {
                const std::size_t base = channel_offset + n_vars * j;
                x[base + 0] = state.Tg[j];
                x[base + 1] = state.Te[j];
                for (std::size_t k = 0; k < n_species; k++)
                {
                    x[base + 2 + k] = state.Y[j * n_species + k];
                }
            }
        }

        std::string decodeSolverFailure(const std::string &message,
                                        std::size_t failed_node,
                                        const PlasmaChannel1D &channel,
                                        std::size_t channel_offset,
                                        std::size_t n_vars,
                                        std::size_t n_points)
        {
            // Failure time [s], if the message carries it
            // ("At t = <x> and h = ...").
            double t_fail = -1.0;
            // Character offset into `message` [-]
            std::size_t t_pos = message.find("At t = ");
            if (t_pos != std::string::npos)
            {
                try
                {
                    t_fail = std::stod(message.substr(t_pos + 7));
                }
                catch (...)
                {
                }
            }
            std::ostringstream out;
            out << "PlasmaChannel1D BDF solve failed to converge";
            if (t_fail >= 0.0)
                out << " near t=" << t_fail << " s";
            out << ".\nThe implicit corrector could not resolve the local "
                   "dynamics and the step size hit its floor. The dominant "
                   "error contributions are:\n";
            std::istringstream in(message);
            std::string line;
            bool in_error_block = false, found_any = false;
            while (std::getline(in, line))
            {
                if (line.find("Components with largest") != std::string::npos)
                {
                    in_error_block = true;
                    continue;
                }
                if (!in_error_block)
                    continue;
                std::size_t colon = line.find(':');
                if (colon == std::string::npos)
                    continue;
                long global_index;
                double error_value;
                try
                {
                    global_index = std::stol(line.substr(0, colon));
                    error_value = std::stod(line.substr(colon + 1));
                }
                catch (...)
                {
                    continue;
                }
                // A single-node system numbers its components from zero and
                // carries no node in the index, so the node comes from which
                // system failed rather than from the message.
                const long base =
                    (failed_node == Cantera::npos)
                        ? static_cast<long>(channel_offset)
                        : 0;
                long local_index = global_index - base;
                if (local_index < 0)
                    continue;
                std::size_t j = (failed_node == Cantera::npos)
                                    ? static_cast<std::size_t>(local_index) / n_vars
                                    : failed_node;
                std::size_t comp = static_cast<std::size_t>(local_index) % n_vars;
                if (j >= n_points)
                    continue;
                out << "  node " << j << " (r=" << channel.z(j) << " m)  "
                    << channel.componentName(comp)
                    << "  (weighted error " << error_value << ")\n";
                found_any = true;
            }
            if (!found_any)
            {
                out << message << "\n";
            }
            else
            {
                // The states the closures refused (TrialStateError and the like)
                // are why the step size fell; dropping them with the component
                // table would leave the caller a convergence failure with no
                // cause.
                const std::string caught_marker = "Exceptions caught during RHS evaluation";
                const std::size_t caught = message.find(caught_marker);
                if (caught != std::string::npos)
                {
                    const std::size_t table =
                        message.find("Components with largest", caught);
                    out << message.substr(caught, table == std::string::npos
                                                      ? std::string::npos
                                                      : table - caught)
                        << "\n";
                }
            }
            out << "Remedies: reduce dt, loosen tolerances, soften the initial "
                   "profiles or the field ramp, or verify the initial state is "
                   "physically consistent at the flagged nodes.";
            return out.str();
        }
 
    } // namespace
 
    // ------------------------------------------------------------------------------ //
    //   Free functions                                                               //
    // ------------------------------------------------------------------------------ //
 
    ChannelHistory solveChannelTransient(const ChannelOptions &opts)
    {
        using std::make_shared;
 
        // Ten steps, in order:
        //    1. build the domain and bracket it into a Sim1D
        //    2. mirror the domain's state into the flat vector the integrator owns
        //    3. the history, and the recorders that append one frame to it
        //    4. the applied field: a prescribed schedule, or a closed circuit
        //    5. reject option combinations that cannot mean what they say
        //    6. the hydro sub-step's own geometry and equation of state
        //    7. the t=0 frame, recorded before anything is advanced
        //    8. the CVODE sub-systems: one coupled, or one per node
        //    9. graded startup, on the path that integrates straight through
        //   10. the frame loop, then the recorded history out
 
        // ---- 1. domain and Sim1D --------------------------------------------
        auto channel = make_shared<PlasmaChannel1D>(opts.cfg);
 
        auto left = make_shared<Cantera::Empty1D>();
        auto right = make_shared<Cantera::Empty1D>();
        std::vector<std::shared_ptr<Cantera::Domain1D>> domains{left, channel, right};
        Cantera::Sim1D sim(domains);
        sim.getInitialSoln(); // Apply PlasmaChannel1D::initialValue
 
        const std::size_t n_species = channel->nSpecies();
        const std::size_t n_points = channel->nPoints();
        const std::size_t n_vars = 2 + n_species; // Per-node: Tg, Te, Y_1..Y_K
        const std::size_t channel_offset =
            channel->loc(); // This domain's offset in the global state
 
        // ---- 2. the flat state the integrator owns --------------------------
        // Flat global state vector, units per-component (Tg/Te [K], Y_k [-]), length
        // channel_offset + n_vars*n_points. The integrator owns it and passes it to
        // ChannelRHS::eval() at every internal step, so it must be kept current
        // with the domain's own state (PlasmaChannel1D::m_state) and vice versa.
        std::vector<double> x(sim.size(), 0.0);
 
        // Initialize x from the populated solution (getInitialSoln wrote it into
        // the internal m_state). NOTE: getState() copies the Newton work buffer
        // m_xnew (empty here), so read the real state via the domain's values()
        // (Domain1D::getValues reads m_state).
        for (std::size_t comp = 0; comp < n_vars; comp++)
        {
            std::vector<double> vals = channel->values(channel->componentName(comp));
            for (std::size_t j = 0; j < n_points; j++)
            {
                x[channel_offset + n_vars * j + comp] = vals[j];
            }
        }
 
        // ---- 3. the history and its recorders -------------------------------
        // Static metadata (grid, species names, sizes) is fixed for the whole run
        // and only needs setting once, up front; only t/Tg/Te/ne/Y grow per frame.
        ChannelHistory history;
        history.npts = n_points;
        history.nsp = n_species;
        history.boundary_species_mass_out.assign(n_species, 0.0);
        for (std::size_t j = 0; j < n_points; j++)
            history.r.push_back(channel->z(j));
        for (std::size_t k = 0; k < n_species; k++)
        {
            history.species.push_back(channel->componentName(2 + k));
        }
 
        // Appends one frame (all nodes) to `history` from the CURRENT flat state
        // `x`, in the flat layout x[channel_offset + n_vars*j + component].
        // Captured by reference so the loop below can call it at its own cadence.
        auto record = [&](double t)
        {
            history.t.push_back(t);
            for (std::size_t j = 0; j < n_points; j++)
            {
                const double Tg = x[channel_offset + n_vars * j + 0];
                const double Te = x[channel_offset + n_vars * j + 1];
                history.Tg.push_back(Tg);
                history.Te.push_back(Te);
                double pressure, electron_density; // [Pa], [1/m^3]
                channel->nodeDiagnostics(Tg, Te,
                                         &x[channel_offset + n_vars * j + 2], j,
                                         pressure, electron_density);
                history.ne.push_back(electron_density);
                for (std::size_t k = 0; k < n_species; k++)
                {
                    history.Y.push_back(x[channel_offset + n_vars * j + 2 + k]);
                }
                history.rho.push_back(channel->density()[j]);
                history.p.push_back(pressure);
                history.u_r.push_back(channel->velocity()[j]);
            }
        };
 
        // ---- 4. the applied field: schedule, or circuit ---------------------
        // Field schedule E(t): an explicit table if provided, else the constant.
        // A malformed one is refused rather than dropped: silently falling back to
        // the constant would run the whole discharge at the wrong field and report
        // nothing.
        PropertyTable Efield;
        // Round-off allowance on the schedule's end [s]: a caller computing
        // t_end + dt and this solver's n_steps * dt + dt can differ in the last
        // bit.
        const double schedule_tolerance =
            1.0e-9 * opts.dt * static_cast<double>(opts.n_steps + 1);
        if (!opts.Et_t.empty() || !opts.Et_v.empty())
        {
            if (opts.Et_t.size() < 2 || opts.Et_t.size() != opts.Et_v.size())
            {
                throw std::invalid_argument(
                    "solveChannel: the E(t) schedule needs at least 2 points and "
                    "Et_t and Et_v of equal length; got " +
                    std::to_string(opts.Et_t.size()) + " times and " +
                    std::to_string(opts.Et_v.size()) + " values.");
            }
            // CVODE's internal step is capped at dt (setMaxStepSize below), so it
            // can evaluate the field up to one step past the last output time.
            // Outside the schedule the table would hold an end value.
            const double run_end = opts.dt * static_cast<double>(opts.n_steps);
            if (opts.Et_t.front() > 0.0 ||
                opts.Et_t.back() < run_end + opts.dt - schedule_tolerance)
            {
                throw std::invalid_argument(
                    "solveChannel: the E(t) schedule Et_t covers [" +
                    significant(opts.Et_t.front()) + ", " +
                    significant(opts.Et_t.back()) + "] s but has to cover [0, " +
                    significant(run_end + opts.dt) +
                    "] s: the run's n_steps * dt plus one step dt, the farthest "
                    "CVODE evaluates the field.");
            }
            Efield = PropertyTable(opts.Et_t, opts.Et_v);
        }
 
        // Circuit coupling (ChannelOptions::plasma_voltage).
        // Inactive without a callback, in which case every method below is a
        // no-op and the applied field comes from the prescribed schedule.
        ChannelCircuit circuit(*channel, opts.gap, opts.plasma_voltage,
                               channel_offset, n_vars);

        auto electricFieldAt = [&](double t)
        {
            if (circuit.active())
            {
                return circuit.electricField(t);
            }
            if (Efield.empty())
            {
                return opts.cfg.electric_field;
            }
            if (t < Efield.first() || t > Efield.last() + schedule_tolerance)
            {
                throw std::runtime_error(
                    "solveChannel: the field is evaluated at t = " + significant(t) +
                    " s, outside the E(t) schedule [" + significant(Efield.first()) +
                    ", " + significant(Efield.last()) + "] s.");
            }
            return Efield.eval(t);
        };

        auto recordCircuit = [&]()
        {
            if (!circuit.active())
                return;
            const CircuitSample circuit_sample = circuit.sample(x);
            history.Rp.push_back(circuit_sample.Rp);
            history.Rp_state.push_back(circuit_sample.Rp_state);
            history.Vp.push_back(circuit_sample.Vp);
            history.I.push_back(circuit_sample.current);
            history.E.push_back(circuit_sample.electric_field);
        };
 
        // ---- 5. refuse option combinations that cannot hold -----------------
        // Under hydrodynamics the Strang step is opts.dt, so a frame boundary
        // that is not a whole number of steps would land inside one. Without the
        // hydro sub-step nothing depends on where the boundary falls, so the
        // check is scoped to the mode that needs it.
        // n_steps == 0 would otherwise still run: n_intervals is floored at 1 and
        // the inner loop takes record_every split steps regardless, so the hydro
        // path would advance time the caller never asked for.
        // std:: exceptions with concatenated text, not CanteraError's {} form:
        // fmt is compiled into cantera_shared rather than linked separately, so
        // formatting here would couple this translation unit to a library
        // CMakeLists.txt only treats as optional (.claude/rules/cantera-extension.md).
        if (opts.independent_nodes && channel->hasTransport())
        {
            throw std::runtime_error(
                "solveChannelTransient: independent_nodes solves each node as its "
                "own ODE system, which is only equivalent to the coupled solve "
                "when no transport term couples the nodes. This channel has "
                "conduction or species diffusion enabled, so its nodes are not "
                "independent. Leave independent_nodes off.");
        }
        if (!(opts.cfg.hydro_substep >= 0.0) ||
            (opts.cfg.hydro_substep > 0.0 && !opts.cfg.hydrodynamics))
        {
            throw std::invalid_argument(
                "solveChannelTransient: hydro_substep must be >= 0 s and is used only "
                "with hydrodynamics; got " + significant(opts.cfg.hydro_substep) +
                " s with hydrodynamics " +
                (opts.cfg.hydrodynamics ? "on." : "off."));
        }
        if (opts.plasma_voltage &&
            (opts.cfg.electric_field != 0.0 || !opts.Et_t.empty() || !opts.Et_v.empty()))
        {
            throw std::invalid_argument(
                "solveChannelTransient: plasma_voltage sets the field from the circuit, "
                "so a prescribed electric_field or Et_t/Et_v schedule would be ignored. "
                "Give one or the other.");
        }
        if (opts.plasma_voltage && opts.n_threads > 1)
        {
            throw std::runtime_error(
                "solveChannelTransient: plasma_voltage is evaluated once per "
                "global step from the solver's own thread and would have to "
                "re-acquire the GIL from each worker otherwise. Use n_threads=1 "
                "with a circuit callback.");
        }
        if (opts.n_threads < 1)
        {
            throw std::runtime_error(
                "solveChannelTransient: n_threads must be at least 1.");
        }
        if (opts.n_threads > 1 && !opts.independent_nodes)
        {
            throw std::runtime_error(
                "solveChannelTransient: n_threads=" +
                std::to_string(opts.n_threads) +
                " needs independent_nodes: the coupled solve is one banded "
                "system with no per-node loop to spread across threads. "
                "Rejected rather than silently ignored.");
        }
        if (opts.rtol <= 0.0 || opts.atol_T <= 0.0 || opts.atol_Y <= 0.0)
        {
            std::ostringstream message;
            message << "solveChannelTransient: the bdf tolerances must all be "
                       "positive; rtol="
                    << opts.rtol << ", atol_T=" << opts.atol_T
                    << ", atol_Y=" << opts.atol_Y << ".";
            throw std::runtime_error(message.str());
        }
 
        // Hydrodynamics and a closed circuit both act once per global step: the
        // hydro half-steps bracket each one, and the circuit freezes R_p over
        // each one. Both therefore advance step by step inside a frame, and a
        // frame has to be a whole number of steps.
        const bool stepped = opts.cfg.hydrodynamics || circuit.active();
        if (stepped && (opts.record_every == 0 || opts.n_steps == 0 ||
                        opts.n_steps % opts.record_every != 0))
        {
            throw std::runtime_error(
                "solveChannelTransient: hydrodynamics and plasma_voltage advance "
                "one global step at a time, so they need n_steps > 0 and "
                "record_every to divide it exactly; n_steps=" +
                std::to_string(opts.n_steps) + " and record_every=" +
                std::to_string(opts.record_every) + ".");
        }
 
        // ---- 6. the hydro sub-step's geometry and EOS -----------------------
        // The geometry and the equation of state the hydro sub-step runs on,
        // and how many split steps make up one frame.
        const std::unique_ptr<CanteraThermo> hydro_thermo =
            opts.cfg.hydrodynamics
                ? std::make_unique<CanteraThermo>(opts.cfg.mech, opts.cfg.phase)
                : nullptr;
        HydroStepConfig hydro_cfg =
            opts.cfg.hydrodynamics
                ? buildHydroConfig(*channel, opts, *hydro_thermo)
                : HydroStepConfig{};
        const std::size_t steps_per_frame = stepped ? opts.record_every : 1;
 
        // ---- 7. the t=0 frame -----------------------------------------------
        // The t=0 frame has had no step driven through it, so the resistance
        // recorded against it is the initial state's and the voltage is zero.
        circuit.seedResistance(x);
        record(/*t=*/0.0);
        recordCircuit();
 
        // Fixed output times, uniform over the total interval
        // @f$dt \cdot n_\text{steps}@f$.
        const double t_total = // Total simulated time [s]
            opts.dt * static_cast<double>(opts.n_steps);
        const std::size_t n_intervals = // Number of recorded frames after t=0 [-]
            std::max<std::size_t>(
                1, opts.n_steps / std::max<std::size_t>(1, opts.record_every));
 
        // ---- 8. the CVODE sub-systems ---------------------------------------
        // One sub-system under "coupled", N of them under "independent". The
        // loops below never ask which: they advance, re-seed and read the
        // counters of whatever is in here. `offset` and `length` say which
        // slice of `x` each one owns, which is the whole domain or one node.
        struct SubSystem
        {
            std::unique_ptr<Cantera::FuncEval> rhs;
            std::unique_ptr<Cantera::Integrator> cvode;
            std::size_t offset = 0; 
            std::size_t length = 0; 
            std::size_t node = Cantera::npos;
        };
        const bool independent = opts.independent_nodes;
        // More threads than nodes would only idle, and each one costs a
        // mechanism parse.
        const int threads = independent
                                ? static_cast<int>(std::min(opts.n_threads, n_points))
                                : 1;
        if (threads > 1)
        {
            channel->setEvaluatorPoolSize(static_cast<std::size_t>(threads));
        }
        std::vector<SubSystem> systems;
        systems.reserve(independent ? n_points : 1);
        if (independent)
        {
            for (std::size_t j = 0; j < n_points; j++)
            {
                SubSystem system;
                system.offset = channel_offset + n_vars * j;
                system.length = n_vars;
                system.node = j;
                system.rhs = std::make_unique<NodeRHS>(
                    channel.get(), x, n_vars, system.offset, j,
                    electricFieldAt);
                systems.push_back(std::move(system));
            }
        }
        else
        {
            SubSystem system;
            system.offset = 0;
            system.length = x.size();
            system.rhs = std::make_unique<ChannelRHS>(
                channel.get(), &sim, x, n_vars, channel_offset, n_points,
                electricFieldAt);
            systems.push_back(std::move(system));
        }
 
        for (SubSystem &system : systems)
        {
            system.cvode.reset(Cantera::newIntegrator("CVODE"));
            Cantera::Integrator *cvode = system.cvode.get();
            cvode->setMethod(Cantera::BDF_Method);
            cvode->setLinearSolverType("BAND");
            // With conduction and diffusion off, every node's residual depends on
            // that node alone, so the Jacobian is block diagonal with blocks of
            // n_vars and its half-bandwidth is n_vars-1 rather than the 2*n_vars
            // a j +/- 1 stencil needs. The difference is not cosmetic: CVODE
            // forms this Jacobian by difference quotients, one right-hand side
            // evaluation per column group, so the narrower band halves the groups
            // and with them the kinetics evaluations that dominate a restart.
            // Measured at 1.83x on a decoupled breakdown run, to the same answer
            // in every digit.
            // A single node's system is n_vars wide and densely coupled within
            // itself, so n_vars-1 is a full band there rather than a narrowing;
            // it is the same number either way.
            const int bandwidth = static_cast<int>(
                channel->hasTransport() ? 2 * n_vars : n_vars - 1);
            cvode->setBandwidth(bandwidth, bandwidth);
            // Per-component absolute tolerances: K-scale for T, mass-fraction
            // floor for Y_k. See ChannelOptions for what each one satisfies.
            std::vector<double> atol(system.length);
            for (std::size_t j = 0; j < (independent ? 1 : n_points); j++)
            {
                atol[n_vars * j + 0] = opts.atol_T; // Tg [K]
                atol[n_vars * j + 1] = opts.atol_T; // Te [K]
                for (std::size_t k = 0; k < n_species; k++)
                {
                    atol[n_vars * j + 2 + k] = opts.atol_Y; // Y_k [-]
                }
            }
            cvode->setTolerances(opts.rtol,
                                 Cantera::span<const double>(atol.data(), atol.size()));
            // Cap the internal step at the output cadence. This is essential
            // for the stiff plasma electron-energy equation (tiny electron
            // heat capacity -> sub-ns Te dynamics): without it CVODE attempts
            // an over-large first step and stalls. The 0D NRP reactor relies
            // on the same max_step setting.
            cvode->setMaxStepSize(opts.dt);
            // A sane step budget: high enough for a long stiff run, low enough
            // that a genuine stall throws a clean (decoded) error instead of
            // spinning at a denormal step until SUNDIALS corrupts (the
            // access-violation we saw).
            cvode->setMaxSteps(500000);
            // A trial state the closures reject throws a recoverable
            // CanteraError (TrialStateError), on which CVODE retries a smaller
            // step. Suppressed rather than logged, so each retry does not print;
            // CVODE folds the collected messages into its own error if the step
            // finally fails.
            system.rhs->suppressErrors(true);
            cvode->initialize(/*t0=*/0.0, *system.rhs);
        }
 
        // Advance the reaction sub-step to `t1`, decoding a CVODE failure
        // into the node and component it belongs to.
        // An exception cannot cross an OpenMP structured block -- that is
        // undefined behaviour and terminates the process -- so a failing node
        // records its decoded message and lets the others finish; the first
        // failure in node order is rethrown once the loop is done. The nodes
        // write disjoint slices of `x`, so nothing else needs guarding.
        //
        // EVERY exception type has to be caught, not just CanteraError, and
        // the catch body must not throw either: `if (thread_count > 1)` only
        // shrinks the team to one thread, it does not remove the construct,
        // so the coupled solve runs inside it too. A std::bad_alloc out
        // of SUNDIALS, or out of decodeSolverFailure's own string building,
        // would otherwise abort the process with no diagnostic instead of
        // raising in Python.
        auto integrateTo = [&](double t1)
        {
            // Plain locals, not captures: MSVC implements OpenMP 2.0, whose
            // clauses cannot name a variable a lambda captured, and whose
            // parallel-for index must be a signed int.
            const int count = static_cast<int>(systems.size());
            const int thread_count = threads;
            std::vector<std::string> failures(systems.size());
            // Allocated before the region so that recording a failure is a
            // byte store that cannot itself throw. The decoded text is a
            // best effort on top of it.
            std::vector<char> failed(systems.size(), 0);
#ifdef _OPENMP
#pragma omp parallel for schedule(dynamic) num_threads(thread_count) if (thread_count > 1)
#endif
            for (int i = 0; i < count; i++)
            {
                SubSystem &system = systems[static_cast<std::size_t>(i)];
                try
                {
                    system.cvode->integrate(t1);
                    auto solution = system.cvode->solution();
                    std::copy(
                        solution.begin(), solution.end(),
                        x.begin() + static_cast<std::ptrdiff_t>(system.offset));
                }
                catch (const std::exception &err)
                {
                    failed[static_cast<std::size_t>(i)] = 1;
                    try
                    {
                        // Format-free std::string
                        failures[static_cast<std::size_t>(i)] =
                            decodeSolverFailure(err.what(), system.node, *channel,
                                                channel_offset, n_vars,
                                                n_points);
                    }
                    catch (...)
                    {
                        // Decoding is a convenience; losing it must not cost
                        // the diagnostic the flag above already recorded.
                    }
                }
                catch (...)
                {
                    failed[static_cast<std::size_t>(i)] = 1;
                }
            }
            for (std::size_t i = 0; i < failed.size(); i++)
            {
                if (!failed[i])
                    continue;
                // Concatenated rather than formatted: the CVODE text can
                // contain braces, and this translation unit does not pull in
                // fmt (.claude/rules/cantera-extension.md).
                throw std::runtime_error(
                    "solveChannelTransient: " +
                    (failures[i].empty()
                         ? "the reaction sub-step failed and the error could "
                           "not be decoded"
                         : failures[i]));
            }
        };
 
        // Advance to `t1` like integrateTo(), and accumulate what conduction
        // and species diffusion moved out through the wall face meanwhile, by
        // the trapezoidal rule over the interval. Both ends are read from `x`
        // as it stands, so an interval that opens after a hydro half-step
        // starts from what that half-step wrote. Only times the caller was
        // already going to stop at are sampled, so CVODE's step sequence is
        // the one integrateTo() alone would produce.
        double t_advanced = 0.0; // Time the reaction sub-step last reached [s]
        auto advanceTo = [&](double t1)
        {
            if (!channel->hasTransport())
            {
                integrateTo(t1);
                t_advanced = t1;
                return;
            }
            const Cantera::span<const double> x_channel(x.data() + channel_offset,
                                                        n_vars * n_points);
            double conduction_start, diffusion_start; // [W/m]
            std::vector<double> species_start, species_end; // [kg/m/s]
            channel->wallEnergyOutflow(x_channel, conduction_start, diffusion_start,
                                       species_start);
            integrateTo(t1);
            double conduction_end, diffusion_end; // [W/m]
            channel->wallEnergyOutflow(x_channel, conduction_end, diffusion_end,
                                       species_end);
            const double half_interval = 0.5 * (t1 - t_advanced); // [s]
            history.boundary_conduction_energy_out +=
                half_interval * (conduction_start + conduction_end);
            history.boundary_diffusion_energy_out +=
                half_interval * (diffusion_start + diffusion_end);
            for (std::size_t k = 0; k < n_species; k++)
            {
                history.boundary_species_mass_out[k] +=
                    half_interval * (species_start[k] + species_end[k]);
            }
            t_advanced = t1;
        };
 
        // ---- 9. graded startup ----------------------------------------------
        // Graded startup. The 2T electron-energy + plasma-chemistry layer is
        // ferociously stiff at t=0 (tiny electron heat capacity ->
        // @f$dT_e/dt \sim 10^{15}@f$ K/s once the field is on). The proven 0D
        // reactor tames this with scipy vode's first_step=1e-15; Cantera's
        // CVODES exposes no initial-step setter, but CVODE caps its first
        // internal step at (tout - t0), so integrating to a sequence of tiny
        // times forces a small, safe first step that then grows under error
        // control -- the same effect, from the caller side.
        // The Strang loop below re-seeds the integrator at every split step
        // and would discard whatever this advanced, so it runs only on the
        // path that integrates straight through.
        // Warm-up integration times [s]
        //
        // With a circuit, close it over the first step BEFORE warming up.
        // electricFieldAt() reads the circuit's two voltages once it is
        // active, and they are still their zero initialisers until the first
        // closeOver() call -- so an unseeded warm-up would advance to
        // 1e-10 s with no applied field at all, which on a 1 ns cadence is a
        // tenth of the opening step integrating the wrong problem. The step
        // loop below skips its own call for the first step rather than
        // querying a stateful circuit twice for the same window.
        //
        // The warm-up stops short of the first time the loop below advances
        // to -- the first step's end, or the first frame's -- because CVODE
        // cannot integrate backwards to a time it has already passed.
        const double first_target =
            stepped ? opts.dt : t_total / static_cast<double>(n_intervals);
        bool circuit_seeded = false;
        if (!opts.cfg.hydrodynamics)
        {
            if (circuit.active())
            {
                circuit.closeOver(/*t0=*/0.0, /*t1=*/opts.dt, x);
                circuit_seeded = true;
            }
            for (double t_warm :
                 {1.0e-15, 1.0e-14, 1.0e-13, 1.0e-12, 1.0e-11, 1.0e-10})
            {
                if (t_warm >= first_target)
                    break;
                advanceTo(t_warm);
            }
        }
 
        // CVODE's own counters for one frame, summed over the sub-systems.
        // Under "independent" these are summed over the nodes, except
        // `last_order`, which is the worst (lowest) order any node finished
        // on: a mean would hide the one node that never climbed back.
        struct SolverCounters
        {
            double steps = 0.0, rhs_evals = 0.0, jac_evals = 0.0;
            double err_test_fails = 0.0, nonlin_conv_fails = 0.0;
            double lin_rhs_evals = 0.0, lin_setups = 0.0;
            double last_order = 0.0;
        };
        auto readSolverCounters = [&]()
        {
            SolverCounters counters;
            counters.last_order = std::numeric_limits<double>::max();
            for (const SubSystem &system : systems)
            {
                const Cantera::AnyMap stats = system.cvode->solverStats();
                counters.steps += static_cast<double>(stats.getInt("steps", 0));
                counters.rhs_evals +=
                    static_cast<double>(stats.getInt("rhs_evals", 0));
                counters.jac_evals +=
                    static_cast<double>(stats.getInt("jac_evals", 0));
                counters.err_test_fails +=
                    static_cast<double>(stats.getInt("err_test_fails", 0));
                counters.nonlin_conv_fails += static_cast<double>(
                    stats.getInt("nonlinear_conv_fails", 0));
                counters.lin_rhs_evals +=
                    static_cast<double>(stats.getInt("lin_rhs_evals", 0));
                counters.lin_setups +=
                    static_cast<double>(stats.getInt("lin_solve_setups", 0));
                counters.last_order = std::min(
                    counters.last_order,
                    static_cast<double>(stats.getInt("last_order", 0)));
            }
            return counters;
        };
        auto recordSolverCounters = [&](const SolverCounters &counters)
        {
            history.solver_steps.push_back(counters.steps);
            history.solver_rhs_evals.push_back(counters.rhs_evals);
            history.solver_jac_evals.push_back(counters.jac_evals);
            history.solver_err_test_fails.push_back(counters.err_test_fails);
            history.solver_nonlin_conv_fails.push_back(counters.nonlin_conv_fails);
            history.solver_last_order.push_back(counters.last_order);
            history.solver_lin_rhs_evals.push_back(counters.lin_rhs_evals);
            history.solver_lin_setups.push_back(counters.lin_setups);
        };
        // The t=0 frame: nothing has been integrated, so every count is zero,
        // and the counters stay one entry per recorded frame.
        recordSolverCounters(SolverCounters{});
 
        // ---- 10. the frame loop --------------------------------------------
        for (std::size_t frame = 1; frame <= n_intervals; frame++)
        {
            const double t_rec = t_total *
                                 static_cast<double>(frame) /
                                 static_cast<double>(n_intervals);
            // Summed over this frame's steps under hydrodynamics, where every
            // CVodeReInit resets SUNDIALS' counters; a running total otherwise.
            SolverCounters frame_counters;
            if (!stepped)
            {
                // Through equal parts no longer than dt, so the wall-outflow
                // quadrature in advanceTo() has a step of at most dt however
                // far apart the frames are. Past the integrator's first call
                // CVODE only interpolates at each stop and steps on as it
                // would have, so the parts do not move the solution.
                const std::size_t n_parts =
                    channel->hasTransport()
                        ? std::max<std::size_t>(
                              1, static_cast<std::size_t>(
                                     std::ceil((t_rec - t_advanced) / opts.dt)))
                        : 1;
                const double t_frame_start = t_advanced; // [s]
                for (std::size_t part = 1; part <= n_parts; part++)
                {
                    advanceTo(part == n_parts
                                  ? t_rec
                                  : t_frame_start + (t_rec - t_frame_start) *
                                                        static_cast<double>(part) /
                                                        static_cast<double>(n_parts));
                }
                frame_counters = readSolverCounters();
            }
            else
            {
                // One global step at a time, NUMERICS.md section 4.2, hydro
                // outermost: half a hydro step, the circuit closed over the
                // step, the whole reaction step, half a hydro step. Conduction
                // and diffusion ride inside the reaction step, because they
                // live in PlasmaChannel1D::eval() and CVODE integrates them
                // together with the chemistry.
                for (std::size_t sub = 0; sub < steps_per_frame; sub++)
                {
                    // Derived from the step index rather than accumulated, so
                    // the step's end time cannot drift off the output grid
                    // over a long run.
                    const std::size_t step = (frame - 1) * steps_per_frame + sub;
                    const double t0 = static_cast<double>(step) * opts.dt;
                    const double t1 = static_cast<double>(step + 1) * opts.dt;
                    if (opts.cfg.hydrodynamics)
                    {
                        hydroHalfStep(*channel, hydro_cfg, x, channel_offset,
                                      n_vars, t0, 0.5 * opts.dt, history);
                    }
                    // After the half-step, so the resistance the circuit closes
                    // on is the expanded channel's. The first step is already
                    // closed, by the seeding the warm-up needed.
                    if (!(step == 0 && circuit_seeded))
                    {
                        circuit.closeOver(t0, t1, x);
                    }
                    if (opts.cfg.hydrodynamics)
                    {
                        // The hydro halves rewrite Tg/Te/Y and the density
                        // behind the integrator's back; getState() reads the
                        // live `x`, so reinitialize() picks up what they wrote.
                        for (SubSystem &system : systems)
                        {
                            system.cvode->reinitialize(t0, *system.rhs);
                        }
                        advanceTo(t1);
                        const SolverCounters step_counters = readSolverCounters();
                        frame_counters.steps += step_counters.steps;
                        frame_counters.rhs_evals += step_counters.rhs_evals;
                        frame_counters.jac_evals += step_counters.jac_evals;
                        frame_counters.err_test_fails += step_counters.err_test_fails;
                        frame_counters.nonlin_conv_fails +=
                            step_counters.nonlin_conv_fails;
                        frame_counters.lin_rhs_evals += step_counters.lin_rhs_evals;
                        frame_counters.lin_setups += step_counters.lin_setups;
                        frame_counters.last_order = step_counters.last_order;
                        hydroHalfStep(*channel, hydro_cfg, x, channel_offset,
                                      n_vars, t0 + 0.5 * opts.dt, 0.5 * opts.dt,
                                      history);
                    }
                    else
                    {
                        advanceTo(t1);
                    }
                }
                if (!opts.cfg.hydrodynamics)
                {
                    frame_counters = readSolverCounters();
                }
            }
            // A stepped run advances to (step+1)*dt, so that is the time the
            // recorded state is at; t_rec divides the total instead and lands
            // an ulp away. Labelling a frame with a time the state is not at
            // would put the recorded history off its own time grid.
            record(stepped
                       ? static_cast<double>(frame * steps_per_frame) * opts.dt
                       : t_rec);
            recordSolverCounters(frame_counters);
            recordCircuit();
            if (opts.loglevel > 0)
            {
                std::fprintf(stderr, "[bdf] frame %zu/%zu  t=%.3e s\n",
                             frame, n_intervals, t_rec);
                std::fflush(stderr);
            }
        }
        history.nframes = history.t.size();
        return history;
    }
 
} // namespace rizer