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

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

• Derived variables

v = getvars(bd, ["Bx", "By", "Bz"])
Bmag = bd["Bmag"]

• 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)

Here is a full list of predefined derived quantities:

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=1 returns rectilinear vtr file, gridType=2 returns structured vts file):

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

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

filetag = "3d_var_1_n00002500"
convertIDLtoVTK(filetag)

• 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

Variables within a box region can be extracted as following:

out, (xl_new, yl_new, zl_new), (xu_new, yu_new, zu_new) =
   extract_field(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

A general plotdata function is provided for quick visualizations using Matplotlib.

• 1D binary

plotdata(bd, "p", plotmode="line")
plotdata(bd, "p", plotmode="linegrid")

• 2D Cartesian (structured)

plotdata(bd, "p bx;by", plotmode="contbar streamover")
plotdata(bd, "p bx;by", plotmode="contbar quiverover")
plotdata(bd, "p bx;by", plotmode="contbar streamover", density=2.0)
plotdata(bd, "p", plotmode="grid")
plotdata(bd, "p", plotmode="contbar", plotrange=[-50., 50., -1., 1.])
plotdata(bd, "p", plotmode="contbar")
plotdata(bd, "p", plotmode="contbarlog")
plotdata(bd, "p", plotmode="surfbar")

• 2D unstructured

plotdata(bd, "rho", plotmode="contbar")
plotdata(bd, "rho", plotmode="trimesh")
plotdata(bd, "rho", plotmode="tricont")

• 2D structured spherical coordinates

plotdata(bd, "rho", plotmode="contbar")

• 3D box

plotdata(bd, "bx", plotmode="contbar", dir="y", sequence=1, level=20)
plotdata(bd, "bx", plotmode="contbar", dir="y", plotrange=[-1.4,-1.1,0.70,0.78])
plt.axis("scaled")

subplot(2,2,(1,3))
cutplot(bd, "Ex"; dir="y", sequence=128, plotrange)

Finding indexes

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

ρe_= findfirst(x->x=="rhoS0", bd.head.wnames)

Multiple dispatch for Matplotlib functions

Using the same plotting functions as in Matplotlib is allowed, and actually recommended. 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.

• line plot

plot(bd, "p", linewidth=2, color="green")
c = plot(bd, "p")
plt.setp(c, linestyle="--", linewidth=2);

• scatter plot

scatter(bd, "p")

• 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")

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.