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.

Physical constantValueMeaning
qₑ-1.60217662e-19electron charge, [C]
mₑ9.10938356e-31electron mass, [kg]
qᵢ1.60217662e-19proton mass, [C]
mᵢ1.673557546e-27proton mass, [kg]
c299792458.speed of light, [m/s]
μ₀4π*1e-7Vacuum permeability, [H/m]
ϵ₀1/(c^2*μ₀)Vacuum permittivity, [F/m]
kB1.38064852e-23Boltzmann constant, [m²kg/(s²K)]
RE6.371e6Earth radius, [m]

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]:

Derived variable nameMeaningRequired variable[2]
Bmagmagnetic field magnitudevg_b_vol
Emagelectric field magnitudevg_e_vol
Vmagbulk speedvg_v
Bhatunit magnetic fieldvg_b_vol
VSsound speedvg_ptensor_diagonal; vg_rho
VAAlfvén speedvg_rho; Bmag
MAAlfvén Mach numberVmag; VA
MSSonic Mach numberVmag; VS
Epar$\mathbf{E}_\parallel$vg_e_vol; Bhat
Eperp$\mathbf{E}_\perp$vg_e_vol; Bhat
Vparbulk velocity $\parallel\mathbf{B}$vg_v; vg_b_vol
Vperpbulk velocity $\perp \mathbf{B}$vg_v; vg_b_vol
Vthproton thermal velocityP; vg_rho
Pscalar thermal pressurevg_ptensor_diagonal
Pparpressure $\parallel\mathbf{B}$vg_ptensor_diagonal; vg_b_vol
Pperppressure $\perp \mathbf{B}$vg_ptensor_offdiagonal; vg_b_vol
Tscalar temperatureP; vg_rho
Tpartemperature $\parallel\mathbf{B}$vg_rho; vg_ptensor_diagonal; vg_b_vol
Tperptemperature $\perp \mathbf{B}$vg_rho; vg_ptensor_offdiagonal; vg_b_vol
Tanisotropy$T_\perp / T_\parallel$Tpar; Tperp
Jcurrent densityvg_b_vol
Jpar$j_\parallel$vg_b_vol
Jperp$j_\perp$vg_b_vol
Protatedpressure tensor with $\widehat{z} \parallel \mathbf{B}$vg_b_vol; vg_ptensor_diagonal; vg_ptensor_offdiagonal
Panisotropy$P_\perp / P_\parallel$ptensor; B
Pramdynamic ram pressurevg_rho; Vmag
Pbmagnetic pressurevg_b_vol
PoyntingPoynting fluxE; B
Betaplasma $\beta$, $P / P_B$P; vg_b_vol
BetaStarmodified $\beta$, $(P+P_{ram})/P_B$P; Pram; vg_b_vol
IonInertialproton inertial lengthvg_rho
Larmorproton Larmor radiusVth; Bmag
Gyroperiodproton gyroperiodBmag
PlasmaPeriodplasma oscillation periodvg_rho
Gyrofrequencyproton gyro-frequencyBmag
Omegapplasma frequency (proton)vg_rho
MagneticTensionmagnetic tension forcevg_b_vol
nproton number densityvg_rho

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

Note

In Vlasiator, the cells inside the inner boundary (which is usually a sphere/circle) are filled with zero density values. This is then used to identify the inner boundary for all other quantities. Therefore, if you are manipulating directly on data, make sure that the nonsense values inside the inner boundary are excluded. One way to do this can be found 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.

Warning

The method call to certain axes is not dispatched, e.g. ax.plot; as an alternative, one needs to pass ax as the third argument to the functions, e.g. plot(meta, "rho", ax). See Matplotlib's two interfaces for the history of the interfaces.

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

Note

This is an experimental feature. We plan to have GUI-based plotting support in the future.

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

Writing new FsGrid variables is not supported. All quantities from the original file is maintained.

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.

Warning

As of ParaView 5.9.1, there are display issues with VTKOverlappingAMR. However, we can read the generated image files directly. There is also an keyword argument for write_vtk called maxamronly: when it is set to true, then only the image file at the highest refinement level is generated. This part is experimental and subject to change in the future.

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!

  • 1For species specific variables, you need to add the species name at the front, separated by a slash. For example, the proton bulk velocity is a string proton/vg_v.
  • 2If a required variable exists in the VLSV file, we try to use it directly instead of calculating from other variables. The interpolated FS grid variables onto DCCRG grid are preferred over original FS grid variables.