rizer.models.nrp.plasma_channel_split_2T1D#

Two transient 1D radial plasma-channel solvers over one native C++ physics core.

Head-to-head comparison (the difference between them is numerics, not physics):

  • ChannelReactor1D.solve_monolithic() — method-of-lines: the whole coupled state [Tg, Te, Y] over all nodes is integrated by a single stiff CVODE/vode (the per-node reaction RHS comes from the C++ Reactor0D; radial conduction couples neighbours). No splitting error; one big stiff system.

  • ChannelReactor1D.solve_split() — operator splitting (Strang): each global step does ½ transport → an independent per-node stiff reaction solve → ½ transport. The OpenSMOKE++/reactPlasFoam pattern; each node is just a 0D reactor.

Both reuse isomass_2T_volume_reactor_cpp (the C++ Reactor0D with composition-resolved sigma = n_e e^2/(m_e nu_eH)), so a difference between the two is numerics, not physics. Constant volume per node (rho fixed), gas heat conduction only (electron conduction is intentionally omitted — a constant kappa_e over the vanishing electron heat capacity at the cold edge blows up), species frozen. The field E(t) is prescribed.

When hydrodynamics is enabled, the transport half-step solves the radially symmetric Euler equations, i.e. the planar ones plus a geometric source $-(alpha/r)$ multiplying $(rho u, rho u^2, (E+p)u)$ with $alpha = 1$ in cylindrical geometry ([LeVeque2002], equation (18.53); [Toro2009], equations (17.1)-(17.2)). The C++ sub-step carries that system in area-weighted form, with r-weighted fluxes over the cell volume and a p geometric source. NUMERICS.md section “The planar limit of a cylindrical release” checks the two forms agree by taking the planar limit, where an exact solution exists, and section “The cylindrical scheme against outside references” checks the cylindrical solution itself against [Toro2009]’s published explosion test and against an independent Lagrangian solver that carries no geometric source term at all.

Classes#

ChannelReactor1D

Radial 1D channel built from per-node C++ 0D reactors + gas conduction.

Module Contents#

class rizer.models.nrp.plasma_channel_split_2T1D.ChannelReactor1D(mech: str, phase: str, grid: numpy.ndarray | list[float], rho: float, mtcf: list[rizer.transport.collision_frequency.MomentumTransferCollisionFrequencyModel], species_names: list[str], reacting: bool = True, kappa: float = 0.0, T_amb: float = 300.0, p_ext: float = ct.one_atm, Te_n: int = 400, hydrodynamics: bool = False, Y_amb: numpy.ndarray | None = None, hydro_substep: float | None = None, outer_bc_mode: str = 'ambient_ghost')#

Radial 1D channel built from per-node C++ 0D reactors + gas conduction.

Parameters:

  • mech (str) – Mechanism (YAML) file path and phase name (a PlasmaPhase with e-).

  • phase (str) – Mechanism (YAML) file path and phase name (a PlasmaPhase with e-).

  • grid (array_like) – Radial node positions [m], increasing from 0 (axis) to the wall.

  • rho (float or array_like) – Mass density [kg/m^3], scalar or one value per node. A scalar is broadcast: it is the uniform field of the u_r = 0 limit (PHYSICS.md section 5.3), not a separate constant-volume model.

  • mtcf (list of MomentumTransferCollisionFrequencyModel) – Per-species momentum-transfer collision-frequency models (as built by get_momentum_transfer_collision_frequencies_list); drives the C++ Reactor0D collision model (see rizer.models.nrp.isomass_2T_volume_reactor_cpp.build_reactor0d()).

  • species_names (list of str) – Species names, in mechanism order.

  • reacting (bool, optional) – Include finite-rate chemistry in the per-node C++ RHS. Default True.

  • kappa (float, optional) – Gas thermal conductivity [W/m/K] for radial conduction (0 disables it). Default 0.0.

  • T_amb (float, optional) – Wall temperature [K] (Dirichlet boundary condition). Default 300.0.

  • p_ext (float, optional) – Ambient pressure [Pa] (only used by the reactor’s V-equation; constant-V here so it only sets the state used to evaluate cv). Default cantera.one_atm.

  • Te_n (int, optional) – Number of points in the C++ collision model’s electron-temperature table. Default 400.

  • hydrodynamics (bool, optional) – Add an explicit hydro half-step to the Strang loop, solving continuity, radial momentum, species advection and the two pool energies (NUMERICS.md section 5). Default False, which leaves the loop bit-identical to the pre-hydrodynamic sequence.

  • Y_amb (numpy.ndarray, optional) – Ambient composition at the outer boundary [-], length nsp; required when hydrodynamics is set. Must sum to 1, or the species fluxes stop telescoping into the mass flux at the wall.

  • hydro_substep (float, optional) – Fixed hydro sub-step size [s]. The default None lets the sub-step follow the CFL bound, which is what a production run wants. Pinning it is what NUMERICS.md section “Strang splitting order” case 18 needs: the sub-cycle count is a ceiling, so across a sweep of output intervals the CFL-chosen sub-step wanders over half a factor and the hydro error floor with it, which is a sawtooth the size of the splitting error being measured. Each half-step must then be an exact integer multiple of this value or the solve raises.

reactor#
names#
nsp#
nv#
r#
npts#
rho#
u_r#
T_amb#
p_ext#
hydrodynamics = False#
boundary_mass_out = 0.0#
boundary_total_energy_out = 0.0#
boundary_electron_energy_out = 0.0#
pressure_work_electron = 0.0#
face_radii() → numpy.ndarray#

Face radii \(r_{i\pm1/2}\), length npts+1 [m].

NUMERICS.md section 2’s conventions: faces at node midpoints, the axis face pinned to zero so it carries no area, and the wall face mirrored, \(r_{N-1/2} = r_{N-1} + (r_{N-1}-r_{N-2})/2\).

Returns:

Face radii [m]. The stored array, built once at construction from a mesh that never moves; callers must not write to it.

Return type:

numpy.ndarray, shape (npts+1,)

cell_volumes() → numpy.ndarray#

Cell cross-sections \(V_i=(r_{i+1/2}^2-r_{i-1/2}^2)/2\) [m^2].

Per unit axial length; the \(2\pi\) of the azimuthal integration cancels throughout and is dropped (NUMERICS.md section 2).

Returns:

Cell volumes [m^2]. The stored array, built once at construction alongside the faces; callers must not write to it.

Return type:

numpy.ndarray, shape (npts,)

cross_section_integral(field: numpy.ndarray) → float#

Integrate a radial field over the channel cross-section.

Evaluates \(\int f \, \mathrm{d}A = 2\pi\sum_i f_i V_i\) on the same cells the conduction operator and the hydro half-step use, so a radially-integrated quantity is not computed on a second grid.

A node trapezoid over the same profile is not equivalent. It weights each node by \(2\pi r\), which vanishes on the axis, so the axis node – the hottest and most conductive point in a channel – carries no weight at all, while the cell sum gives it the area \(\pi r_{1/2}^2\). On the reference discharge of tests/models/nrp/test_channel_circuit.py the two differ by 1.6 to 1.8%, the cell sum giving the larger conductance.

Parameters:

field (numpy.ndarray, shape (npts,)) – Per-node values of the quantity to integrate.

Returns:

The integral over the cross-section, in the field’s units times m^2.

Return type:

float

pressure_profile(Tg: numpy.ndarray, Te: numpy.ndarray, Y: numpy.ndarray) → numpy.ndarray#

Total pressure per node from the local state [Pa].

A diagnostic for the recorded history, not an input to anything: the reaction sub-step already sees the local pressure because it is handed the node’s density and derives p from (Tg, Te, Y, rho) itself. Once rho is solved, p_ext is the local pressure only in the initial uniform state, so a history that reported p_ext would be wrong everywhere the density profile has developed.

Parameters:

  • Tg (numpy.ndarray, shape (npts,)) – Heavy-pool and electron temperature per node [K].

  • Te (numpy.ndarray, shape (npts,)) – Heavy-pool and electron temperature per node [K].

  • Y (numpy.ndarray, shape (npts, nsp)) – Mass fractions per node [-].

Returns:

Total pressure per node [Pa].

Return type:

numpy.ndarray, shape (npts,)

conductivity_profile(Tg: numpy.ndarray, Te: numpy.ndarray, Y: numpy.ndarray) → numpy.ndarray#

Radial electrical-conductivity profile \(\sigma(r)\) [S/m].

Evaluates the composition-resolved Mitchner conductivity \(\sigma = n_e e^2/(m_e \bar\nu_{eH})\) (C++ collision model) at every node. The total current and plasma resistance follow as

\[I = E \int_0^{R} \sigma(r)\,2\pi r\,dr, \qquad R_p = \frac{\text{gap}}{\int_0^{R} \sigma(r)\,2\pi r\,dr}.\]
Parameters:

  • Tg (numpy.ndarray, shape (npts,)) – Gas and electron temperature per node [K].

  • Te (numpy.ndarray, shape (npts,)) – Gas and electron temperature per node [K].

  • Y (numpy.ndarray, shape (npts, nsp)) – Mass fractions per node.

Returns:

Conductivity per node [S/m].

Return type:

numpy.ndarray, shape (npts,)

solve_monolithic(t_out: numpy.ndarray, E_func: Callable[[float], float], Tg0: numpy.ndarray, Te0: numpy.ndarray, Y0: numpy.ndarray, rtol: float = 1e-07, atol: float = 1e-13) → dict[str, Any]#

Method-of-lines solve: one stiff integrator over the full coupled state.

Flattens the radial field into \(y = [T_g, T_e, Y]_j\) for \(j = 0\dots N-1\) and integrates

\[\frac{dy}{dt} = \underbrace{f_\text{react}(y_j, E(t))}_{\text{per node}} + \underbrace{\alpha\,(L\,T_g + w)}_{\text{radial conduction on }T_g}\]

with a single scipy vode/BDF over all nodes (radial conduction couples neighbours; electron conduction omitted). No splitting error; the reference against which the operator split is checked.

Parameters:

  • t_out (numpy.ndarray) – Output times [s] (t_out[0] is the initial time).

  • E_func (typing.Callable) – E_func(t) -> float, the prescribed field [V/m].

  • Tg0 (numpy.ndarray, shape (npts,)) – Initial gas / electron temperatures [K].

  • Te0 (numpy.ndarray, shape (npts,)) – Initial gas / electron temperatures [K].

  • Y0 (numpy.ndarray, shape (npts, nsp)) – Initial mass fractions.

  • rtol (float) – Integrator tolerances.

  • atol (float) – Integrator tolerances.

Returns:

History with keys t, Tg, Te, Y, ne (lists over frames).

Return type:

dict

solve_split(t_out: numpy.ndarray, E_func: Callable[[float], float], Tg0: numpy.ndarray, Te0: numpy.ndarray, Y0: numpy.ndarray, rtol: float = 1e-07, atol: float = 1e-13) → dict[str, Any]#

Strang operator-split solve (prescribed field).

Each global step \([t_n, t_{n+1}]\) of size \(\Delta t\) applies a nested Strang sequence, hydrodynamics outermost:

\[\tfrac12\Delta t\ \text{hydro},\quad \tfrac12\Delta t\ \text{cond},\quad \Delta t\ \text{react (per node)},\quad \tfrac12\Delta t\ \text{cond},\quad \tfrac12\Delta t\ \text{hydro} ,\]

i.e. a half hydro step (_hydro_half_step, a no-op unless hydrodynamics is set), a half conduction step (_transport_step), an independent per-node stiff reaction solve (each node a 0D reactor, _react_node), then the two half steps again in reverse. The external field \(E(t)\) is evaluated continuously inside the reaction substep (it is forcing, not a split operator), so the split carries no field-sampling error and matches solve_monolithic() to integrator tolerance.

Second-order in \(\Delta t\) only where the reaction side is not stiff and conduction is negligible: the self-difference falls by 3.994 per halving with the field off and kappa inert (NUMERICS.md section “Strang splitting order”, case 18). Under Joule heating the measured rate drops to between 1.84 and 2.57, the order reduction splitting suffers when one operator carries a transient faster than the step. With conduction moving the solution it drops to 1.95, because _transport_step is one backward-Euler solve, first order in a step size that scales with \(\Delta t\), so its error does not cancel the way the pinned hydro sub-step’s does.

Parameters:

  • t_out (numpy.ndarray) – Output / split times [s]; the step size is set by consecutive entries.

  • E_func (typing.Callable) – E_func(t) -> float field [V/m].

  • Tg0 (numpy.ndarray, shape (npts,)) – Initial temperatures [K].

  • Te0 (numpy.ndarray, shape (npts,)) – Initial temperatures [K].

  • Y0 (numpy.ndarray, shape (npts, nsp)) – Initial mass fractions.

  • rtol (float) – Per-node integrator tolerances.

  • atol (float) – Per-node integrator tolerances.

Returns:

History with keys t, Tg, Te, Y, ne.

Return type:

dict

solve_split_circuit(t_out: numpy.ndarray, circuit: rizer.electrical_model.circuit.transmission_line_circuit.TransmissionLineResistiveCircuit, gap: float, Tg0: numpy.ndarray, Te0: numpy.ndarray, Y0: numpy.ndarray, rtol: float = 1e-07, atol: float = 1e-13) → dict[str, Any]#

Operator-split solve closing the NRP circuit on the channel’s own R_p.

Each global step applies the same nested Strang sequence as solve_split() – half hydro, half conduction, the per-node reaction, then both half-steps again in reverse – with the circuit closed between them.

The field is set by the NRP circuit closed on the channel’s plasma resistance R_p = gap / integral(sigma dA), taken over the solver’s own cells by cross_section_integral(), so as the channel broadens its falling R_p loads the line and pulls the field down (unlike a prescribed E(t)). R_p is frozen over each global step (explicit circuit coupling); the generator waveform V_g(t) still varies continuously within the step via a 2-point interpolation of V_p (which keeps the stateful circuit history monotonic). Records R_p, V_p, total current I = V_p/R_p, and E = V_p/gap.

Parameters:

  • t_out (numpy.ndarray) – Output / split times [s]; the step size is set by consecutive entries.

  • circuit (TransmissionLineResistiveCircuit) – Electric circuit object (generator + cable). Same single-source circuit used by Isomass2TVolumeReactor.

  • gap (float) – Inter-electrode gap [m], used for R_p = gap / integral(sigma dA) and E = V_p / gap.

  • Tg0 (numpy.ndarray, shape (npts,)) – Initial temperatures [K].

  • Te0 (numpy.ndarray, shape (npts,)) – Initial temperatures [K].

  • Y0 (numpy.ndarray, shape (npts, nsp)) – Initial mass fractions.

  • rtol (float, optional) – Per-node integrator tolerances.

  • atol (float, optional) – Per-node integrator tolerances.

Returns:

History with keys t, Tg, Te, Y, ne, plus the circuit diagnostics Rp [Ohm], Vp [V], I [A], E [V/m].

The state arrays are the state at t. The four circuit arrays describe the interval ending at t instead: Rp is the resistance held across it, evaluated from the state at its start, and I is the current that held resistance carried at its end. Plotted against t the resistance therefore leads the current by one interval. Deriving Rp from the end state instead would remove the lag and make I a current the circuit never drove – during breakdown the resistance falls far enough within one step that the peak would read 141 A against the 47 A actually applied – so the lag is kept and stated rather than traded for that.

Rp_state is the other one: the resistance of the frame’s own recorded fields, which is what recomputing gap / cross_section_integral(conductivity_profile(...)) from the history gives. Use Rp with I and Vp; use Rp_state to compare against a model whose resistance is read after its own advance, such as the 0D reactor’s plasma_resistance.

Return type:

dict