Manual

Here we demonstrate some basic usages of processing Vlasiator output. A quick self-contained demo is shown with Pluto. For complete description of the arguments, please refer to the API documents or type ?function_name to display help message in the REPL.

Common physical constants

A bunch of physical constants are predefined in Vlasiator.jl. To use them, you need to import explicitly, e.g. using Vlasiator: RE or prepend the module name like Vlasiator.RE.

Loading VLSV data

A quick way to obtain small test data can be found in F&Q. Larger open-access data can be found from the references in Vlasiator publications, e.g. Takahashi+ 2020. In this tutorial we are using the 2D test file bulk.2d.vlsv from vlsv_data if not specified.

• Read meta data

file = "bulk.2d.vlsv"
meta = load(file)

The VLSV meta data contains information of file name, variable names, ordinary cell ID list, mesh sizes, species, and velocity cell structures. It is often the first argument for methods defined in Vlasiator.jl.

• Read parameter

For convenience we support the do-block syntax that automatically closes the file stream.

t = load(file) do meta
   readparameter(meta, "time")
end

• Read variable meta data

readvariablemeta(meta, "proton/vg_rho")

A list of utility functions has been implemented for checking variable status. See here for the full list.

• Read variable

data = meta["CellID"]
# Or equivalently
data = readvariable(meta, "CellID")

The variable reading is designed for cells, which takes cell ID(s) as the 3rd argument if specified. The same interface works for both DCCRG grid (for storing cell centered quantities like plasma moments) and FS grid (for storing field solver related quantities on a uniform high resolution mesh) variables. By default the returned DCCRG variable array is sorted by cell IDs. If in any case you want the original unsorted version as being stored in the file, use readvariable(meta, var, false).

• Get variable at a given location

loc = [2.0, 0.0, 0.0]
id = getcell(meta, loc)
readvariable(meta, "CellID", id)

• Get variable along a line between two points

using Vlasiator: RE # Earth radii
point1 = [12RE, 0, 0]
point2 = [15RE, 0, 0]
cellids, distances, coords = getcellinline(meta, point1, point2)
var_extract = readvariable(meta, "VA", cellids)

• Extract variable at a static cell ID from a sequence of files under the same grid

files = ["bulk.2d.vlsv"]
var = "CellID"
id = 1
extractsat(files, var, id)

• Downsample field solver variable to DCCRG grid

Vlasiator.downsample_fg(meta, "fg_e")

• Upsample DCCRG variable to field solver grid

data = Vlasiator.read_variable_as_fg(meta, "proton/vg_rho")
# Equivalent to the above, but faster
data = let
   tmp = Vlasiator.fillmesh(meta, ["proton/vg_rho"]; maxamronly=true)[1][1][1]
   reshape(tmp, size(tmp)[2:end])
end

This is useful when corresponding DCCRG variables are not saved, or a uniform mesh is required for further analysis.

• Compare VLSV files

One may want to check if two vlsv files are identical. This is tricky because

1. the structure of VLSV format does not guarantee the writing order in parallel processing;
2. numerical error accumulates with floating point representation, especially with fastmath option.

The issame method does not check quantities that are related to the MPI writing sequence: for some reasons, even file sizes may vary depending on the number of MPI processes!

issame(file1, file2)

There is an optional third argument to issame for setting the relative difference tolerance, with default being 1e-4. In practice relative difference works better for "large" numbers, and absolute difference works better for "small" numbers.

Computing derived quantities

Vlasiator.jl is capable of computing plasma moments and some predefined derived quantities and saving them directly into VLSV files. To avoid confusion about variable names, the conventions are

• raw quantities are all lowercases;
• all predefined derived variable names start with a capital letter;
• exceptions are for aliases (e.g. "n").

To obtain a derived quantity, use either a key of string or symbol,

beta = meta["Beta"]
VA = meta[:VA]

Here is a full list of available quantities^[1]:

which can also be found as keys of dictionary in vlsvvariables.jl.

Velocity space moments

We can also calculate plasma moments from the saved VLSV velocity space distributions.

meta = load("bulk.1d.vlsv")
cellid = 5
# VDF cell indexes and values, with sparsity
vcellids, vcellf = readvcells(meta, cellid; species="proton")

getdensity(meta, vcellf)

getvelocity(meta, vcellids, vcellf)
# pressure tensor components Pxx, Pyy, Pzz, Pyz, Pzx, Pxy
getpressure(meta, vcellids, vcellf)
# heat flux components qⱼⱼᵢ
getheatfluxvector(meta, vcellids, vcellf)

To obtain the original ordering of velocity cells,

vcellids_original = Vlasiator.reorder(meta.meshes["proton"], vcellids)

Non-Maxwellianity represents the deviation from a Maxwellian distribution. Currently we have implemented a monitor quantity named "Maxwellianity", which is defined as $-ln \big[ 1/(2n) \int |f(v) - g(v)| dv \big]$, where n is the density, f(vᵢ) is the actual VDF value at velocity cell i, and g(vᵢ) is the analytical Maxwellian (or strictly speaking, normal) distribution with the same density, bulk velocity and scalar pressure as f.

getmaxwellianity(meta, vcellids, vcellf)

The value ranges from [0, +∞], with 0 meaning not Maxwellian-distributed at all, and +∞ a perfect Maxwellian distribution. An alternative measure is the KL-divergence borrowed from statistics, where we select the reference distribution g to be again a Maxwellian derived from the actual distribution f:

getKLdivergence(meta, vcellids, vcellf)

In this case, 0 indicates perfect Maxwellian while +∞ indicates largest deviation. There is no unique definition of non-Maxwellianity, and we are still trying to see which one works better for describing plasma behaviors.

Sometimes it may be useful to recover the full 3D array of VDFs:

f = Vlasiator.reconstruct(meta.meshes["proton"], vcellids, vcellf)

However, usually in practice there would be only about 1% nonzero values. The moments and maxwellianity calculations above all have an alternative form of using reconstructed VDFs as inputs.

Plotting

Vlasiator.jl does not include any plotting library as explicit dependency, but it offers plotting recipes/wrappers once the target plotting package is used.

Currently PyPlot provides the most complete and fine-tuned plotting capabilities. Plotting with PyPlot by accepting MetaVLSV as the first argument is supported via VlasiatorPyPlot. Plots is a collection of plotting libraries with a uniform frontend, but it lacks fine-tuned supports and consistent APIs between different backends. Makie, a native Julia plotting library, is also supported via package extension after Julia 1.9 and a package VlasiatorMakie.jl before Julia 1.9. Without generating an system image from PackageCompiler, it would take ~20s for the first plot on Julia 1.9. However, Makie has made nice progress in layouts, widgets, docs, and all the tiny things, which makes it a strong candidate for the de facto plotting library in the future.

More examples of customized plots can be found in the repository.

PyPlot Backend

To trigger Matplotlib plotting that supports MetaVLSV, using VlasiatorPyPlot. All the functions with identical names as in Matplotlib accept all possible keyword arguments supported by Matplotlib, e.g. font width, font size, colormap, etc. For detailed adjustment of plots, please refer to the Matplotlib documentation for details.

• Scalar colored contour from 2D simulation

pcolormesh(meta, "n")

• Vector z-component colored contour from 2D simulation in a manually set range

pcolormesh(meta, "n", comp=:z, colorscale=Log, axisunit=EARTH, vmin=1e6, vmax=2e6)

• Vz colored contour from 2D simulation with prescribed colormap

pcolormesh(meta, "proton/vg_v", comp=:z, colorscale=Linear, cmap=matplotlib.cm.RdBu_r)

• Derived quantity colored contour from 2D simulation (as long as the input variable is in the predefined dictionary)

pcolormesh(meta, "Bmag", comp=:z, colorscale=Linear, axisunit=SI)

• Streamline from 2D simulation

streamplot(meta, "proton/vg_v", comp="xy")

• Quiver from 2D simulation

quiver(meta, "proton/vg_v", comp="xy")

The comp option is used to specify the two vector components.

You can choose to use a linear/log/symlog color scale by setting keyword colorscale to Linear, Log, or SymLog, plot vector components by setting keyword op to :x, :y, :z, 1, 2, 3, 0 or :mag, and set axisunit to EARTH or SI etc.

• Mesh denoted by cell centers

plotmesh(meta; projection="z", color="k")

• Cut slice colored contour from 3D simulation

meta = load("bulk.amr.vlsv")
pcolormesh(meta, "proton/vg_rho", normal=:y, origin=0.0)

• Velocity distribution slice plot near a given spatial location

meta = load("bulk.1d.vlsv")
coordinates = [0.0, 0.0, 0.0]
vdfslice(meta, coordinates)

• Quick interactive REPL-based function for data inspection

pui(meta)

Or pass filename directly like pui(file).

For a full list available optional arguments, please refer to the doc for each method

Plots Backend

To trigger Plots.jl plotting, using Plots. This backend supports all available attributes provided by Plots.jl. By default it uses GR, but many other plotting libraries are also supported.

• Scaler colored contour from 2D simulation

var = "n"
heatmap(meta, var, aspect_ratio=:equal, c=:turbo)

• Scaler colored contour with lines from 2D simulation

var = "n"
contourf(meta, var)

• VDF projected slice in a normal direction

meta = load("bulk.1d.vlsv")
location = [0.0, 0.0, 0.0]
vdfslice(meta, location)

Makie Backend

Makie plotting is supported via package extension after Julia 1.9. The extension is automatically loaded when both Vlasiator and one of the Makie backends (e.g. GLMakie, CairoMakie). Before Julia 1.9, a standalone package VlasiatorMakie.jl is designed for plotting with Makie. To trigger Makie plotting with OpenGL, using VlasiatorMakie, GLMakie.

You can either use intrinsic Makie plotting methods like

lines(meta, var)   # 1D
heatmap(meta, var) # 2D

or use full recipes for quick inspection of data

• 2D slices of 3D AMR data

vlslice(meta, var; normal=:x)

• Orthognal slices of 3D AMR data

vlslices(meta, var)

• 2D slice of VDFs at a spatial cell

fig, ax = vdfslice(meta, location)
fig

• Orthognal slices of VDFs at a spatial cell

vdfslices(meta, location)

• 3D scatter of VDFs at a spatial cell

fig, ax = vdfvolume(meta, location)
fig

The interactive plots are available through the OpenGL backend of Makie GLMakie. For noninteractive high fidelity plots, we can also use the Cairo backend of Makie CairoMakie. Other options can be found at Makie Ecosystem.

using Vlasiator, GLMakie

file = "bulk.2d.vlsv"
meta = load(file)

heatmap(meta, "proton/vg_rho")

3D isosurface:

fig = volume(meta, "fg_b", EARTH, 3; algorithm=:iso, isovalue=0.0, isorange=1e-9)

Single figure contour plot:

fig = Figure(size=(700, 600), fontsize=18)
ga = fig[1,1] = GridLayout()
ax = Axis(fig[1,1],
   aspect = DataAspect(),
   title = "t = $(round(meta.time, digits=1))s",
   xlabel = L"x [$R_E$]",
   ylabel = L"y [$R_E$]"
)
hmap = heatmap!(meta, "proton/vg_rho", colormap=:turbo)
cbar = Colorbar(fig, hmap, label=L"$\rho$ [amu/cc]", width=13,
                ticksize=13, tickalign=1, height=Relative(1))
fig[1,2] = cbar
colgap!(ga, 1)

Multi-figure contour plots:

fig = Figure(size=(1100, 800), fontsize=18)

axes = []
v_str = ["CellID", "proton/vg_rho", "proton/vg_v",
   "vg_pressure", "vg_b_vol", "vg_e_vol"]
c_str = ["", L"$\rho$ [amu/cc]", "[m/s]", "[Pa]", "[T]", "[V/m]"]
c = 1

for i in 1:2, j in 1:2:5
   ax = Axis(fig[i,j], aspect=DataAspect(),
      xgridvisible=false, ygridvisible=false,
      title = v_str[c],
      xlabel = L"x [$R_E$]",
      ylabel = L"y [$R_E$]")
   hmap = heatmap!(meta, v_str[c], colormap=:turbo)
   cbar = Colorbar(fig, hmap, label=c_str[c], width=13,
                ticksize=13, tickalign=1, height=Relative(1))
   fig[i, j+1] = cbar
   c += 1
   push!(axes, ax) # just in case you need them later.
end

fig[0, :] = Label(fig, "t = $(round(meta.time, digits=1))s")

Adjusting axis limits:

location = [0, 0, 0]
fig = vdfslice(meta, location)
xlims!(fig.content[1], -1000, 1000)
ylims!(fig.content[1], -1000, 1000)
limits!(fig.content[1], 0, 10, 0, 10) # xmin, xmax, ymin, ymax

Saving figure:

fig = vdfvolume(meta, location)
save("output.png", fig)

The resolution is a property of the Figure object returned from the function.

Appending to VLSV

We are able to compute derived quantities from an original VLSV file and generate a new VLSV output with new quantities included.

meta = load("bulk.1d.vlsv")
vmag = readvariable(meta, "Vmag")
pa = readvariable(meta, "Panisotropy")
vars = Vector{Tuple{VecOrMat, String, VarInfo}}(undef, 0)
using LaTeXStrings
push!(vars, (vmag, "vmag", VarInfo("m/s", L"$\mathrm{m}/mathrm{s}$", L"$V$", "")))
push!(vars, (pa, "panisotropy", VarInfo("", "", "", "")))

write_vlsv("bulk.1d.vlsv", "bulk.1d_new.vlsv", vars)

Converting to VTK

We can convert VLSV files into VTK format. Since DCCRG is Cartesian based with uniform spacing, each level of mesh refinement corresponds to a VTK image file, and the cell refinement relationships are stored in vtkGhostType as well as the vthb file.

To convert a VLSV file into VTK,

write_vtk(file)

This function accepts either string of file names or MetaVLSV.

To see the full list of options, please refer to the documentation in API Reference. Example usage can be found here.

Tracking log files

The runtime performance per iteration can be monitored through log files:

file = "logfile.txt"
timestamps, speed = readlog(file)

Here is a live demo.

Examples

More examples are provided about

• Plotting with PyPlot
• Plotting with Plots
• Extracing variable along a line
• Field line tracing
• Simulation log file tracking
• Converting VLSV to VTK format
• Parallel post-processing
• Finding X-points and O-points in 2D reconnections

Feel free to check those out and try on your data!