HydroStep.cpp

Go to the documentation of this file.

#include "models/nrp/HydroStep.h"
 
#include "numeric/EulerFlux.h"
#include "numeric/Format.h"
 
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <limits>
#include <stdexcept>
#include <string>
#include <vector>
 
namespace rizer
{
 
    namespace
    {
 
        // -------------------------------------------------------------------------- //
        //   File-local helpers                                                       //
        // -------------------------------------------------------------------------- //

        struct Workspace
        {
            std::size_t n_cells = 0; // N + 2
            std::size_t n_species = 0;
            std::vector<double> rho, u_r, Tg, Te, Y;
            std::vector<double> p_heavy, p_electron, p_total, sound_speed;
            std::vector<double> U_rho, U_mom, U_Y, U_total_energy, U_electron_energy;
 
            void resize(std::size_t cells, std::size_t species)
            {
                n_cells = cells;
                n_species = species;
                rho.assign(cells, 0.0);
                u_r.assign(cells, 0.0);
                Tg.assign(cells, 0.0);
                Te.assign(cells, 0.0);
                Y.assign(cells * species, 0.0);
                p_heavy.assign(cells, 0.0);
                p_electron.assign(cells, 0.0);
                p_total.assign(cells, 0.0);
                sound_speed.assign(cells, 0.0);
                U_rho.assign(cells, 0.0);
                U_mom.assign(cells, 0.0);
                U_Y.assign(cells * species, 0.0);
                U_total_energy.assign(cells, 0.0);
                U_electron_energy.assign(cells, 0.0);
            }
        };
 
        void requireSize(const std::vector<double> &v, std::size_t expected,
                         const char *name)
        {
            if (v.size() != expected)
            {
                throw std::invalid_argument(
                    std::string("hydroStep: ") + name + " has length " +
                    std::to_string(v.size()) + ", expected " +
                    std::to_string(expected) + ".");
            }
        }

        void requireFinite(const std::vector<double> &v, const char *name)
        {
            for (std::size_t i = 0; i < v.size(); i++)
            {
                if (!std::isfinite(v[i]))
                {
                    throw std::invalid_argument(
                        std::string("hydroStep: ") + name + "[" + std::to_string(i) +
                        "] is " + std::to_string(v[i]) + ", which is not finite.");
                }
            }
        }

        void requireFinite(double value, const char *name)
        {
            if (!std::isfinite(value))
            {
                throw std::invalid_argument(
                    std::string("hydroStep: ") + name + " is " +
                    std::to_string(value) + ", which is not finite.");
            }
        }

        void requireUnitSum(const std::vector<double> &Y, const char *name)
        {
            double sum = 0.0;
            for (const double Y_k : Y)
                sum += Y_k;
            // Loose enough to pass a normalized vector summed in any order, tight
            // enough that nothing a caller would consider "close" slips through.
            constexpr double tolerance = 1.0e-10;
            if (!(std::abs(sum - 1.0) <= tolerance))
            {
                throw std::invalid_argument(
                    std::string("hydroStep: ") + name +
                    " sums to " + std::to_string(sum) +
                    ", which is not 1 to within 1e-10.");
            }
        }
 
    } // namespace
 
    // ------------------------------------------------------------------------------ //
    //   Free functions                                                               //
    // ------------------------------------------------------------------------------ //
 
    HydroStepResult hydroStep(const HydroStepConfig &cfg, HydroStepState &state)
    {
        // Four steps, the third sub-cycled:
        //   1. validate every input, before touching the state
        //   2. widen the interior into ghost-extended working arrays
        //   3. sub-cycle to the CFL limit, each pass being the seven stages
        //      banner-marked below (ghosts, conservative, step size, boundary
        //      regime, faces, forward Euler, back to primitive)
        //   4. narrow the working arrays back onto the caller's state
 
        // ---- 1. validate the inputs -----------------------------------------
        // Check the inputs before doing any work, so a caller cannot get a half-done
        // state and have to guess which of the many possible preconditions failed.
        const std::size_t n_points = state.rho.size();
        const std::size_t n_species = state.n_species;
        if (n_points < 2)
        {
            throw std::invalid_argument(
                "hydroStep: needs at least 2 radial nodes, got " +
                std::to_string(n_points) +
                "; a single node has no face to carry a flux, and no node "
                "spacing to place the ghost cells one beyond each end.");
        }
        if (n_species == 0)
        {
            throw std::invalid_argument("hydroStep: the state has no species.");
        }
        requireSize(state.u_r, n_points, "state.u_r");
        requireSize(state.Tg, n_points, "state.Tg");
        requireSize(state.Te, n_points, "state.Te");
        requireSize(state.Y, n_points * n_species, "state.Y");
        requireSize(cfg.r_node, n_points, "cfg.r_node");
        requireSize(cfg.r_face, n_points + 1, "cfg.r_face");
        requireSize(cfg.V_cell, n_points, "cfg.V_cell");
        requireFinite(state.rho, "state.rho");
        requireFinite(state.u_r, "state.u_r");
        requireFinite(state.Tg, "state.Tg");
        requireFinite(state.Te, "state.Te");
        requireFinite(state.Y, "state.Y");
        requireFinite(cfg.r_node, "cfg.r_node");
        requireFinite(cfg.r_face, "cfg.r_face");
        requireFinite(cfg.V_cell, "cfg.V_cell");
        requireFinite(cfg.Y_amb, "cfg.Y_amb");
        requireFinite(cfg.p_amb, "cfg.p_amb");
        requireFinite(cfg.T_amb, "cfg.T_amb");
        requireFinite(cfg.t0, "cfg.t0");
        requireFinite(cfg.dt_total, "cfg.dt_total");
        requireFinite(cfg.cfl, "cfg.cfl");
        if (cfg.thermo == nullptr)
        {
            throw std::invalid_argument("hydroStep: cfg.thermo is null.");
        }
        if (cfg.dt_total < 0.0)
        {
            throw std::invalid_argument("hydroStep: dt_total must be >= 0.");
        }
        if (!(cfg.cfl > 0.0 && cfg.cfl <= 1.0))
        {
            throw std::invalid_argument(
                "hydroStep: cfl must be in (0, 1], got " + std::to_string(cfg.cfl));
        }
        if (cfg.electron_index < -1 || cfg.electron_index >= static_cast<int>(n_species))
        {
            throw std::invalid_argument(
                "hydroStep: electron_index must be -1 or a valid species index.");
        }
        for (std::size_t i = 0; i < n_points; i++)
        {
            if (!(cfg.V_cell[i] > 0.0))
            {
                throw std::invalid_argument(
                    "hydroStep: V_cell[" + std::to_string(i) + "] must be > 0.");
            }
            if (cfg.r_face[i + 1] < cfg.r_face[i])
            {
                throw std::invalid_argument(
                    "hydroStep: r_face must be non-decreasing, but r_face[" +
                    std::to_string(i + 1) + "] < r_face[" + std::to_string(i) +
                    "]. Equal face radii are the planar entry point, where every "
                    "face carries unit area.");
            }
            if (!(cfg.r_face[i] >= 0.0))
            {
                throw std::invalid_argument(
                    "hydroStep: r_face[" + std::to_string(i) + "] must be >= 0.");
            }
        }
        const bool inner_face_has_area = cfg.r_face[0] > 0.0;
        if (inner_face_has_area && !cfg.ghost_state)
        {
            throw std::invalid_argument(
                "hydroStep: r_face[0] > 0, so the inner face has real area and needs "
                "a boundary state, but no ghost_state callback was given. "
                "outer_bc_mode governs the outer face only.");
        }
        const bool needs_ambient = !cfg.ghost_state;
        if (needs_ambient)
        {
            // Both outer modes close the ghost density from p_amb, so both need it
            // and T_amb; AmbientGhost needs the composition as well.
            if (!(cfg.p_amb > 0.0) || !(cfg.T_amb > 0.0))
            {
                throw std::invalid_argument(
                    "hydroStep: outer_bc_mode needs p_amb > 0 and T_amb > 0, got "
                    "p_amb = " +
                    std::to_string(cfg.p_amb) + " Pa, T_amb = " +
                    std::to_string(cfg.T_amb) +
                    " K. A zero closes the ghost "
                    "density at zero, which is a vacuum wall rather than an error.");
            }
            if (cfg.outer_bc_mode == OuterBc::AmbientGhost)
            {
                requireSize(cfg.Y_amb, n_species, "cfg.Y_amb");
                requireUnitSum(cfg.Y_amb, "cfg.Y_amb");
            }
        }
 
        const bool has_electron = cfg.electron_index >= 0;
        const std::size_t electron =
            has_electron ? static_cast<std::size_t>(cfg.electron_index) : 0;
        const double cv_electron = has_electron ? cfg.thermo->electronCv() : 0.0;
 
        // Two different length scales, and they are not interchangeable.
        //
        // `face_span` is r_{i+1/2} - r_{i-1/2}, which is what appendix A.1's
        // geometric-source integral needs: A.3's well-balancedness cancels it
        // against the flux difference p*(r_east - r_west) EXACTLY, so it has to be
        // the same subtraction, not an algebraically equal rearrangement.
        //
        // `cell_width` is the physical extent the CFL bound is measured over,
        // V_i divided by the mean face area. In cylindrical mode that reduces to
        // the same r_{i+1/2} - r_{i-1/2}; in the planar entry point, where every
        // face radius is 1 and the face span is zero, it correctly returns V_i,
        // which is the cell width there. Using the face span for the CFL would
        // hand the planar mode a zero step.
        std::vector<double> face_span(n_points, 0.0);
        std::vector<double> cell_width(n_points, 0.0);
        for (std::size_t i = 0; i < n_points; i++)
        {
            face_span[i] = cfg.r_face[i + 1] - cfg.r_face[i];
            const double mean_face_area =
                0.5 * (cfg.r_face[i] + cfg.r_face[i + 1]);
            if (!(mean_face_area > 0.0))
            {
                throw std::invalid_argument(
                    "hydroStep: cell " + std::to_string(i) +
                    " has zero mean face area, so no CFL length can be derived.");
            }
            cell_width[i] = cfg.V_cell[i] / mean_face_area;
        }
 
        // ---- 2. widen into the ghost-extended arrays ------------------------
        // Create the working arrays, which hold the interior cells plus one ghost on
        // each side. The ghosts are in the same arrays so the face loop below has no
        // special case at either end.
        Workspace work;
        work.resize(n_points + 2, n_species);
        const std::size_t n_faces = n_points + 1;
        std::vector<double> flux_rho(n_faces, 0.0);
        std::vector<double> flux_mom(n_faces, 0.0);
        std::vector<double> flux_Y(n_faces * n_species, 0.0);
        std::vector<double> flux_total_energy(n_faces, 0.0);
        std::vector<double> flux_electron_energy(n_faces, 0.0);
        std::vector<double> velocity_flux(n_faces, 0.0);
 
        // Initialize the result.
        HydroStepResult result;
        result.boundary_species_mass_out.assign(n_species, 0.0);
        result.rho_min = std::numeric_limits<double>::infinity();
        result.p_min = std::numeric_limits<double>::infinity();
 
        // Copy the interior state into the working arrays, leaving the ghosts to be
        // filled in below. The caller's state is the interior only, so the ghosts
        // are not copied back to it at the end.
        for (std::size_t i = 0; i < n_points; i++)
        {
            const std::size_t c = i + 1;
            work.rho[c] = state.rho[i];
            work.u_r[c] = state.u_r[i];
            work.Tg[c] = state.Tg[i];
            work.Te[c] = state.Te[i];
            for (std::size_t k = 0; k < n_species; k++)
            {
                work.Y[c * n_species + k] = state.Y[i * n_species + k];
            }
        }
 
        // Radius of each extended cell, handed to the thermo accessor and to the
        // ghost-state callback: the NODE radius inside (not the cell midpoint,
        // which differs from it on a stretched grid and at the axis cell), and
        // for a ghost the adjacent interior node reflected across the boundary.
        //
        // A ghost given the FACE radius instead would carry its state -- and
        // under a prescribed equation of state its pressure -- half a cell
        // inside where the ghost cell is, which is an O(1) error in the boundary
        // cell's momentum divergence that refinement does not remove.
        //
        // The reflection is built from the node array rather than as
        // 2*r_face - r_node: in the planar entry point of NUMERICS.md section
        // 11.5 every face radius is 1, an area factor and not a position at all.
        // On a cylindrical grid the two agree, because r_face is mirrored half a
        // spacing beyond the end node.
        const double inner_spacing = cfg.r_node[1] - cfg.r_node[0];
        const double outer_spacing =
            cfg.r_node[n_points - 1] - cfg.r_node[n_points - 2];
        // The axis ghost mirrors the axis CELL onto itself: that cell's node sits
        // on the axis, so reflecting it across r = 0 leaves it there. Only a face
        // with real area has a ghost a full spacing beyond it.
        const double r_ghost_inner =
            inner_face_has_area ? cfg.r_node[0] - inner_spacing : cfg.r_node[0];
        const double r_ghost_outer = cfg.r_node[n_points - 1] + outer_spacing;
 
        // An annulus whose inner radius is smaller than its node spacing
        // reflects its first node past the axis. r_face is checked for sign
        // above, but this radius is derived from r_node, which is not, and it
        // is handed to the ghost_state callback and to the thermo accessor as
        // the POSITION to evaluate a prescribed field at. A field going as
        // 1/r, log r or sqrt(r) returns a sign-flipped or non-finite state
        // there and the step carries on, so the condition is refused here
        // rather than left to surface as a wrong boundary cell.
        //
        // Cylindrical geometry only. In the planar entry point of NUMERICS.md
        // section "Full system, planar - Sod shock tube" the coordinate is a position on a line and r_face is an
        // area factor of 1 everywhere, so the inner face always "has area",
        // the first node sits at x = 0, and a ghost at x = -dx is exactly
        // where it belongs rather than an error.
        if (cfg.cylindrical_source && inner_face_has_area && r_ghost_inner < 0.0)
        {
            throw std::invalid_argument(
                "hydroStep: the inner ghost cell lands at r = " +
                std::to_string(r_ghost_inner) +
                " m, past the axis. The annulus inner radius (" +
                std::to_string(cfg.r_node[0]) +
                " m) is smaller than its first node spacing (" +
                std::to_string(inner_spacing) +
                " m), so there is no radius a ghost one spacing inside the "
                "first node could occupy; use more nodes or a larger r_min.");
        }
 
        auto cellRadius = [&](std::size_t c)
        {
            if (c == 0)
                return r_ghost_inner;
            if (c == n_points + 1)
                return r_ghost_outer;
            return cfg.r_node[c - 1];
        };
 
        // No sum check here: this writes both CALLER-supplied ghosts and ghosts
        // the sub-step derives from its own interior cells. The sum belongs on
        // the first kind only -- see requireUnitSum's own comment, and the
        // explicit calls at the ghost_state call sites below. An interior-derived
        // ghost carries exactly what the interior carries, including the
        // round-off the scheme accumulates in sum_s Y_s over many sub-cycles, and
        // the interior itself is never held to an input tolerance. Checking the
        // copy but not the original turns the solver's own round-off into a hard
        // abort on fine grids, where the smaller CFL step means more sub-cycles
        // to accumulate it.
        auto writeGhost = [&](std::size_t c, const NodeState &ghost)
        {
            if (ghost.Y.size() != n_species)
            {
                throw std::invalid_argument(
                    "hydroStep: a ghost state carries " +
                    std::to_string(ghost.Y.size()) +
                    " mass fractions, expected " + std::to_string(n_species) + ".");
            }
            work.rho[c] = ghost.rho;
            work.u_r[c] = ghost.u_r;
            work.Tg[c] = ghost.Tg;
            work.Te[c] = ghost.Te;
            for (std::size_t k = 0; k < n_species; k++)
            {
                work.Y[c * n_species + k] = ghost.Y[k];
            }
        };
 
        auto cellAsNodeState = [&](std::size_t c)
        {
            NodeState node;
            node.rho = work.rho[c];
            node.u_r = work.u_r[c];
            node.Tg = work.Tg[c];
            node.Te = work.Te[c];
            const auto begin = work.Y.begin() + static_cast<std::ptrdiff_t>(c * n_species);
            node.Y.assign(begin, begin + static_cast<std::ptrdiff_t>(n_species));
            return node;
        };
 
        // ---- 3. sub-cycle to the CFL limit ----------------------------------
        // The caller's time is the start of the sub-step, and the total interval is
        // the Strang half-step. The sub-cycles below advance t in place, so
        // time_left is the remaining interval to reach the end of the sub-step.
        // The ambient ghost: the far field really is this state, so imposing all
        // of it is accurate rather than over-specified, and it anchors the domain
        // to ambient as the flow settles. It is the same state on every
        // sub-cycle, so it is built once.
        NodeState ambient_ghost;
        if (!cfg.ghost_state && cfg.outer_bc_mode == OuterBc::AmbientGhost)
        {
            ambient_ghost.u_r = 0.0;
            ambient_ghost.Tg = cfg.T_amb;
            ambient_ghost.Te = cfg.T_amb;
            ambient_ghost.Y = cfg.Y_amb;
            ambient_ghost.rho = cfg.thermo->densityFromPressure(
                cfg.p_amb, cfg.T_amb, cfg.T_amb, ambient_ghost.Y.data());
        }
 
        double t = cfg.t0;
        double time_left = cfg.dt_total;
        const std::size_t outer = n_points + 1;
 
        // One pass per sub-cycle, plus a pass that stops after the step-size
        // stage when dt_total is zero: that probe then reports the dt_cfl the
        // first sub-cycle would apply, from the same ghosts and face speeds,
        // rather than from a second copy of the bound.
        for (;;)
        {
            // ---- boundary ghosts -------------------------------------------
            if (cfg.ghost_state)
            {
                // Caller-supplied, so the sum is held to the input tolerance here.
                auto checkGhost = [&](const NodeState &ghost)
                {
                    requireSize(ghost.Y, n_species, "a ghost_state callback's Y");
                    requireUnitSum(ghost.Y, "a ghost_state callback's Y");
                    // The sum above already rejects a NaN in Y, since a NaN sum
                    // equals nothing. The other four slots have no such check and
                    // are the one place arbitrary caller code reaches the inner
                    // loop, so they get an explicit one.
                    requireFinite(ghost.rho, "a ghost_state callback's rho");
                    requireFinite(ghost.u_r, "a ghost_state callback's u_r");
                    requireFinite(ghost.Tg, "a ghost_state callback's Tg");
                    requireFinite(ghost.Te, "a ghost_state callback's Te");
                };
                if (inner_face_has_area)
                {
                    // For the inner face, the caller's callback is the only source of a
                    // ghost state. Only used for a uniform grid with r_face[0] > 0,
                    // which is a planar entry point with a real boundary, so the caller
                    // has to supply a state. The sub-step does not know what the
                    // physics is there, so it cannot derive a ghost from the interior.
                    const NodeState inner =
                        cfg.ghost_state(r_ghost_inner, t, true, cellAsNodeState(1));
                    checkGhost(inner);
                    writeGhost(0, inner);
                }
                const NodeState outer_ghost = cfg.ghost_state(
                    r_ghost_outer, t, false, cellAsNodeState(n_points));
                checkGhost(outer_ghost);
                writeGhost(outer, outer_ghost);
            }
            else if (cfg.outer_bc_mode == OuterBc::AmbientGhost)
            {
                writeGhost(outer, ambient_ghost);
            }
            else
            {
                // PressureOutlet: exactly ONE quantity imposed, p = p_amb;
                // u_r, Tg, Te and Y_s extrapolated with zero gradient from the
                // adjacent interior cell; rho closed from the equation of state at
                // that p_amb. Extrapolating rho as well would already fix p at the
                // face and make p = p_amb one condition too many. Getting that
                // count right is still not enough to stop reflections: pinning p
                // is a pressure-release surface, so this mode reflects far more
                // than AmbientGhost does (NUMERICS.md section "Outgoing acoustic pulse, the boundary test", case 21).
                NodeState ghost = cellAsNodeState(n_points);
                ghost.rho = cfg.thermo->densityFromPressure(
                    cfg.p_amb, ghost.Tg, ghost.Te, ghost.Y.data());
                writeGhost(outer, ghost);
            }
            if (!inner_face_has_area)
            {
                // The axis face has zero area, so nothing is ever read through it.
                // The mirror only keeps the flux expression finite.
                NodeState axis = cellAsNodeState(1);
                axis.u_r = -axis.u_r;
                writeGhost(0, axis);
            }
 
            // ---- primitive -> conservative, plus the face-blend inputs -------
            for (std::size_t c = 0; c < work.n_cells; c++)
            {
                const double *Y_c = &work.Y[c * n_species];
                const double r_c = cellRadius(c);
                // One call for all four state functions, so an accessor with a
                // costly state-setting (CanteraThermo) sets it once per cell.
                const StateProperties properties = cfg.thermo->stateProperties(
                    r_c, t, work.rho[c], work.Tg[c], work.Te[c], Y_c);
                work.p_heavy[c] = properties.p_heavy;
                work.p_electron[c] = properties.p_electron;
                work.p_total[c] = work.p_heavy[c] + work.p_electron[c];
                // Each later pass starts from the state the previous sub-cycle
                // produced, so its interior pressures are that sub-cycle's
                // minimum candidates; the final state is scanned after the loop.
                if (result.n_subcycles_used > 0 && c >= 1 && c <= n_points)
                {
                    result.p_min = std::min(result.p_min, work.p_total[c]);
                }
                work.sound_speed[c] = properties.sound_speed;
 
                // Construct U=[ρ, ρu, ρY, ρE, ρe_e] from the primitive state.
                work.U_rho[c] = work.rho[c];
                work.U_mom[c] = work.rho[c] * work.u_r[c];
                for (std::size_t k = 0; k < n_species; k++)
                {
                    work.U_Y[c * n_species + k] = work.rho[c] * Y_c[k];
                }
                // One conservative TOTAL energy row and one sourced ELECTRON
                // internal-energy row; the heavy pool has no row of its own and
                // is recovered below as total minus kinetic minus electron, which
                // hands it -p_H div(u) exactly (NUMERICS.md section 5.4).
                // Splitting the shared kinetic term BETWEEN two conservative pool
                // rows instead leaves each one a spurious
                // Y_pool u.grad(p) - u.grad(p_pool); no weight removes it, because
                // the one that would is grad(p_pool)/grad(p) and varies in space.
                const double Y_electron = has_electron ? Y_c[electron] : 0.0;
                const double rho_electron = work.rho[c] * Y_electron;
                work.U_electron_energy[c] = rho_electron * cv_electron * work.Te[c];
                work.U_total_energy[c] =
                    work.rho[c] * properties.heavy_energy +
                    work.U_electron_energy[c] +
                    0.5 * work.rho[c] * work.u_r[c] * work.u_r[c];
            }
 
            // ---- step size ---------------------------------------------------
            // Section 3.3 requires the local CFL to be computed every step, so the
            // bound is re-derived here rather than fixed once from the initial
            // state: a rarefaction accelerates the flow within a single half-step.
            //
            // The speed that bounds a cell is the one its two FACES apply, which
            // is maxWaveSpeed over the cell and each neighbour -- the ghosts
            // included, since the boundary faces blend against them. Scanning the
            // interior cells alone would derive the bound from a speed smaller
            // than the dissipation actually uses at those faces whenever a ghost
            // is faster than every interior cell, which an ambient hotter than the
            // channel or a ghost_state callback imposing an inflow both produce.
            double dt_cfl = std::numeric_limits<double>::infinity();
            for (std::size_t i = 0; i < n_points; i++)
            {
                const std::size_t c = i + 1;
                const double speed = std::max(
                    numeric::maxWaveSpeed(work.u_r[c - 1], work.sound_speed[c - 1],
                                          work.u_r[c], work.sound_speed[c]),
                    numeric::maxWaveSpeed(work.u_r[c], work.sound_speed[c],
                                          work.u_r[c + 1], work.sound_speed[c + 1]));
                if (speed > 0.0)
                {
                    dt_cfl = std::min(dt_cfl, cell_width[i] / speed);
                }
            }
            if (result.n_subcycles_used == 0)
            {
                result.dt_cfl = dt_cfl;
            }
            if (!(cfg.dt_total > 0.0))
            {
                // Zero-length probe: nothing to advance, but the caller still
                // gets the positivity minima of the state it handed in.
                for (std::size_t c = 1; c <= n_points; c++)
                {
                    result.rho_min = std::min(result.rho_min, work.rho[c]);
                    result.p_min = std::min(result.p_min, work.p_total[c]);
                }
                break;
            }
 
            double dt = 0.0;
            if (cfg.n_subcycles > 0)
            {
                // Equal fixed sub-steps, a test-only control for the cases where a
                // known step size matters more than the stability margin.
                dt = cfg.dt_total / static_cast<double>(cfg.n_subcycles);
                if (cfg.refuse_above_cfl && dt > dt_cfl)
                {
                    throw std::invalid_argument(
                        "hydroStep: the pinned sub-step " + significant(dt) +
                        " s exceeds the CFL bound " + significant(dt_cfl) +
                        " s at this state, where the explicit update is unstable. "
                        "Pin a smaller hydro_substep.");
                }
            }
            else
            {
                if (!std::isfinite(dt_cfl))
                {
                    throw std::invalid_argument(
                        "hydroStep: every wave speed is zero, so no CFL step can be "
                        "derived; pass n_subcycles explicitly.");
                }
                dt = cfg.cfl * dt_cfl;
            }
            if (dt > time_left)
                dt = time_left;
            if (!(dt > 0.0))
            {
                throw std::invalid_argument(
                    "hydroStep: the derived sub-step is not positive.");
            }
 
            // ---- the boundary regime the characteristic condition assumes -----
            if (!cfg.ghost_state && cfg.outer_bc_mode == OuterBc::PressureOutlet)
            {
                const double u_wall = work.u_r[n_points];
                const double c_wall = work.sound_speed[n_points];
                // The guard is for a RESOLVED regime change, so the inflow half is
                // thresholded on a Mach number rather than on the sign. A boundary
                // cell at rest produces a u_r of order 1e-16 m/s from round-off
                // alone, and the wake behind an outgoing acoustic wave leaves a few
                // times 1e-5 of the sound speed pointing inward; neither is the
                // subsonic-inflow regime, which needs almost the whole state
                // supplied rather than a single imposed pressure. Both supersonic
                // regimes are rejected outright, with no tolerance.
                constexpr double inflow_mach_threshold = 1.0e-3;
                if (u_wall < -inflow_mach_threshold * c_wall || u_wall >= c_wall)
                {
                    throw std::invalid_argument(
                        "hydroStep: outer_bc_mode=pressure_outlet assumes SUBSONIC "
                        "OUTFLOW (0 <= u_r < c) at R_max, but u_r = " +
                        std::to_string(u_wall) + " m/s against c = " +
                        std::to_string(c_wall) +
                        " m/s. Inflow, supersonic outflow and supersonic inflow "
                        "each change the characteristic count in a different "
                        "direction and need their own treatment.");
                }
            }
 
            // ---- faces --------------------------------------------------------
            for (std::size_t f = 0; f < n_faces; f++)
            {
                const std::size_t left = f;
                const std::size_t right = f + 1;
                const double r_f = cfg.r_face[f];
                const double wave = numeric::maxWaveSpeed(
                    work.u_r[left], work.sound_speed[left], work.u_r[right],
                    work.sound_speed[right]);
 
                // Construct F=[., ., ., ., .] using the Rusanov fluxes, which are the
                // same for every pool and every species.
                flux_rho[f] = numeric::rusanovScalarFlux(
                    r_f, work.U_rho[left], work.U_rho[right], work.u_r[left],
                    work.u_r[right], wave);
                flux_mom[f] = numeric::rusanovMomentumFlux(
                    r_f, work.rho[left], work.u_r[left], work.p_total[left],
                    work.rho[right], work.u_r[right], work.p_total[right], wave);
                for (std::size_t k = 0; k < n_species; k++)
                {
                    flux_Y[f * n_species + k] = numeric::rusanovScalarFlux(
                        r_f, work.U_Y[left * n_species + k],
                        work.U_Y[right * n_species + k], work.u_r[left],
                        work.u_r[right], wave);
                }
                flux_total_energy[f] = numeric::rusanovEnergyFlux(
                    r_f, work.U_total_energy[left], work.p_total[left],
                    work.u_r[left], work.U_total_energy[right],
                    work.p_total[right], work.u_r[right], wave);
                flux_electron_energy[f] = numeric::rusanovScalarFlux(
                    r_f, work.U_electron_energy[left],
                    work.U_electron_energy[right], work.u_r[left],
                    work.u_r[right], wave);
                velocity_flux[f] =
                    numeric::faceVelocityFlux(r_f, work.u_r[left], work.u_r[right]);
            }
 
            // ---- forward Euler, NUMERICS.md (10), (11), (13), (14), (15) --------
            // Once the fluxes are computed, the interior cells are updated with a single
            // forward-Euler step. The ghosts are never updated, so the fluxes at the
            // domain boundaries are always computed from the same ghost state.
            for (std::size_t i = 0; i < n_points; i++)
            {
                const std::size_t c = i + 1;
                const double inv_volume = 1.0 / cfg.V_cell[i];
 
                work.U_rho[c] -= dt * (flux_rho[i + 1] - flux_rho[i]) * inv_volume;
 
                double momentum_rate = -(flux_mom[i + 1] - flux_mom[i]) * inv_volume;
                if (cfg.cylindrical_source)
                {
                    momentum_rate +=
                        numeric::momentumGeometricSourceIntegral(work.p_total[c],
                                                                 face_span[i]) *
                        inv_volume;
                }
                work.U_mom[c] += dt * momentum_rate;
 
                for (std::size_t k = 0; k < n_species; k++)
                {
                    work.U_Y[c * n_species + k] -=
                        dt * (flux_Y[(i + 1) * n_species + k] - flux_Y[i * n_species + k]) * inv_volume;
                }
 
                // The total row is a pure flux difference: the compression work
                // rides inside the enthalpy rusanovEnergyFlux() transports, so
                // every interior face appears twice with opposite signs and
                // cancels, and total energy changes only by what crosses the two
                // domain boundaries.
                work.U_total_energy[c] -=
                    dt * (flux_total_energy[i + 1] - flux_total_energy[i]) * inv_volume;
 
                // The electron row keeps its physical -p_e div(u) source and is
                // NOT conservative, deliberately: only the total carries a jump
                // condition, while how a shock partitions energy between the pools
                // is a modelling statement rather than a Rankine-Hugoniot one. The
                // divergence is differenced over the VELOCITY flux, never the mass
                // flux (appendix A.2).
                const double divergence =
                    (velocity_flux[i + 1] - velocity_flux[i]) * inv_volume;
                work.U_electron_energy[c] +=
                    dt * (-(flux_electron_energy[i + 1] - flux_electron_energy[i]) * inv_volume -
                          work.p_electron[c] * divergence);
                result.pressure_work_electron +=
                    dt * work.p_electron[c] * divergence * cfg.V_cell[i];
            }
 
            // Mass that crossed the two domain-boundary faces, so section 10's
            // identity can be checked as "domain mass change equals accumulated
            // boundary flux" rather than the weaker "mass constant".
            result.boundary_mass_out += dt * (flux_rho[n_points] - flux_rho[0]);
            for (std::size_t k = 0; k < n_species; k++)
            {
                result.boundary_species_mass_out[k] +=
                    dt * (flux_Y[n_points * n_species + k] - flux_Y[k]);
            }
            result.boundary_total_energy_out +=
                dt * (flux_total_energy[n_points] - flux_total_energy[0]);
            result.boundary_electron_energy_out +=
                dt * (flux_electron_energy[n_points] - flux_electron_energy[0]);
 
            // ---- conservative -> primitive ------------------------------------
            for (std::size_t i = 0; i < n_points; i++)
            {
                const std::size_t c = i + 1;
                // The density the scheme produced, and the only one: a
                // non-positive value is not repaired, because the kinetic energy
                // the total-energy row carries is 0.5*(rho u)^2/rho, which a
                // non-positive mass row leaves undefined.
                const double rho = work.U_rho[c];
                if (!(rho > 0.0))
                {
                    throw std::runtime_error(
                        "hydroStep: density went non-positive at node " +
                        std::to_string(i) +
                        " (rho = " +
                        std::to_string(rho) +
                        " kg/m^3). Lower cfl or refine the grid.");
                }
                // The kinetic energy the total-energy row actually carries. A
                // collapsed cell drives internal_heavy negative and
                // heavyTemperature() raises, which is the honest outcome.
                const double carried_kinetic = 0.5 * work.U_mom[c] * work.U_mom[c] / rho;
                result.rho_min = std::min(result.rho_min, rho);
 
                double *Y_c = &work.Y[c * n_species];
                for (std::size_t k = 0; k < n_species; k++)
                {
                    // No renormalization. The species fluxes sum to the mass flux
                    // exactly, so sum_s Y_s stays 1 to round-off on its own;
                    // rescaling here would hide a broken flux rather than reveal it.
                    Y_c[k] = work.U_Y[c * n_species + k] / rho;
                }
 
                // The electron row holds internal energy directly, so it inverts
                // with no kinetic bookkeeping of its own.
                const double Y_electron = has_electron ? Y_c[electron] : 0.0;
                const double rho_electron = rho * Y_electron;
                if (has_electron && !(rho_electron > 0.0))
                {
                    // Te is 0/0 here (or, for rho_e < 0, has no meaning), and no
                    // value carried forward would mean anything: electrons
                    // chemistry later creates in this cell would start at it. In
                    // the native solver a cell only gets here when it and its
                    // neighbours hold no electrons, which a seeded start rules
                    // out, or when a reaction sub-step hands in rho_e < 0: neither
                    // CVODE here nor the Python split's scipy vode enforces
                    // Y >= 0.
                    throw std::runtime_error(
                        "hydroStep: node " + std::to_string(i) +
                        " has a non-positive electron density (rho_e = " +
                        std::to_string(rho_electron) +
                        " kg/m^3), so its electron temperature is undefined. "
                        "Either the node started without electrons (seed a positive "
                        "electron fraction at every node) or the reaction sub-step "
                        "returned a negative electron mass fraction.");
                }
                if (has_electron)
                {
                    // Divide by rho_e, NOT rho: e_e and c_v,e are both per kg of
                    // ELECTRONS. Dividing by rho would undershoot Te by a factor
                    // Y_e, around 1e-10 in a weakly ionized channel. Both rows
                    // scale with rho_e -- the only source on the electron row is
                    // the compression work, proportional to p_e -- so no
                    // threshold is applied (NUMERICS.md section 9).
                    work.Te[c] =
                        work.U_electron_energy[c] / (rho_electron * cv_electron);
                }
                // The heavy pool has no row of its own: it is what the total holds
                // once the kinetic energy and the electron pool are taken out,
                // which gives it -p_H div(u) exactly. The momentum row's O(dr^2)
                // geometric-source quadrature lands here too, as it does in any
                // conservative Euler scheme recovering e = E - K.
                const double internal_heavy =
                    work.U_total_energy[c] - carried_kinetic -
                    rho_electron * cv_electron * work.Te[c];
 
                work.Tg[c] = cfg.thermo->heavyTemperature(
                    rho, Y_c, internal_heavy / rho, work.Tg[c]);
 
                work.rho[c] = rho;
                work.u_r[c] = work.U_mom[c] / rho;
 
                // The state is primitive again, which is where a NaN born inside
                // the step first becomes nameable: this is the only point that has
                // the node and the sub-cycle index to report it with. Y is left
                // out deliberately -- it is the expensive half at fifty-odd
                // species, and a NaN in one mass fraction reaches the thermo call
                // above within this same iteration, so it surfaces through Tg.
                if (!std::isfinite(rho) || !std::isfinite(work.u_r[c]) ||
                    !std::isfinite(work.Tg[c]) || !std::isfinite(work.Te[c]))
                {
                    throw std::runtime_error(
                        "hydroStep: node " + std::to_string(i) +
                        " went non-finite on sub-cycle " +
                        std::to_string(result.n_subcycles_used) + " (rho = " +
                        std::to_string(rho) + " kg/m^3, u_r = " +
                        std::to_string(work.u_r[c]) + " m/s, Tg = " +
                        std::to_string(work.Tg[c]) + " K, Te = " +
                        std::to_string(work.Te[c]) + " K).");
                }
            }
 
            result.n_subcycles_used++;
            t += dt;
            time_left -= dt;
            // A pinned count ends on the count, not on the accumulated time:
            // M subtractions of dt_total/M leave a rounding residual of order
            // M * eps * dt_total, which for a large M would otherwise exceed the
            // cut-off below and run one extra, round-off-sized sub-cycle.
            if (time_left < 1.0e-12 * cfg.dt_total ||
                (cfg.n_subcycles > 0 && result.n_subcycles_used == cfg.n_subcycles))
            {
                time_left = 0.0;
            }
            if (!(time_left > 0.0))
            {
                break;
            }
        }
 
        if (result.n_subcycles_used > 0)
        {
            // The final state's pressures: every earlier one was scanned at the
            // start of the sub-cycle that advanced from it.
            for (std::size_t i = 0; i < n_points; i++)
            {
                const std::size_t c = i + 1;
                double p_heavy = 0.0;
                double p_electron = 0.0;
                cfg.thermo->pressures(cfg.r_node[i], t, work.rho[c], work.Tg[c],
                                      work.Te[c], &work.Y[c * n_species], p_heavy,
                                      p_electron);
                result.p_min = std::min(result.p_min, p_heavy + p_electron);
            }
        }
 
        // ---- 4. narrow back onto the caller's state ------------------------
        for (std::size_t i = 0; i < n_points; i++)
        {
            const std::size_t c = i + 1;
            state.rho[i] = work.rho[c];
            state.u_r[i] = work.u_r[c];
            state.Tg[i] = work.Tg[c];
            state.Te[i] = work.Te[c];
            for (std::size_t k = 0; k < n_species; k++)
            {
                state.Y[i * n_species + k] = work.Y[c * n_species + k];
            }
        }
        return result;
    }
 
} // namespace rizer