Examples

IDL format output loader

Read data

file = "1d_bin.out";
bd = load(file);
bd = load(file, verbose=true);
bd = load(file, npict=1);

3D structured spherical coordinates

file = "3d_structured.out";
bd = load(file, verbose=false);

log file

logfilename = "shocktube.log";
head, data = readlogdata(logfilename)

Data Extraction

Checking variable range

get_var_range(bd, "rho")

Raw variables

Note that the variable names for queries must be in lowercase!

ρ = getvar(bd, "rho")
bd["rho"]

Extracting data at a given location

loc = Float32[0.0, 0.0] # The type determines the output type
d = interp1d(bd, "rho", loc)

Extracting data along a given line

point1 = Float32[-10.0, -1.0]
point2 = Float32[10.0, 1.0]
w = interp1d(bd, "rho", point1, point2)

bd["rho"][X=-10 .. 10, Y=Near(0.0)]

Extracting data using DimensionalData

We can also use Selectors from DimensionalData for extracting data. Note that the Selectors need to be imported from Batsrus.jl; alternatively you can simply using DimensionalData.

bd["rho"][X=At(0.0), Y=At(0.0)]
bd["rho"][X=-10 .. 10, Y=Near(0.5)]
bd["rho"][X=-10 .. 10, Y=-0.5 .. 0.5]

Derived variables

We provide utility methods get_magnitude, get_magnitude2, and fill_vector_from_scalars for vector processing:

Bmag = get_magnitude(bd, :B)
B2 = get_magnitude2(bd, :B)
Bvec = Batsrus.fill_vector_from_scalars(bd, :B)
paniso0 = get_anisotropy(bd, 0)

These are built upon get_vectors. fill_vector_from_scalars is slower than get_vectors since it involves additional array allocations. Here is a full list of predefined derived quantities in get_vectors:

Derived variable name	Meaning	Required variable
:B	Magnetic field vector	Bx, By, Bz
:E	Electric field vector	Ex, Ey, Ez
:U	Velocity vector	Ux, Uy, Uz
:U0	Electron velocity vector	UxS0, UyS0, UzS0
:U1	Proton velocity vector	UxS1, UyS1, UzS1

Output format conversion

We can convert 2D/3D BATSRUS outputs *.dat to VTK formats. It uses the VTK XML format writer writeVTK to generate files for Paraview and Tecplot. The default converted filename is out.vtu.

ASCII Tecplot file (supports both tec and tcp) and binary Tecplot file (set DOSAVETECBINARY=TRUE in BATSRUS PARAM.in):

file = "x=0_mhd_1_n00000050.dat"
convertTECtoVTU(file)

3D structured IDL file (gridType=:vti returns image vti file, gridType=:vtr returns rectilinear vtr file, gridType=:vts returns structured vts file):

file = "3d_structured.out"
convertIDLtoVTK(file, gridType=:vti)

3D unstructured IDL file together with header and tree file:

filetag = "3d_var_1_n00002500"
convertIDLtoVTK(filetag)

Note

The file suffix should not be provided for this to work correctly!

Multiple files:

dir = "./"
filenames = filter(file -> startswith(file, "3d") && endswith(file, ".dat"), readdir(dir))
filenames = dir .* filenames

for filename in filenames
   convertTECtoVTU(filename, filename[1:end-4])
end

Processing multiple files with threads in parallel:

dir = "./"
filenames = filter(file -> startswith(file, "3d") && endswith(file, ".dat"), readdir(dir))
filenames = dir .* filenames

Threads.@threads for filename in filenames
   println("filename=$filename")
   convertTECtoVTU(filename, filename[1:end-4])
end

More examples can be found in examples.

HDF format output loader

filename = "3d__var_1_n00006288.h5"
file = BatsrusHDF5Uniform(filename)

Field extraction

Variable var can be extracted in the whole domain:

var, (xl_new, yl_new, zl_new), (xu_new, yu_new, zu_new) = extract_var(file, "bx")

where (xl_new, yl_new, zl_new) and (xu_new, yu_new, zu_new) return the lower and upper bound, respectively.

Variables within a box region can be extracted as following:

var, (xl_new, yl_new, zl_new), (xu_new, yu_new, zu_new) =
   extract_var(file, "bx"; xmin, xmax, ymin, ymax, zmin, zmax)

Data visualization

We provide plot recipes for Plots.jl, Makie.jl, and wrappers for PyPlot.jl.

The recipes for Plots.jl and Makie.jl will work on all kinds of plots given the correct dimensions, e.g.

using Plots

plot(bd, "p")
contourf(bd, "Mx", xlabel="x")

See the official documentation for Plots.jl for more information.

On the other hand, most common 1D and 2D plotting functions are wrapped over their Matplotlib equivalences through PyPlot.jl. To trigger the wrapper, using PyPlot. Check out the documentation for more details.

Quick exploration of data

Using the same plotting functions as in Matplotlib is allowed, and actually recommended. This takes advantage of multiple dispatch mechanism in Julia. Some plotting functions can be directly called as shown below, which allows for more control from the user. using PyPlot to import the full capability of the package, etc. adding colorbar, changing line colors, setting colorbar range with clim.

For 1D outputs, we can use plot or scatter.

line plot

plot(bd, "p", linewidth=2, color="tab:red", linestyle="--", linewidth=2)

scatter plot

scatter(bd, "p")

For 2D outputs, we can select the following functions:

contour
contourf
imshow
pcolormesh
plot_surface
plot_tricontour
plot_tricontourf
plot_trisurf
tripcolor

with either quiver or streamplot. By default the linear colorscale is applied. If you want to switch to logarithmic, set argument colorscale=:log.

contour

contour(bd, "p")

filled contour

contourf(bd, "p")
contourf(bd, "p"; levels, plotrange=[-10,10,-Inf,Inf], plotinterval=0.1)

surface plot

plot_surface(bd, "p")

triangle surface plot

plot_trisurf(bd, "p")

triangle filled contour plot

tricontourf(bd, "p")

streamline

streamplot(bd, "bx;bz")
streamplot(bd, "bx;bz"; density=2.0, color="k", plotinterval=1.0, plotrange=[-10,10,-Inf,Inf])

quiver (currently only for Cartesian grid)

quiver(bd, "ux;uy"; stride=50)

streamline + contourf

file = "y.out"
bd = load(file)

DN = matplotlib.colors.DivergingNorm
cmap = matplotlib.cm.RdBu_r

contourf(bd, "uxS0", 50; plotrange=[-3,3,-3,3], plotinterval=0.05, norm=DN(0), cmap)
colorbar()
streamplot(bd, "uxS0;uzS0"; density=2.0, color="g", plotrange=[-3,3,-3,3])
xlabel("x"); ylabel("y"); title("Ux [km/s]")

contourf(bd, "uxS0", 50; plotinterval=0.05, norm=DN(0), cmap)
colorbar()
axis("scaled")
xlabel("x"); ylabel("y"); title("uxS0")

For 3D outputs, we may use cutplot for visualizing on a sliced plane, or streamslice to plot streamlines on a given slice.

Finding indexes

To get the index of a certain quantity, e.g. electron number density

ρe_= findfirst(x->x=="rhos0", bd.head.wname)

Get variable range

wmin, wmax = get_var_range(bd, var)

Tracing

The built-in streamplot function in Matplotlib is not satisfactory for accurately tracing. Instead we recommend FieldTracer.jl for tracing fieldlines and streamlines.

An example of tracing in a 2D cut and plot the field lines over contour:

file = "test/y=0_var_1_t00000000_n00000000.out"
bd = load(file)

bx = bd.w[:,:,5]
bz = bd.w[:,:,7]
x  = bd.x[:,1,1]
z  = bd.x[1,:,2]

seeds = select_seeds(x, z; nSeed=100) # randomly select the seeding points

for i in 1:size(seeds)[2]
   xs = seeds[1,i]
   zs = seeds[2,i]
   # Tracing in both direction. Check the document for more options.
   x1, z1 = trace2d_eul(bx, bz, xs, zs, x, z, ds=0.1, maxstep=1000, gridType="ndgrid")
   plot(x1, z1, "--")
end
axis("equal")

Currently the select_seeds function uses pseudo random number generator that produces the same seeds every time.

AMReX Particle Data

We provide support for reading and analyzing AMReX particle data.

Loading Data

To load AMReX particle data:

data = AMReXParticle("path/to/data_directory")

This will parse the header and prepare for lazy loading of particle data.

Phase Space Plotting

We can calculate and plot the phase space density distribution of particles.

First, load the PyPlot extension. The plot_phase function automatically calculates the density and plots it.

using PyPlot

plot_phase(data, "x", "vx"; 
   bins=100, 
   x_range=(-10, 10), 
   y_range=(-5, 5),
   log_scale=true,
   plot_zero_lines=true,
   normalize=true
)

We can also calculate the phase space density histogram directly without plotting:

# 1D density
hist1d = get_phase_space_density(data, "vx")

# 2D density
hist2d = get_phase_space_density(data, "x", "vx"; bins=(100, 50))

# 3D density
hist3d = get_phase_space_density(data, "vx", "vy", "vz"; bins=50)

# Weighted 2D density (automatically detects "weight" component if present)
# or passes weights explicitly if not in data
hist2d_w = get_phase_space_density(data, "v_parallel", "v_perp")

We can also apply coordinate transformations to the particle data.

Transformation with only B field. This decomposes velocity into parallel and perpendicular components relative to B.

transform_b = get_particle_field_aligned_transform([1.0, 0.0, 0.0])

plot_phase(data, "v_parallel", "v_perp"; 
   transform=transform_b,
   bins=50,
   log_scale=true
)

Transformation with both B and E fields. This creates an orthonormal basis ($v_B$, $v_E$, $v_{B \times E}$), where $v_B$ is along B, $v_E$ is along the perpendicular component of E, and $v_{B \times E}$ is along the ExB drift direction.

transform_eb = get_particle_field_aligned_transform([1.0, 0.0, 0.0], [0.0, 1.0, 0.0])

plot_phase(data, "v_B", "v_BxE"; 
   transform=transform_eb,
   bins=50,
   log_scale=true
)

Particle Classification

You can classify particles into Core Maxwellian and Suprathermal populations using classify_particles. This function allows for specifying a spatial region and handling velocity distributions in 1D, 2D, or 3D.

# range can be specified by keywords x_range, y_range, z_range locally
core, halo = classify_particles(data; 
    x_range=(-1.0, 1.0), 
    y_range=(-1.0, 1.0), 
    z_range=(-1.0, 1.0),
    vdim=3,          # Velocity dimension (1, 2, or 3)
    vth=1.0,         # Core thermal velocity (required)
    nsigma=3.0,      # Separation threshold
    bulk_vel=nothing # Auto-detect if nothing
)

The function returns two matrices containing the classified particles. If bulk_vel is not provided, it is automatically estimated from the peak of the velocity distribution.