Plotting Examples

There are a variety of methods that can be used to visualize a network.

import pynucastro as pyna

Standard view

Let’s build a quick network that represents He, C, and O burning

library = pyna.ReacLibLibrary()

sub = library.linking_nuclei(["p", "n", "he4", "c12", "n13", "n14", "o16",
                              "ne20", "ne21", "na23", "mg23", "mg24", "al27",
                              "si27", "si28"])
rc = pyna.RateCollection(libraries=sub)

/opt/hostedtoolcache/Python/3.11.10/x64/lib/python3.11/site-packages/pynucastro/networks/rate_collection.py:576: UserWarning: ReacLib neutron decay rate (<n_to_p_weak_wc12>) does not account for degeneracy at high densities. Consider using tabular rate from Langanke.
  warnings.warn(msg)

The standard network plot of this is:

fig = rc.plot()

Rotated view

Since everything lies close to the diagonal, it is hard to see the structure. We can instead look at a rotated view of the network.

While we’re at it, we’ll also color the rates by reaction rate strength

rho = 1.e4
T = 1.e8
comp = pyna.Composition(rc.get_nuclei())
comp.set_solar_like()

fig = rc.plot(rho, T, comp, rotated=True)

This looks a bit crowded, since it is showing all of the links to \({}^4\mathrm{He}\) from reactions of the form \(A(x, \alpha)B\). We can hide those links with hide_xalpha=True

fig = rc.plot(rho, T, comp, rotated=True, hide_xalpha=True)

We can also plot the edges curved, to help distinguish the forward and reverse rates

fig = rc.plot(rho, T, comp, rotated=True, hide_xalpha=True,
              curved_edges=True)

We can also label edges if desired

edge_labels = {(pyna.Nucleus("he4"), pyna.Nucleus("c12")):
               r"$\alpha(\alpha\alpha,\gamma){}^{12}\mathrm{C}$"}

fig = rc.plot(rho, T, comp, rotated=True, hide_xalpha=True,
              edge_labels=edge_labels, size=(1200, 700),
              curved_edges=True)

Filters

We can use a rate filter to show only certain links. For example, to show only links involving proton capture, we can do:

fig = rc.plot(rho, T, comp, rotated=True, hide_xalpha=True,
              rate_filter_function=lambda r: pyna.Nucleus("p") in r.reactants)

We can also highlight rates using the highlight_filter_function mechanism. This takes a function that operates on a rate and returns True if we wish to highlight the rate (i.e., as if we used a highlighter marker over the line).

Here’s an example where we highlight the rates that have \(Q > 0\)

fig = rc.plot(rotated=True, hide_xalpha=True, size=(1500, 600),
              highlight_filter_function=lambda r: r.Q > 0,
              curved_edges=True)

Plotting the Jacobian

For implicit integration, we solve:

\[\dot{\bf Y} = {\bf f}({\bf Y})\]

and the Jacobian, \({\bf J} = \partial {\bf f} / \partial {\bf Y}\) plays an important role in solving for the update.

We can visualize the Jacobian of our network easily. By default, we will show all elements where

\[\frac{\mathrm{max}_{i,j} (|J_{ij}|)}{|J_{ij}|} > 10^{10}\]

but this ratio can be set via the rate_scaling keyword argument.

fig = rc.plot_jacobian(rho, T, comp, rate_scaling=1.e50)

Plot nuclides on a grid

Nuclides in a network may also be visualized as cells on a grid of Z vs. N, colored by some quantity. This can be more interpretable for large networks. Calling gridplot without any arguments will just plot the grid - to see anything interesting we need to supply some conditions. Here is a plot of nuclide mass fraction on a log scale, with a 36 square inch figure:

rate_names = ["c12(p,g)n13",
              "c13(p,g)n14",
              "n13(,)c13",
              "n13(p,g)o14",
              "n14(p,g)o15",
              "n15(p,a)c12",
              "o14(,)n14",
              "o15(,)n15"]
rates = library.get_rate_by_name(rate_names)
rc = pyna.RateCollection(rates=rates)

rho = 1.e4
T = 1.e8
comp = pyna.Composition(rc.get_nuclei())
comp.set_solar_like()

fig = rc.gridplot(comp=comp, color_field="X", scale="log", area=36)

The plot is configurable through a large number of keyword arguments. Here we want to look at the rates at which nuclides are being created or destroyed, so we color by \(\dot{Y}\), the rate of change of molar abundance. Density and temperature need to be supplied to evaluate the rates. A full list of valid keyword arguments can be found in the API documentation.

fig = rc.gridplot(comp=comp, rho=rho, T=T, color_field="ydot", area=36,
                  cmap="RdGy", cbar_bounds=(-0.1, 0.1))

Unlike the network plot, this won’t omit hydrogen and helium by default. To just look at the heavier nuclides, we can define a function to filter by proton number. Here we also plot activity, which is the sum of all rates contributions to Ydot ignoring the sign.

ff = lambda nuc: nuc.Z > 2
fig = rc.gridplot(comp=comp, rho=rho, T=T, color_field="activity", scale="log",
                  filter_function=ff, area=20, cmap="Blues")

Network Chart

A network chart shows how each rate affects each nucleus

fig = rc.plot_network_chart(rho, T, comp)