Network

Network is the entry point for every JAFF workflow and the hub the other objects hang off. One constructor call reads a network file, parses it into typed Species, Reactions, and Elements catalogues, rewrites every rate into a common symbolic form, validates the chemistry, and exposes the symbolic ODE/flux/radiation expressions that downstream code generation consumes.

from jaff import Network

net = Network("networks/h_photoionization/h_photo.jet")

What loading does

A single Network(...) call runs a fixed pipeline. Knowing the phases explains where every attribute comes from:

1. Read — the file is resolved and, if it is a binary .jaff, deserialized directly; otherwise the text parser auto-detects the format (KROME, PRIZMO, UDFA, KIDA, UCLChem, …).
2. Parse & assemble — reactions are read one by one. Each new species name becomes a Specie appended to net.species; each line becomes a Reaction appended to net.reactions. Any accompanying .jfunc auxiliary file is loaded here and its custom rates, heating/cooling, and radiation terms are wired in.
3. Standardize symbols — every rate and energy expression is rewritten so that number densities use the common nden[i] representation (shorthands like nh, ntot, n_X are expanded).
4. Validate — mass/charge conservation, sinks/sources, recombinations, isomers, and duplicate reactions are checked.
5. Finalize — integer stoichiometry matrices are built and the unique Elements set is derived from the species.

After this, the network is a fully assembled, queryable model.

Constructor

Network(
    fname,
    errors=False,
    label=None,
    funcfile=None,
    replace_nH=True,
    rad_bands=[],
    rad_powerlaw_index=0,
    rad_energy_density=False,
    c=29979245800.0,
)

Examples

# Minimal load
net = Network("networks/GOW/GOW.jet")

# Strict validation — abort on the first conservation failure
net = Network("networks/GOW/GOW.jet", errors=True)

# Custom label and explicit auxiliary-function file
net = Network(
    "networks/GOW/GOW.jet",
    label="GOW-2017",
    funcfile="networks/GOW/GOW.jfunc",
)

# With radiation bands (needed for photochemistry during code generation)
net = Network(
    "networks/h_photoionization/h_photo.jet",
    rad_bands=[13.6, "inf"],
    rad_energy_density=False,
)

Key attributes

net = Network("networks/h_photoionization/h_photo.jet")

net.label                   # 'h_photo'
net.species.count           # 3
net.reactions.count         # 2
net.reactant_matrix.shape   # (2, 3)   — 2 reactions × 3 species

The three catalogues are the same typed collections documented on Species, Reactions and Elements — anything there works straight off the network:

net.species.names()                 # ['H', 'H+', 'e-']
net.species.charged("name")         # ['H+', 'e-']
net.reactions.verbatim()            # ['H -> H+ + e-', 'H+ + e- -> H']
net.reactions.photo_reactions()     # [ReactionObject(H -> H+ + e-)]
net.elements.symbols()              # ['H']

Stoichiometry matrices

reactant_matrix[i, j] counts how many times species j appears as a reactant in reaction i (and product_matrix likewise for products). Row order is reaction index; column order is species index.

net.reactant_matrix
# [[1 0 0]      reaction 0: H -> H+ + e-   consumes 1× H (species 0)
#  [0 1 1]]     reaction 1: H+ + e- -> H   consumes 1× H+ and 1× e-

net.product_matrix
# [[0 1 1]
#  [1 0 0]]

The net stoichiometry product_matrix - reactant_matrix is the change in each species per reaction — the backbone of the ODE system below.

Number densities: the nden representation

Inside JAFF every species number density is a reference into one SymPy MatrixSymbol("nden", n_species, 1). Species X is nden[idx_X, 0], where the index is the species' position in net.species. This is why rate and ODE expressions print terms like nden[1, 0] rather than a named density — it keeps the whole network on one indexable vector the code generator can emit as an array.

During loading (phase 3), convenience shorthands in rate expressions are expanded into this form:

For example, a cosmic-ray ionization rate written with nh shorthand comes out as a weighted nden sum (note 2*nden[3, 0] — H₂ contributes two H atoms):

g = Network("networks/GOW/GOW.jet")
g.reactions.with_rtype("cosmic_ray")[0].rate
# crate*(1.5*nden[0, 0]/(nden[0, 0] + nden[1, 0] + 2*nden[3, 0] + ...) + ...)

Pass replace_nH=False to keep nh / nhe as opaque free symbols instead of expanding them — useful when an external driver supplies those densities directly.

Symbolic expressions

Once loaded, the network exposes its dynamics as SymPy expressions.

Species ODEs

Network.sodes() returns one expression per species — the net rate of change of its number density (cm⁻³ s⁻¹), summed over every reaction it takes part in:

where \(f_i\) is the flux of the \(i^{\text{th}}\) reaction.

for specie, expr in zip(net.species, net.sodes()):
    print(f"d[{specie}]/dt = {expr}")

# d[H]/dt  = 1.659e-10*nden[1, 0]*nden[2, 0]/tgas**0.7 - photorates(1, 13.6, 1.0e+99)*nden[0, 0]
# d[H+]/dt = -1.659e-10*nden[1, 0]*nden[2, 0]/tgas**0.7 + photorates(1, 13.6, 1.0e+99)*nden[0, 0]
# d[e-]/dt = -1.659e-10*nden[1, 0]*nden[2, 0]/tgas**0.7 + photorates(1, 13.6, 1.0e+99)*nden[0, 0]

Fluxes and radiation moments

net.sfluxes()       # list[Expr] — per-reaction flux  k_i * nden[r1] * nden[r2] ...
net.sradodes(0)     # list[Expr] — radiation moment ODEs (order 0), if rad_bands set

A flux is just the reaction rate times its reactant densities; the species ODEs are signed sums of these. See sfluxes and sradodes for details.

Validation

JAFF validates the network during loading. With errors=False (default) each problem is logged as a warning; with errors=True the process exits on the first one. The individual checks can also be re-run by hand — each takes an errors flag and logs its findings (it returns None, it does not return a list):

net.check_sink_sources(errors=False)     # species never produced, or never consumed
net.check_recombinations(errors=False)   # cations lacking an electron recombination
net.check_isomers(errors=False)          # species sharing one atomic composition
net.check_unique_reactions(errors=False) # duplicate reactions (same species, type, T-range)

Mass and charge conservation are validated per reaction at construction time — see Reactions.

Comparing networks

compare_reactions and compare_species diff two loaded networks, logging what is shared and what is unique to each. verbosity=1 (default) also prints the full sets:

net2 = Network("networks/GOW/GOW.jet")

net.compare_reactions(net2)             # reactions in one but not the other
net.compare_species(net2, verbosity=0)  # counts only, no per-item listing

Comparison is by serialized form, so isomer/duplicate handling is consistent with the rest of JAFF.

Export and caching

Rate tables

Reaction rates can be tabulated against temperature and written to HDF5 or text for use by external solvers:

net.to_hdf5("rates.hdf5", T_min=10, T_max=1e4, nT=200, err_tol=1e-3)
net.to_txt("rates.txt", T_min=10, T_max=1e4, nT=200)
net.write_table("rates.hdf5", format="auto")  # general form; format from extension (.hdf5/.hdf/.txt)

err_tol drives adaptive temperature sampling (set it to None for a fixed grid). See to_hdf5 and to_txt.

JAFF binary format

to_jaff serializes an already-loaded network to a gzip-compressed JSON file (conventional extension .jaff). Re-loading skips parsing and the expensive SymPy assembly, so it is the recommended cache for large networks:

net = Network("networks/GOW/GOW.jet")
net.to_jaff("gow.jaff")        # save

net2 = Network("gow.jaff")     # reload — same species, reactions, and ODEs