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.
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 constant | Value | Meaning |
|---|
| qₑ | -1.60217662e-19 | electron charge, [C] |
| mₑ | 9.10938356e-31 | electron mass, [kg] |
| qᵢ | 1.60217662e-19 | proton mass, [C] |
| mᵢ | 1.673557546e-27 | proton mass, [kg] |
| c | 299792458. | speed of light, [m/s] |
| μ₀ | 4π*1e-7 | Vacuum permeability, [H/m] |
| ϵ₀ | 1/(c^2*μ₀) | Vacuum permittivity, [F/m] |
| kB | 1.38064852e-23 | Boltzmann constant, [m²kg/(s²K)] |
| RE | 6.371e6 | Earth radius, [m] |
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.
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.
For convenience we support the do-block syntax that automatically closes the file stream.
t = load(file) do meta
+ readparameter(meta, "time")
+end
readvariablemeta(meta, "proton/vg_rho")
A list of utility functions has been implemented for checking variable status. See here for the full list.
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.
One may want to check if two vlsv files are identical. This is tricky because
- the structure of VLSV format does not guarantee the writing order in parallel processing;
- 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.
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:
| Derived variable name | Meaning | Required variable |
|---|
| Bmag | magnetic field magnitude | vg_b_vol |
| Emag | electric field magnitude | vg_e_vol |
| Vmag | bulk speed | vg_v |
| Bhat | unit magnetic field | vg_b_vol |
| VS | sound speed | vg_ptensor_diagonal; vg_rho |
| VA | Alfvén speed | vg_rho; Bmag |
| MA | Alfvén Mach number | Vmag; VA |
| MS | Sonic Mach number | Vmag; VS |
| Epar | $\mathbf{E}_\parallel$ | vg_e_vol; Bhat |
| Eperp | $\mathbf{E}_\perp$ | vg_e_vol; Bhat |
| Vpar | bulk velocity $\parallel\mathbf{B}$ | vg_v; vg_b_vol |
| Vperp | bulk velocity $\perp \mathbf{B}$ | vg_v; vg_b_vol |
| Vth | proton thermal velocity | P; vg_rho |
| P | scalar thermal pressure | vg_ptensor_diagonal |
| Ppar | pressure $\parallel\mathbf{B}$ | vg_ptensor_diagonal; vg_b_vol |
| Pperp | pressure $\perp \mathbf{B}$ | vg_ptensor_offdiagonal; vg_b_vol |
| T | scalar temperature | P; vg_rho |
| Tpar | temperature $\parallel\mathbf{B}$ | vg_rho; vg_ptensor_diagonal; vg_b_vol |
| Tperp | temperature $\perp \mathbf{B}$ | vg_rho; vg_ptensor_offdiagonal; vg_b_vol |
| Tanisotropy | $T_\perp / T_\parallel$ | Tpar; Tperp |
| J | current density | vg_b_vol |
| Jpar | $j_\parallel$ | vg_b_vol |
| Jperp | $j_\perp$ | vg_b_vol |
| Protated | pressure tensor with $\widehat{z} \parallel \mathbf{B}$ | vg_b_vol; vg_ptensor_diagonal; vg_ptensor_offdiagonal |
| Panisotropy | $P_\perp / P_\parallel$ | ptensor; B |
| Pram | dynamic ram pressure | vg_rho; Vmag |
| Pb | magnetic pressure | vg_b_vol |
| Poynting | Poynting flux | E; B |
| Beta | plasma $\beta$, $P / P_B$ | P; vg_b_vol |
| BetaStar | modified $\beta$, $(P+P_{ram})/P_B$ | P; Pram; vg_b_vol |
| IonInertial | proton inertial length | vg_rho |
| Larmor | proton Larmor radius | Vth; Bmag |
| Gyroperiod | proton gyroperiod | Bmag |
| PlasmaPeriod | plasma oscillation period | vg_rho |
| Gyrofrequency | proton gyro-frequency | Bmag |
| Omegap | plasma frequency (proton) | vg_rho |
| MagneticTension | magnetic tension force | vg_b_vol |
| n | proton number density | vg_rho |
which can also be found as keys of dictionary in vlsvvariables.jl.
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.
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.
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.
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.
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).
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
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 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
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.
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)
Writing new FsGrid variables is not supported. All quantities from the original file is maintained.
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.
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.
The runtime performance per iteration can be monitored through log files:
file = "logfile.txt"
+timestamps, speed = readlog(file)
Here is a live demo.
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!
