Output & visualization¶

To enable the runtime output of the simulation data, configure the code with the -D output=ON flag. As a backend Entity uses the open-source ADIOS2 library compiled in-place. The output is written in the ADIOS2 format called BP5, but HDF5 is also available (but not recommended).

The output is configured using the following configurations in the input file:

[simulation]
  name   = "MySimulation" # (5)!

  # ...
[[particles.species]]
  tracking = true # (15)!

  # ...
[output]
  format = "BPFile" # (2)!
  interval = 100 # (3)!
  interval_time = 0.1 # (8)!

  [output.fields]
    quantities = ["B", "E", "Rho_1_2", "..."] # (1)!
    stride = 2 # (9)!
    mom_smooth = 2 # (4)!

  [output.particles]
    species = [1, 2, 4] # (7)!
    stride = 10 # (6)!

  [output.spectra]
    e_min = 1e-2 # (12)!
    e_max = 1e3
    log_bins = true # (13)!

  [output.stats]
    quantities = ["N", "Npart", "ExB", "J.E"] # (14)!

  [output.debug]
    as_is = false # (10)!
    ghosts = false # (11)!

fields to write
output format (current supported: "BPFile"/"HDF5", or "disabled" for no output)
output interval (in the number of time steps)
smoothing stencil size for moments (in the number of cells) [defaults to 1]
title is used for the output filename
stride used for particle output (write every prtl_stride-th particle) [defaults to 100]
particle species to output
output interval in time units (overrides interval if specified)
stride used for field output (write every fields_stride-th cell) [defaults to 1]
write the field quantities as-is (without conversion/interpolation) [defaults to false]
write the ghost cells [defaults to false]
Min/max energies for binning the energy distribution [default to 1e-3 -> 1e3]
whether to use logarithmic energy bins or linear
box reduced quantities to output as stats
enable tracking for a given particle species

For the full list, please look at the input.example.toml file or refer to the following section.

Following is the list of all supported fields

Field name	Description	Normalization
`E`	Electric field (all components)	$B_0$
`B`	Magnetic field (all components)	$B_0$
`D`	GR: electric field (all components)	$B_0$
`H`	GR: aux. magnetic field (all components)	$B_0$
`J`	Current density (all components)	$4\pi q_0 n_0$
`Rho`	Mass density	$m_0 n_0$
`Charge`	Charge density	$q_0 n_0$
`N`	Number density	$n_0$
`V` 1.2.0	Mean 3-velocity	dimensionless
`Nppc`	Raw number of particles per cell	dimensionless
`Tij`	Energy-momentum tensor (all components)	$m_0 n_0$
`divE` 1.2.0	Divergence of $E$	arb. units
`divD` 1.2.0	GR: divergence of $D$	arb. units
`A`	GR: 2D vector potential $A_\varphi$	arb. units

and particle quantities

Particle quantity	Description	Units
`X`	Coordinates (all components)	physical
`U`	Four-velocities (all components)	dimensionless
`W`	Weights	dimensionless
`PLDR` 1.3.0	Real-valued payloads	arbitrary
`PLDI` 1.3.0	Integer-valued payloads	arbitrary
`RNK` 1.3.0	Meshblock rank the particle was created (if MPI is ON)	--
`IDX` 1.3.0	Index of the particle on the given rank	--

1.2.0 The code also has an output of box-averaged stats into a .csv file, which are simply scalars per each output timestep. The following quantities can be computed

Box-reduced quantity	Description	Units
`E^2`	Mean $E^2$	$B_0^2$
`B^2`	Mean $B^2$	$B_0^2$
`ExB`	Mean $\bm{E}\times \bm{B}$	$B_0^2$
`J.E`	Mean $\bm{J}\cdot \bm{E}$	$4\pi q_0 n_0 B_0$
`N`	Mean $n$	$n_0$
`Npart`	Total # of particles	dimensionless
`Rho`	Mean mass density	$m_0 n_0$
`Charge`	Mean charge density	$q_0 n_0$
`Tij`	Mean energy-momentum tensor (all components)	$m_0 n_0$

"Mean" in this context refers to volume-averaging: i.e., $\langle E_x^2 \rangle = V^{-1}\int \sqrt{h} d^3 \bm{x}~ E_x^2 $, or $\langle T^{ij}\rangle \equiv V^{-1} \int d^3\bm{u} \sqrt{h} d^3 \bm{x} ~(u^i u^j / u^0) f(\bm{u}) $, where $V\equiv \int \sqrt{h} d^3\bm{x}$. As such, these values (except for Npart) are insensitive to the resolution of the grid or the number of particles per cell.

Refining moments for the output

One can specify particular components to output for the Tij fields/stats: T0i will output the T00, T01, and T02 components, while Tii will output only the diagonal components: T11, T22, and T33, and Tij will output all the 6 components. For quantities computed from particles (moments of the distribution), one can also specify the particle species which will be used to compute the moments: Rho_1 (density of species 1), N_2_3 (number density of species 2 and 3), Tij_1_3 (energy-momentum tensor for species 1 and 3), etc.

All of the vector fields are interpolated to cell centers before the output, and converted to orthonormal basis. The particle-based moments are smoothed with a stencil (specified in the input file; mom_smooth) for each particle.

In addition, one can write custom user-defined field quantities to the output with the fields or stats. Refer to the following section for more details.

Can one track particles at different times?

1.3.0 Yes! Simply enable particle tracking for a particular species. Then each particle is uniquely identified by a combination of IDX and RNK (if no MPI is used, then only IDX is sufficient). nt2py already automatically combines the variables producing a unique id for each particle (for the species where tracking is enabled). However, keep in mind, that the simulations are not reproducible and will unfortunately never be due to limitations imposed by the nature of GPU computations.

`nt2py`¶

We provide the nt2py python package to help easily access and manipulate the simulation data. nt2py package uses the dask and xarray libraries together with adios2 and/or h5py to lazily load the output data and provide a convenient interface for the data analysis and quick visualization.

To start using nt2py, it is recommended to create a python virtual environment and install the required packages:

python3 -m venv .venv
source .venv/bin/activate # (1)!
pip install nt2py # (2)!

Now all the packages will be installed in the .venv directory which you can remove at any time without affecting the system.
If you plan to use jupyter you might also need to run the following pip install jupyterlab ipykernel.

Usage¶

Simply pass the location to the data when initializing the main Data object:

import nt2

data = nt2.Data("path/to/data")

The data is stored in specialized containers which can be accessed via corresponding attributes:

data.fields      # < xr.Dataset
data.particles   # < special object which returns a pd.DataFrame when .load() is called
data.spectra     # < xr.Dataset
data.diagnostics # < pd.DataFrame

If using Jupyter notebook, you can quickly preview the loaded metadata by simply running a cell with just data in it (or in regular python, by doing print(data)).

Note, that by default, the hdf5 support is disabled in nt2py (i.e., only ADIOS2 format is supported). To enable it, install the package as pip install "nt2py[hdf5]" instead of simply pip install nt2py.

Accessing the data¶

Fields and spectra are stored as lazily loaded xarray datasets (a collection of equal-sized arrays with shared axis coordinates). You may access the coordinates in each dimension using .coords:

data.fields.coords
data.spectra.coords

Individual arrays can be requested by simply using, e.g., data.fields.Ex etc. One can also use slicing/selecting via the coordinates, i.e.,

data.fields.sel(t=5, method="nearest")

accesses all the fields at time t=5 (using method="nearest" means it will take the closest time to value 5). You may also access by index in each coordinate:

data.fields.isel(x=-1)

accesses all the fields in the last position along the x coordinate.

Note that all these operations do not load the actual data into memory; instead, the data is only loaded when explicitly requested (i.e., when plotting or explicitly calling .values or .load().

Particles are stored in a special lazy container which acts very similar to xarray; you can still make selections using specific queries. For instance,

data.particles.sel(sp=[1, 2, 4]).isel(t=-1)

selects all the particles of species 1, 2, and 4 on the last timestep. The loading of the data itself is done by calling: .load() method, which returns a simple pandas dataframe.

Plotting¶

Plot a field (in Cartesian coordinates) at a specific time (or output step):

data.fields.Ex.sel(t=10.0, method="nearest").plot() # time ~ 10
data.fields.Ex.isel(t=5).plot()                     # output step = 5

Plot a slice or time-averaged field quantities:

data.fields.Bz.mean("t").plot()
data.fields.Bz.sel(t=10.0, x=0.5, method="nearest").plot()

Plot in spherical coordinates (+ combine several fields):

e_dot_b = (data.fields.Er * data.fields.Br +\
           data.fields.Eth * data.fields.Bth +\
           data.fields.Eph * data.fields.Bph)
bsqr = data.fields.Br**2 + data.fields.Bth**2 + data.fields.Bph**2
# only plot radial extent of up to 10
(e_dot_b / bsqr).sel(t=50.0, method="nearest").sel(r=slice(None, 10)).polar.pcolor()

You can also quickly plot the fields at a specific time using the handy .inspect accessor:

data.fields\
    .sel(t=3.0, method="nearest")\
    .sel(x=slice(-0.2, 0.2))\
    .inspect.plot(only_fields=["E", "B"])
# Hint: use `<...>.plot?` to see all options

Or if no time is specified, it will create a quick movie (need to also provide a name in that case):

data.fields\
    .sel(x=slice(-0.2, 0.2))\
    .inspect.plot(name="inspect", only_fields=["E", "B", "N"])

You can also create a movie of a single field quantity (can be custom):

(data.fields.Ex * data.fields.Bx).sel(x=slice(None, 0.2)).movie.plot(name="ExBx")

For particles, one can also make 2D phase-space plots:

data.particles.sel(sp=1).sel(t=1.0, method="nearest").phase_plot(
    x_quantity=lambda f: f.x,
    y_quantity=lambda f: f.ux,
    xy_bins=(np.linspace(0, 60, 100), np.linspace(-2, 2, 100)),
)

or a spectrum plot:

data.particles.sel(sp=[1, 2]).sel(t=1.0, method="nearest").spectrum_plot()

You may also combine different quantities and plots (e.g., fields & particles) to produce a more customized movie:

def plot(t, data):
    fig, ax = plt.subplots()
    data.fields.Ex.sel(t=t, method="nearest").sel(x=slice(None, 0.2)).plot(
        ax=ax, vmin=-0.001, vmax=0.001, cmap="BrBG"
    )
    prtls = data.particles.sel(t=t, method="nearest").load()
    ax.scatter(prtls.x, prtls.y, c="r" if prtls.sp == 1 else "b")
    ax.set_aspect(1)
data.makeMovie(plot)

You may also access the movie-making functionality directly in case you want to use it for other things:

import nt2.plotters.export as nt2e

def plot(t):
  ...

#             this will be the array of `t`-s passed to `plot`
#                           |
#                           V
nt2e.makeFrames(plot, np.arange(100), "myAnim")
nt2e.makeMovie(
    input="myAnim/", output="myAnim.mp4", number=5, overwrite=True
)

# or combined together
nt2e.makeFramesAndMovie(
    name="myAnim", plot=plot, times=np.arange(100)
)

Raw readers¶

In case you want to access the raw data without using nt2py's xarray/dask lazy-loading, you may do so by using the readers. For example, for ADIOS2 output data format:

import nt2.readers.adios2 as nt2a

# define a reader
reader = nt2a.Reader()

# get all the valid steps for particles
valid_steps = reader.GetValidSteps("path/to/sim", "particles")

# get all variable names which have prefix "p" at the first valid step
variable_names = reader.ReadCategoryNamesAtTimestep(
    "path/to/sim", "particles", "p", valid_steps[0]
)

# convert the variable set into a list and take the first element
variable = list(variable_names)[0]

# read the actual array from the file
reader.ReadArrayAtTimestep(
    "path/to/sim", "particles", variable, valid_steps[0]
)

There are many more functions available within the reader. For hdf5, you can simply change the import to nt2.readers.hdf5, and the rest should remain the same.

CLI¶

Since version 1.0.0, nt2py also offers a command-line interface, accessed via nt2 command. To view all the options, simply run:

nt2 --help

The plotting routine is pretty customizable. For instance, if the data is located in myrun/mysimulation, you can inspect the content of the data structure using:

nt2 show myrun/mysimulation

Or if you want to make a quick plot (a-la inspect discussed above) of the specific quantities, you may simply run:

nt2 plot myrun/mysimulation --fields "E.*;B.*" --isel "t=5" --sel "x=slice(-5, None); z=0.5"

This plots the 6-th snapshot (t=5) of all the E and B field components, sliced for x > -5, and at z = 0.5 (notice, that you can use both --isel and --sel). If instead, you prefer to make a movie, simply do not specify the time:

nt2 plot myrun/mysimulation --fields "E.*;B.*" --sel "x=slice(-5, None); z=0.5"

If you want to only install the CLI, without the library itself, you may do that via pipx: pipx install nt2py.

Field name	Description	Normalization
`E`	Electric field (all components)	\(B_0\)
`B`	Magnetic field (all components)	\(B_0\)
`D`	GR: electric field (all components)	\(B_0\)
`H`	GR: aux. magnetic field (all components)	\(B_0\)
`J`	Current density (all components)	\(4\pi q_0 n_0\)
`Rho`	Mass density	\(m_0 n_0\)
`Charge`	Charge density	\(q_0 n_0\)
`N`	Number density	\(n_0\)
`V` 1.2.0	Mean 3-velocity	dimensionless
`Nppc`	Raw number of particles per cell	dimensionless
`Tij`	Energy-momentum tensor (all components)	\(m_0 n_0\)
`divE` 1.2.0	Divergence of \(E\)	arb. units
`divD` 1.2.0	GR: divergence of \(D\)	arb. units
`A`	GR: 2D vector potential \(A_\varphi\)	arb. units

Box-reduced quantity	Description	Units
`E^2`	Mean \(E^2\)	\(B_0^2\)
`B^2`	Mean \(B^2\)	\(B_0^2\)
`ExB`	Mean \(\bm{E}\times \bm{B}\)	\(B_0^2\)
`J.E`	Mean \(\bm{J}\cdot \bm{E}\)	\(4\pi q_0 n_0 B_0\)
`N`	Mean \(n\)	\(n_0\)
`Npart`	Total # of particles	dimensionless
`Rho`	Mean mass density	\(m_0 n_0\)
`Charge`	Mean charge density	\(q_0 n_0\)
`Tij`	Mean energy-momentum tensor (all components)	\(m_0 n_0\)