+ ${display_result} +
+This section contains thorough examples of using Batsrus.jl.
|
+ 1D line animation using pyplot |
+
|
+ Subplot 1D line animation using pyplot |
+
|
+ 2D contour animation using pyplot |
+
Settings
This document was generated with Documenter.jl version 1.7.0 on Tuesday 12 November 2024. Using Julia version 1.11.1.
This example shows how to create 1D space line animation from series of SWMF outputs.
using Batsrus, PyPlot, Printf
+
+"""
+ animate1d(files::Vector{String}, vars::Vector{String}; kwargs...)
+
+Saving series of plots of `vars` from SWMF output `files`.
+
+# Keywords
+- `filedir::String="./"`: input SWMF file directory.
+- `outdir::String="out/"`: output directory.
+- `output_interval::Int=1`: Timestep interval for output files.
+- `vmin=-1`: plot value lower bound.
+- `vmax=1`: plot value upper bound.
+- `refloc=1e5`: reference vertical line initial location.
+"""
+function animate1d(files::Vector{String}, vars::Vector{String};
+ filedir::String="./", outdir::String="out/", output_interval::Int=1,
+ vmin=-1, vmax=1, refloc=1e5)
+ nfile = length(files)
+ x = let
+ bd = load(joinpath(filedir, files[1]))
+ bd.x[:,1,1]
+ end
+
+ fig = plt.figure(figsize=(12,5), constrained_layout=true)
+ ax = plt.axes(xlim=extrema(x), ylim=(vmin, vmax))
+
+ line, = ax.plot([], [], lw=1)
+ # Move reference line forward in the canvas
+ vl = ax.axvline(0, ls="-", color="tab:brown", lw=1, zorder=10)
+
+ color = "tab:blue"
+ ax.set_xlabel("x [km]"; fontsize=14)
+ ax.set_ylabel(var*" [nT]"; fontsize=14, color)
+ ax.tick_params(axis="y", labelcolor=color)
+
+ ax2 = ax.twinx()
+ ax2.set_ylim(-600, 0)
+
+ color = "tab:red"
+ ax2.set_ylabel("Ux [km/s]"; fontsize=14, color)
+
+ line2, = ax2.plot([], []; color, alpha=0.7)
+ ax2.tick_params(axis="y", labelcolor=color)
+
+
+ for (i, file) in enumerate(files)
+ @info "$i in $nfile"
+ outname = outdir*lpad(i, 4, '0')*".png"
+ isfile(outname) && continue
+
+ bd = load(joinpath(filedir, file))
+
+ d = bd[vars[1]][:,1]
+ title_str = @sprintf "t = %4.1f s" bd.head.time
+ line.set_data(x, d)
+ d2 = bd[vars[2]][:,1]
+ line2.set_data(x, d2)
+
+ ax.set_title(title_str)
+
+ refloc -= output_interval * VSW
+ if refloc <= 0
+ refloc += 1e5
+ end
+ vl.set_xdata([refloc, refloc])
+
+ savefig(outname, bbox_inches="tight", dpi=200)
+ end
+
+ close()
+end
+
+####################
+# Constants
+const VSW = 500.0 # Solar wind speed, [km/s]
+
+# Data directory
+filedir = "./"
+# Plot variables
+vars = ["By", "uxs1"]
+
+# Find simulation data
+files = let
+ pick_file = file -> startswith(file, "z") && endswith(file, ".out")
+ filter(pick_file, readdir(filedir))
+end
+
+animate1d(files, vars)This page was generated using DemoCards.jl.
Settings
This document was generated with Documenter.jl version 1.7.0 on Tuesday 12 November 2024. Using Julia version 1.11.1.
This example shows how to create multi-panel 1D space line animation from a series of FLEKS outputs.
Since currently we don't have true 1D outputs, in this demo we show the way to extract the first column from 2D cuts.
using Batsrus, Printf, PyPlot
+
+function create_figure()
+ fig, axs = subplots(5, 1, figsize=(12, 8),
+ sharex=true, sharey=false, constrained_layout=true)
+
+ lw1 = 1.0
+ lw2 = 1.0
+ lw3 = 1.0
+
+ xmin, xmax = -6300.0, -1300.0
+
+ ρemin, ρemax = 0.0, 0.05 # mi/me = 400
+ ρimin, ρimax = 0.0, 20.0
+ vmin, vmax = -600.0, 300.0
+ pmin, pmax = 1e-3, 1.0
+ Bmin, Bmax = -3.0, 30.0
+ Emin, Emax = -5.0, 5.0
+
+ l11, = axs[1].plot([], [], lw=lw1)
+ l211, = axs[2].plot([], [], lw=lw1, label="Uex")
+ l212, = axs[2].plot([], [], lw=lw1, label="Uey")
+ l213, = axs[2].plot([], [], lw=lw1, label="Uez")
+ l311, = axs[3].plot([0.0, 1.0], [pmin, pmax], lw=lw1, alpha=0.8, label=L"$P_{exx}$")
+ l312, = axs[3].plot([0.0, 1.0], [pmin, pmax], lw=lw1, alpha=0.8, label=L"$P_{eyy}$")
+ l313, = axs[3].plot([0.0, 1.0], [pmin, pmax], lw=lw1, alpha=0.8, label=L"$P_{ezz}$")
+ l41, = axs[4].plot([], [], lw=lw2, label="x")
+ l42, = axs[4].plot([], [], lw=lw2, label="y")
+ l43, = axs[4].plot([], [], lw=lw2, label="z")
+ l51, = axs[5].plot([], [], lw=lw2, label="x")
+ l52, = axs[5].plot([], [], lw=lw2, label="y")
+ l53, = axs[5].plot([], [], lw=lw2, label="z")
+
+ axs[1].set_xlim(xmin, xmax)
+ axs[1].set_ylim(ρemin, ρemax)
+ axs[2].set_ylim(vmin, vmax)
+ axs[3].set_yscale("log")
+ axs[3].set_ylim(pmin, pmax)
+ axs[4].set_ylim(Bmin, Bmax) # B
+ axs[5].set_ylim(Emin, Emax) # E
+
+ for ax in axs
+ ax.grid(true)
+ end
+
+ axs[4].set_ylabel("B [nT]"; fontsize)
+ axs[5].set_ylabel("E [mV/m]"; fontsize)
+
+ axs[end].set_xlabel("x [km]"; fontsize)
+
+ axs[1].set_ylabel(L"$\rho_e$ [amu/cc]"; fontsize, color="tab:blue")
+ axs[1].tick_params(axis="y", labelcolor="tab:blue")
+
+ ax12 = axs[1].twinx()
+ ax12.set_ylim(ρimin, ρimax)
+
+ ax12.set_ylabel(L"$\rho_i$ [amu/cc]"; fontsize, color="tab:red")
+
+ l12, = ax12.plot([], []; lw=lw1, color="tab:red", alpha=0.8)
+ ax12.tick_params(axis="y", labelcolor="tab:red")
+
+ axs[2].set_ylabel(L"$U_e$ [km/s]"; fontsize)
+
+ ax22 = axs[2].twinx()
+ ax22.set_ylim(vmin, vmax)
+
+ ax22.set_ylabel(L"$U_i$ [km/s]"; fontsize)
+
+ l221, = ax22.plot([], [], lw=lw3, label="Uix", color="tab:red")
+ l222, = ax22.plot([], [], lw=lw3, label="Uiy", color="tab:purple")
+ l223, = ax22.plot([], [], lw=lw3, label="Uiz", color="tab:brown")
+
+ axs[3].set_ylabel(L"$P_e$ [nT]"; fontsize)
+
+ ax32 = axs[3].twinx()
+ ax32.set_ylabel(L"$P_i$ [nT]"; fontsize)
+
+ fake_range = [0.0, 1.0]
+ p_range = [pmin, pmax]
+ p_alpha = 0.8
+ l321, = ax32.plot(fake_range, p_range, lw=lw3, alpha=p_alpha, label=L"$P_{ixx}$",
+ color="tab:red")
+ l322, = ax32.plot(fake_range, p_range, lw=lw3, alpha=p_alpha, label=L"$P_{iyy}$",
+ color="tab:purple")
+ l323, = ax32.plot(fake_range, p_range, lw=lw3, alpha=p_alpha, label=L"$P_{izz}$",
+ color="tab:brown")
+
+ ax32.set_yscale("log")
+ ax32.set_ylim(p_range...)
+
+ leg21 = axs[2].legend(;loc=(0.34, 0.05), ncols=3, frameon=false)
+ leg22 = ax22.legend(;loc=(0.0, 0.05), ncols=3, frameon=false)
+ leg31 = axs[3].legend(;loc=(0.34, -0.05), ncols=3, frameon=false)
+ leg32 = ax32.legend(;loc=(0.0, -0.05), ncols=3, frameon=false)
+ leg4 = axs[4].legend(;loc="upper left", ncols=3, frameon=false)
+ leg5 = axs[5].legend(;loc="upper left", ncols=3, frameon=false)
+
+ # set the linewidth of each legend object
+ legs = (leg21, leg22, leg31, leg32, leg4, leg5)
+ for leg in legs
+ for legobj in leg.legend_handles
+ legobj.set_linewidth(1.5)
+ end
+ end
+
+ ls = (l11, l12, l211, l212, l213, l221, l222, l223, l311, l312, l313, l321, l322, l323,
+ l41, l42, l43, l51, l52, l53)
+
+ axs, ls
+end
+
+function slice1d_avg(bd, var, dir::Int=2)
+ mean(bd[var], dims=dir) |> vec
+end
+
+function animate(files::Vector{String}, axs, ls; outdir="out/", overwrite::Bool=false,
+ icut::Int=1, nbox::Int=1, dir=2, doAverage::Bool=false)
+ l11, l12, l211, l212, l213, l221, l222, l223, l311, l312, l313, l321, l322, l323,
+ l41, l42, l43, l51, l52, l53 = ls
+
+ for (i, file) in enumerate(files)
+ @info "$i in $(length(files))"
+ outname = outdir*lpad(i, 4, '0')*".png"
+ if !overwrite
+ isfile(outname) && continue
+ end
+
+ bd = load(joinpath(filedir, file))
+
+ x = @views bd.x[:,1,1]
+
+ if doAverage
+ slice = let bd = bd, dir = dir, nbox = nbox
+ var -> moving_average(slice1d_avg(bd, var, dir), nbox)
+ end
+ else
+ slice = let bd = bd, dir = dir, nbox = nbox, icut = icut
+ var -> moving_average(slice1d(bd, var, icut, dir), nbox)
+ end
+ end
+
+ d = slice("rhos0")
+ l11.set_data(x, d)
+ d = slice("rhos1")
+ l12.set_data(x, d)
+ d = slice("Uxs0")
+ l211.set_data(x, d)
+ d = slice("Uys0")
+ l212.set_data(x, d)
+ d = slice("Uzs0")
+ l213.set_data(x, d)
+ d = slice("Uxs1")
+ l221.set_data(x, d)
+ d = slice("Uys1")
+ l222.set_data(x, d)
+ d = slice("Uzs1")
+ l223.set_data(x, d)
+ d = slice("PXXS0")
+ l311.set_data(x, d)
+ d = slice("PYYS0")
+ l312.set_data(x, d)
+ d = slice("PZZS0")
+ l313.set_data(x, d)
+ d = slice("PXXS1")
+ l321.set_data(x, d)
+ d = slice("PYYS1")
+ l322.set_data(x, d)
+ d = slice("PZZS1")
+ l323.set_data(x, d)
+
+ d = slice_data("Bx")
+ l41.set_data(x, d)
+ d = slice_data("By")
+ l42.set_data(x, d)
+ d = slice_data("Bz")
+ l43.set_data(x, d)
+ d = slice_data("Ex") ./ 1000
+ l51.set_data(x, d)
+ d = slice_data("Ey") ./ 1000
+ l52.set_data(x, d)
+ d = slice_data("Ez") ./ 1000
+ l53.set_data(x, d)
+
+ title_str = @sprintf "t = %4.1f s" bd.head.time
+ axs[1].set_title(title_str; fontsize)
+
+ savefig(outname, bbox_inches="tight", dpi=200)
+ end
+
+ return
+end
+
+#################
+
+# Data directory
+filedir = "./"
+
+const fontsize = 16
+
+pick_file = file -> startswith(file, "z") && endswith(file, ".out")
+
+files = filter(pick_file, readdir(filedir))
+
+axs, ls = create_figure()
+animate(files, axs, ls; icol=1, outdir="figures/", overwrite=true)
+
+close()This page was generated using DemoCards.jl.
Settings
This document was generated with Documenter.jl version 1.7.0 on Tuesday 12 November 2024. Using Julia version 1.11.1.
This example shows how to create 2D colored contour animation from series of SWMF outputs using Matplotlib.
using Batsrus, PyPlot, Printf, PyCall
+
+const fontsize = 14
+
+"""
+ animate(files::Vector{String}; kwargs...)
+
+Generate figures of colored contour, optionally with streamlines.
+
+# Keywords
+- `var`: variable to plot with pcolormesh.
+- `vmin`: minimum plotting value.
+- `vmax`: maximum plotting value.
+- `plotrange`: 2D plotting spatial range.
+- `plotinterval`: spatial sampling interval.
+- `cmap`: colormap accepted by Matplotlib.
+- `streamvars`: if set, add streamlines to the colored contour.
+- `density`: streamline density.
+- `orientation`: colorbar orientation, "horizontal" or "vertical".
+- `extend`: colorbar extension, "neither", "both", "min", "max".
+- `use_conventional_notation`: if true, denote the vertical axis as "z".
+- `broken_streamlines`: if true, streamlines can be broken.
+- `overwrite`: if true, overwrite the existing image files.
+- `filedir`: source data directory.
+- `outdir`: output image directory.
+"""
+function animate(files::Vector{String}; var="Uy", vmin=-Inf, vmax=Inf,
+ plotrange=[-Inf, Inf, -Inf, Inf], plotinterval=Inf,
+ cmap=matplotlib.cm.RdBu_r, streamvars=nothing, density=1, orientation="horizontal",
+ extend="neither", use_conventional_notation=false, broken_streamlines=false,
+ overwrite=false, filedir="./", outdir="out/")
+
+ bd = load(filedir*files[1])
+
+ nfile = length(files)
+ if isnan(vmin+vmax)
+ vmin = isinf(vmin) ? minimum(v) : vmin
+ vmax = isinf(vmax) ? maximum(v) : vmax
+ end
+
+ norm = matplotlib.colors.Normalize(vmin, vmax)
+
+ filetype =
+ if occursin("region", files[1])
+ :PC
+ else
+ :GM
+ end
+
+ if use_conventional_notation # switch y and z axis
+ if var[2] == 'y'
+ cb_str =
+ if startswith(var, "U")
+ "Uz [km/s]"
+ elseif startswith(var, "B")
+ "Bz [nT]"
+ elseif startswith(var, "E")
+ L"Ez [$\mu\mathrm{V}/\mathrm{m}$]"
+ else
+ var
+ end
+ elseif var[2] == 'z'
+ cb_str =
+ if startswith(var, "U")
+ "Uy [km/s]"
+ elseif startswith(var, "B")
+ "By [nT]"
+ elseif startswith(var, "E")
+ L"Ey [$\mu\mathrm{V}/\mathrm{m}$]"
+ else
+ var
+ end
+ end
+ else
+ cb_str =
+ if startswith(var, "U")
+ var * " [km/s]"
+ elseif startswith(var, "B")
+ var * " [nT]"
+ elseif startswith(var, "E")
+ var * L" [\mu\mathrm{V}/\mathrm{m}]"
+ else
+ var
+ end
+ end
+
+ fig = plt.figure(figsize=(4, 8), constrained_layout=true)
+ ax = plt.axes()
+ ax.set_xlabel(L"x [$\mathrm{R}_\mathrm{E}$]"; fontsize)
+ ax.set_ylabel(L"y [$\mathrm{R}_\mathrm{E}$]"; fontsize)
+
+ @info "1 out of $nfile"
+ x, y, w = Batsrus.getdata2d(bd, var, plotrange, plotinterval; innermask=false)
+ c = ax.pcolormesh(x, y, w; norm, cmap)
+ ax.set_aspect("equal", adjustable="box", anchor="C")
+
+ if orientation == "horizontal"
+ cb = colorbar(c; ax, orientation, location="top", aspect=50, extend)
+ else
+ cb = colorbar(c; ax, orientation, extend)
+ end
+ cb.ax.set_ylabel(cb_str; fontsize)
+ #cb.ax.set_title(cb_str; fontsize)
+ title_str = @sprintf "t = %4.1f s" bd.head.time
+ ax.set_title(title_str)
+
+ if bd.head.gencoord
+ if !isnothing(streamvars)
+ xi, yi, v1, v2 = get_vector(bd, streamvars, plotrange, plotinterval)
+ st = ax.streamplot(xi, yi, v1, v2; color="gray", density)
+ end
+ else
+ if !isnothing(streamvars)
+ st = streamplot(bd, streamvars, ax; color="gray", density, broken_streamlines)
+ end
+ end
+
+ savefig(outdir*lpad(1, 4, '0')*".png", bbox_inches="tight", dpi=200)
+
+ if !isnothing(streamvars) # Clean up streamlines
+ st.lines.remove()
+ for art in ax.get_children()
+ if !pybuiltin(:isinstance)(art, matplotlib.patches.FancyArrowPatch)
+ continue
+ end
+ art.remove()
+ end
+ end
+
+ for i in 2:nfile
+ @info "$i out of $nfile"
+ outname = outdir*lpad(i, 4, '0')*".png"
+ if !overwrite
+ isfile(outname) && continue
+ end
+ bd = load(joinpath(filedir, files[i]))
+ if bd.head.gencoord
+ _, _, wi = getdata2d(bd, var, plotrange, plotinterval; innermask=false)
+ else
+ wi = bd[var]'
+ end
+
+ c.set_array(wi)
+
+ if !isnothing(streamvars)
+ if bd.head.gencoord
+ xi, yi, v1, v2 = get_vector(bd, streamvars, plotrange, plotinterval)
+ st = ax.streamplot(xi, yi, v1, v2; color="gray", density)
+ else
+ st = streamplot(bd, "Bx;By", ax; color="gray", density, broken_streamlines)
+ end
+ end
+
+ title_str = @sprintf "t = %4.1f s" bd.head.time
+ ax.set_title(title_str)
+
+ savefig(outname, bbox_inches="tight", dpi=200)
+ if !isnothing(streamvars) # Clean up streamlines
+ st.lines.remove()
+ for art in ax.get_children()
+ if !pybuiltin(:isinstance)(art, matplotlib.patches.FancyArrowPatch)
+ continue
+ end
+ art.remove()
+ end
+ end
+ end
+
+ close()
+end
+
+function get_vector(bd, var, plotrange, plotinterval)
+ x, w = bd.x, bd.w
+ X, Y = eachslice(x, dims=3)
+ X, Y = vec(X), vec(Y)
+ varstream = split(var, ";")
+ var1_ = findfirst(x->lowercase(x)==lowercase(varstream[1]), bd.head.wnames)
+ var2_ = findfirst(x->lowercase(x)==lowercase(varstream[2]), bd.head.wnames)
+ # Create grid values first.
+ xi = range(Float64(plotrange[1]), stop=Float64(plotrange[2]), step=plotinterval)
+ yi = range(Float64(plotrange[3]), stop=Float64(plotrange[4]), step=plotinterval)
+
+ tr = matplotlib.tri.Triangulation(X, Y)
+ Xi, Yi = Batsrus.meshgrid(xi, yi)
+
+ W1 = @views w[:,:,var1_] |> vec
+ interpolator = matplotlib.tri.LinearTriInterpolator(tr, W1)
+ v1 = interpolator(Xi, Yi)
+
+ W2 = @views w[:,:,var2_] |> vec
+ interpolator = matplotlib.tri.LinearTriInterpolator(tr, W2)
+ v2 = interpolator(Xi, Yi)
+
+ xi, yi, v1, v2
+end
+
+#######################
+filedir = "GM/"
+
+files = filter(file -> startswith(file, "z") && endswith(file, ".out"), readdir(filedir))
+
+var = "Bz"
+vmin = -10
+vmax = 10
+cmap = matplotlib.cm.RdBu_r
+streamvars = "Bx;By"
+density = 1
+orientation = "vertical"
+outdir = "out/"
+extend = "both"
+use_conventional_notation = true
+
+animate(files; filedir, outdir, var, vmin, vmax, orientation, cmap, streamvars, density,
+ extend, use_conventional_notation)This page was generated using DemoCards.jl.
Settings
This document was generated with Documenter.jl version 1.7.0 on Tuesday 12 November 2024. Using Julia version 1.11.1.
This package is still under development, so be careful for any future breaking changes!
BATSRUS and SWMF data reading, converting, visualizing and analyzing in Julia.
This package provides the following functionalities:
The ultimate goal is to build a convenient tool of reading and analyzing simulation outputs which is easy to install and easy to use.
Feel free to contact the author for any help or collaboration!
Install VisAna from the julia REPL prompt with
using Pkg
+Pkg.add("Batsrus")Or in the Pkg REPL
julia> ]
+pkg> add BatsrusData loading speed of a 2.4GB 3D binary file, 317MB 3D binary file, and 65KB 2D binary file on Macbook Pro with quad core 2.2 GHz Intel i7 and 16 GB 1600 MHz DDR3:
| 2.4GB | tmax [s] | tmean [s] |
|---|---|---|
| Julia | 2.73 | 1.32 |
| Python | 1.35 | 1.34 |
| IDL | 6.18 | 6.08 |
| MATLAB | 16.02 | 10.60 |
| 317MB | tmean [ms] |
|---|---|
| Julia | 180.8 |
| Python | 179.5 |
| IDL | 453.5 |
| MATLAB | 698.4 |
| 65KB | tmean [μs] |
|---|---|
| Julia | 163.36 |
| Python | 4390.95 |
| IDL | 1970.29 |
| MATLAB | 19273.25 |
The Julia, IDL, and MATLAB version all shares the same kernel design. The timings are obtained for Julia v1.3.1, Python 3.7.6 + Numpy 1.18.1, IDL 8.5, and MATLAB R2018b. For dynamic languages with JIT, the first time when function gets executed is also the slowest due to runtime compilation, as can be seen from tmax in the tables. spacepy reaches the same level of performance as Batsruls.jl because of the well-optimized numpy library written in C. However, for small data sizes Batsrus.jl is much faster than packages written in other languages.
In Python, you can easily take advantage of this package with the aid of JuliaCall or PyJulia.
With JuliaCall:
from juliacall import Main as jl
+jl.seval("using Batsrus")
+file = 'test/example.out'
+data = Batsrus.load(file)With PyJulia:
from julia import Batsrus
+file = 'test/1d__raw_2_t25.60000_n00000258.out'
+data = Batsrus.load(file)PyPlot package backend may be affected by the settings of PyJulia dependencies. If you want to set it back properly, you need to recompile the PyCall package in Julia.
This package inherits the ideas and code structures from its predecessor in IDL (developed by Gábor Tóth) and MATLAB.
Batsrus.jl is developed and maintained by Hongyang Zhou.
Settings
This document was generated with Documenter.jl version 1.7.0 on Tuesday 12 November 2024. Using Julia version 1.11.1.
Batsrus.BATS — TypeBatsrus data container. Dim is the dimension of output, and Dimp1 is an extra parameter for working around derived types.
Batsrus.Batl — TypeBATSRUS output high-level struct.
Batsrus.BatsHead — TypeBatsrus file head information.
Batsrus.FileList — TypeType for Batsrus file information.
Batsrus.Head — TypeBATSRUS output standalone header information.
Batsrus._fill_vector_from_scalars — MethodConstruct vectors from scalar components.
Batsrus.adjust_plotrange! — MethodAdjust 2D plot ranges.
Batsrus.allocateBuffer — MethodCreate buffer for x and w.
Batsrus.convertIDLtoVTK — MethodconvertIDLtoVTK(filename; gridType=1, verbose=false)Convert 3D BATSRUS *.out to VTK. If gridType==1, it converts to the rectilinear grid; if gridType==2, it converts to the structured grid. If filename does not end with "out", it tries to find the ".info" and ".tree" file with the same name tag and generates 3D unstructured VTU file.
Batsrus.convertTECtoVTU — FunctionconvertTECtoVTU(file::AbstractString, outname="out")
+convertTECtoVTU(head, data, connectivity, outname="out")Convert unstructured Tecplot data to VTK. Note that if using voxel type data in VTK, the connectivity sequence is different from Tecplot: the 3D connectivity sequence in Tecplot is the same as the hexahedron type in VTK, but different with the voxel type. The 2D connectivity sequence is the same as the quad type, but different with the pixel type. For example, in 3D the index conversion is:
# PLT to VTK voxel index_ = [1 2 4 3 5 6 8 7]
+for i in 1:2
+ connectivity = swaprows!(connectivity, 4*i-1, 4*i)
+endBatsrus.create_pvd — Methodcreate_pvd(filepattern)Generate PVD file for a time series collection of VTK data.
Example
create_pvd("*.vtu)Batsrus.cutdata — Methodcutdata(data, var; plotrange=[-Inf,Inf,-Inf,Inf], dir="x", sequence=1)Get 2D plane cut in orientation dir for var out of 3D box data within plotrange. The returned 2D data lies in the sequence plane from - to + in dir.
Batsrus.fillCellNeighbors! — MethodfillCellNeighbors!(batl, iCell_G, DiLevelNei_III, iNodeNei_III, nBlock_P)Fill neighbor cell indexes for the given block. The faces, edges, and vertices are ordered from left (-) to right (+) in x-y-z sequentially.
Vertices: Edges: (10,11 ignored)
7 ––- 8 . –4– .
5 ––- 6 . . –3– . 12 . . . . . . . . . 3 ––- 4 9 . –2– . . - . - . 5 . 6 1 ––- 2 . –1– .
Only tested for 3D.
Batsrus.find_grid_block — Methodfind_grid_block(batl, xyz_D)Return processor local block index that contains a point. Input location should be given in Cartesian coordinates.
Batsrus.find_neighbor_for_anynode — MethodFind neighbors for any node in the tree. Only for Cartesian now.
Batsrus.find_tree_node — Methodfind_tree_node(batl, Coord_D)Find the node that contains a point. The point coordinates should be given in generalized coordinates normalized by the domain size.
Batsrus.findindex — MethodFind variable index in the BATSRUS data.
Batsrus.getConnectivity — MethodGet cell connectivity list.
Batsrus.getRotationMatrix — MethodgetRotationMatrix(axis::AbstractVector, angle::Real) --> SMatrix{3,3}Create a rotation matrix for rotating a 3D vector around a unit axis by an angle in radians. Reference: Rotation matrix from axis and angle
Example
using LinearAlgebra
+v = [-0.5, 1.0, 1.0]
+v̂ = normalize(v)
+θ = deg2rad(-74)
+R = getRotationMatrix(v̂, θ)Batsrus.getSibling — MethodReturn sibling index (1-8) for the given block node matrix.
Batsrus.get_var_range — MethodReturn value range of var in bd.
Batsrus.getascii! — MethodRead ascii format coordinates and data values.
Batsrus.getbinary! — MethodRead binary format coordinates and data values.
Batsrus.getfilehead — Methodgetfilehead(fileID::IoStream, filelist::FileList) -> NameTupleObtain the header information from BATSRUS output file of type linked to fileID.
Input arguments
fileID::IOStream: file identifier.filelist::FileList: file information.Batsrus.getfilesize — Methodgetfilesize(fileID::IOStream, lenstr::Int32, ::Val{FileType})Return the size in bytes for one snapshot.
Batsrus.getfiletype — MethodObtain file type.
Batsrus.getvar — Methodgetvar(bd::BATS, var::AbstractString) -> ArrayReturn variable data from string var. This is also supported via direct indexing,
Examples
bd["rho"]Batsrus.getview — MethodReturn view of variable var in bd.
Batsrus.ibits — MethodLogical shifts as the Fortran instrinsic function.
Batsrus.interp1d — Methodinterp1d(bd::BATS, var::AbstractString, loc::AbstractVector{<:AbstractFloat})Interpolate var at spatial point loc in bd.
Batsrus.interp1d — Methodinterp1d(bd::BATS, var::AbstractString, point1::Vector, point2::Vecto)Interpolate var along a line from point1 to point2 in bd.
Batsrus.interp2d — Methodinterp2d(bd::BATS, var::AbstractString, plotrange=[-Inf, Inf, -Inf, Inf],
+ plotinterval=Inf; kwargs...)Return 2D interpolated slices of data var from bd. If plotrange is not set, output data resolution is the same as the original.
Keyword Arguments
innermask=false: Whether to mask the inner boundary with NaN.rbody=1.0: Radius of the inner mask. Used when the rbody parameter is not found in the header.useMatplotlib=true: Whether to Matplotlib (somehow faster) or NaturalNeighbours for scattered interpolation.Batsrus.interpolate2d_generalized_coords — MethodPerform Triangle interpolation of 2D data W on grid X, Y.
Batsrus.load — Methodload(filename; npict=1, verbose=false)Read BATSRUS output files. Stores the npict snapshot from an ascii or binary data file into the arrays of coordinates x and data w.
Batsrus.meshgrid — MethodGenerating consistent 2D arrays for passing to plotting functions.
Batsrus.nodeToGlobalBlock — MethodReturn global block index for the node.
Batsrus.order_children! — Methodorder_children!(batl::Batl, iNode, iMorton::Int, iNodeMorton_I::Vector{Int32})Recursively apply Morton ordering for nodes below a root block. Store result into iNodeMortonI and iMortonNodeA using the iMorton index.
Batsrus.order_tree — Methodorder_tree(batl)Return maximum AMR level in the used block and the Morton curve order. Set iNodeMorton_I indirect index arrays according to
Batsrus.readhead — MethodReturn header from info file. Currently only designed for 2D and 3D.
Batsrus.readlogdata — MethodRead information from log file.
Batsrus.readtecdata — Methodreadtecdata(file; verbose=false)Return header, data and connectivity from BATSRUS Tecplot outputs. Both 2D and 3D binary and ASCII formats are supported.
Examples
file = "3d_ascii.dat"
+head, data, connectivity = readtecdata(file)Batsrus.readtree — MethodReturn BATL tree structure.
Batsrus.rotateTensorToVectorZ — MethodrotateTensorToVectorZ(tensor::AbstractMatrix, v::AbstractVector) -> SMatrix{3,3}Rotate tensor with a rotation matrix that aligns the 3rd direction with vector, which is equivalent to change the basis from (i,j,k) to (i′,j′,k′) where k′ ∥ vector. Reference: Tensor rotation
Batsrus.setunits — Methodsetunits(head, type; distance=1.0, mp=1.0, me=1.0)Set the units for the output files. If type is given as "SI", "CGS", "NORMALIZED", "PIC", "PLANETARY", "SOLAR", set typeunit = type, otherwise try to guess from the fileheader. Based on typeunit set units for distance [xSI], time [tSI], density [rhoSI], pressure [pSI], magnetic field [bSI] and current density [jSI] in SI units. Distance unit [rplanet | rstar], ion and electron mass in [amu] can be set with optional distance, mp and me.
Also calculate convenient constants ti0, cs0 ... for typical formulas. This function is currently not used anywhere!
Batsrus.showhead — Functionshowhead(file, head)Displaying file header information.
Batsrus.showhead — Methodshowhead(data)
+showhead(io, data)Display file information of data.
Batsrus.slice1d — Functionslice1d(bd, var, icut::Int=1, dir::Int=2)Return view of variable var in bd along 1D slice. icut is the index along axis dir. dir == 1 means align with the 2nd (y) axis, dir == 2 means align with the 1st (x) axis.
Batsrus.squeeze — MethodSqueeze singleton dimensions for an array A.
Batsrus.subsurface — Methodsubsurface(x, y, data, limits)
+subsurface(x, y, u, v, limits)Extract subset of 2D surface dataset in ndgrid format. See also: subvolume.
Batsrus.subvolume — Methodsubvolume(x, y, z, data, limits)
+subvolume(x, y, z, u, v, w, limits)Extract subset of 3D dataset in ndgrid format. See also: subsurface.
Batsrus.swaprows! — MethodReturn matrix X with swapped rows i and j.
Batsrus.UnitfulBatsrus.@bu_str — Macromacro bu_str(unit)String macro to easily recall Batsrus units located in the UnitfulBatsrus package. Although all unit symbols in that package are suffixed with _bu, the suffix should not be used when using this macro. Note that what goes inside must be parsable as a valid Julia expression.
Example
julia> 1.0bu"u"
+1.0 u
+julia> 1.0bu"amucc"
+1.0 amu/ccBatsrus.HDF — ModuleModule for BATSRUS HDF5 file processing.
Batsrus.HDF.BatsrusHDF5Uniform — TypeBATSRUS HDF5 file with uniform Cartesian mesh.
Batsrus.HDF.HDF5Common — TypeBATSRUS hdf5 file wrapper.
The data are stored in blocks, i.e., each field component is stored in a 4D array in the order (iblock, iz, iy, ix). This is a generic wrapper and does not assume grid type, i.e., uniform, stretched nonuniform, or AMR, etc. Classes to handle data with different grids can be derived from this class.
Batsrus.HDF.coord2index — MethodReturn lower corner index.
Batsrus.HDF.extract_var — Methodextract_var(file::BatsrusHDF5Uniform, var::String; kwargs...)Extract variable var from HDF5 file.
Keywords
xmin: minimum extracted coordinate in x.xmax: maximum extracted coordinate in x.stepx: extracted stride in x.ymin: minimum extracted coordinate in y.ymax: maximum extracted coordinate in y.stepy: extracted stride in y.zmin: minimum extracted coordinate in z.zmax: maximum extracted coordinate in z.stepz: extracted stride in z.verbose::Bool=true: display type and size information of output variable.Batsrus.HDF.global_slice_to_local_slice — Methodglobal_slice_to_local_slice(file::BatsrusHDF5Uniform, dim, gslc, ib)Convert global slice gslc to local slice lslc on a given block index ib.
Batsrus.HDF.prep_extract — Methodprep_extract(file::BatsrusHDF5Uniform, vmin=-Inf, vmax=Inf, step=1)Get info for data extraction in 1D.
Keywords
vmin, vmax: requested coordinate range (corner values).step: stride.Returns:
gslc: global slice.vmin_new, vmax_new: adjusted coordinate range (corner values).ibmin:ibmax: block range.Batsrus.HDF.prep_extract_per_block — Methodprep_extract_per_block(file::BatsrusHDF5Uniform, dim, gslc, ib)Get info for data extraction on a single block.
Arguments
gslc: global slice from prep_extract.ib::Int: block index.Returns
lslc: range to be used on the current block.ix0:ix1: index range in the global array.Batsrus.HDF.prepslice — Methodprepslice(file::BatsrusHDF5Uniform; dim::Int, vmin, vmax, step=1)Return range that covers [vmin, vmax) along dimension dim.
Returns
slc_new: trimed slice. If the object's Nx==1, then 1:1 will be returned.xl_new: adjusted lower corner coordinate matching slc_new.start.xu_new: adjusted lower corner coordinate matching slc_new.stop.Batsrus.HDF.trimslice — Methodtrimslice(start, stop, step, stop_max)Set slice's start to be nonnegative and start/stop to be within bound. Reverse slicing is not handled.
Settings
This document was generated with Documenter.jl version 1.7.0 on Tuesday 12 November 2024. Using Julia version 1.11.1.
All the workflows here is not restricted to one type of model output. After being familiar with new ideas and new models, one can easily make use of existing samples and create reader of their own. Because of the embarrassing parallelism nature of postprocessing, it is quite easy to take advantage of parallel approaches to process the data.
For the plotting, streamline tracing and particle tracing, a common problem is the grid and related interpolation process. Now I have FieldTracer.jl and TestParticle.jl designed specifically for these tasks.
If you don't have SWMF data at hand, Batsrus.jl provides some test data for you to begin with.
using LazyArtifacts
+
+datapath = artifact"testdata" # where you can find multiple test data filesThese are also used in the standard test. These will be automatically downloaded from batsrus_data if you run the package test locally.
vtkOverlappingAMR implements a somewhat strict Berger-Collela AMR scheme:
Or in other words,
previous level.
You can directly use vtkUniformGridAMR, which does not impose any restrictions. Most filters should work for this class - there just wouldn't be any specialized filters such as the dual-grid contour / clip ones for the vtkOverlappingAMR.
The vtkAMRInformation documentation consists only of
Sample 2D AMR Dataset with two levels and refinement ratio, RL=4. The root level (L0) consists of a single grid shown in black wireframe while the next level (L1) consists of two grids, depicted in green wireframe and red wireframe respectively. The two grids at L1 are projected from the root level to illustrate that the cells underneath are “hidden.”
In VTK, the collection of AMR grids is stored in a vtkHierarchicalBoxDataSet data-structure. Each grid, G(Li,k), is represented by a vtkUniformGrid data structure where the unique key pair (Li,k) denotes the corresponding level (Li) and the grid index within the level (k) with respect to the underlying hierarchical structure. An array historically known as IBLANK, stored as a cell attribute in vtkUniformGrid, denotes whether a cell is hidden or not. The blanking array is subsequently used by the mapper to hide lower resolution cells accordingly when visualizing the dataset.
To enable the execution of data queries without loading the entire dataset in memory, metadata information is employed. The metadata stores a minimal set of geometric information for each grid in the AMR hierarchy. Specifically, the AMR metadata, B(Li,k), corresponding to the grid G(Li,k), is represented using a vtkAMRBox object and it consists of the following information:
Given the metadata information stored in the AMR box of each grid, the refinement ratio at each level can be easily computed using relationship (1) from Table 1. Further, the cartesian bounds the corresponding grid covers and the number of points and cells is also available (see relationships 2-4 in Table 1). Notably, geometric queries such as determining which cell contains a given point, or if a grid intersects a user-supplied slice plane, can be answered using just the metadata.
There is a vtkAMRDualExtractionFilter, which constructs a dual-mesh (i.e., the mesh constructed by connecting the cell-centers) over the computational domain. If we can directly tell ParaView that the mesh we have is a dual-mesh, then the initial trial with multi-block data may work directly.
AMRGaussianPulseSource
See Multi-Resolution Rendering with Overlapping AMR for the implementation of C++ reader in VTK.
Settings
This document was generated with Documenter.jl version 1.7.0 on Tuesday 12 November 2024. Using Julia version 1.11.1.
file = "1d_bin.out";
+bd = load(file);
+bd = load(file, verbose=true);
+bd = load(file, npict=1);file = "3d_structured.out";
+bd = load(file, verbose=false);logfilename = "shocktube.log";
+head, data = readlogdata(logfilename)get_var_range(bd, "rho")ρ = getvar(bd, "rho")
+bd["rho"]loc = Float32[0.0, 0.0] # The type determines the output type
+d = interp1d(bd, "rho", loc)point1 = Float32[-10.0, -1.0]
+point2 = Float32[10.0, 1.0]
+w = interp1d(bd, "rho", point1, point2)Bmag = bd["Bmag"]Here is a full list of predefined derived quantities:
| Derived variable name | Meaning | Required variable |
|---|---|---|
| B2 | Magnetic field magnitude squared | Bx, By, Bz |
| E2 | Electric field magnitude squared | Ex, Ey, Ez |
| U2 | Velocity magnitude squared | Ux, Uy, Uz |
| Bmag | Magnetic field magnitude | Bx, By, Bz |
| Emag | Electric field magnitude | Ex, Ey, Ez |
| Umag | Velocity magnitude | Ux, Uy, Uz |
| Uemag | Electron Velocity magnitude | UxS0, UyS0, UzS0 |
| Uimag | Ion Velocity magnitude | UxS1, UyS1, UzS1 |
| B | Magnetic field vector | Bx, By, Bz |
| E | Electric field vector | Ex, Ey, Ez |
| U | Velocity vector | Ux, Uy, Uz |
| Anisotropy[0-1] | Pressure/Temperature anisotropy | B, P tensor |
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.
tec and tcp) and binary Tecplot file (set DOSAVETECBINARY=TRUE in BATSRUS PARAM.in):file = "x=0_mhd_1_n00000050.dat"
+convertTECtoVTU(file)gridType=1 returns rectilinear vtr file, gridType=2 returns structured vts file):file = "3d_structured.out"
+convertIDLtoVTK(file, gridType=1)filetag = "3d_var_1_n00002500"
+convertIDLtoVTK(filetag)The file suffix should not be provided for this to work correctly!
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])
+enddir = "./"
+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])
+endMore examples can be found in examples.
filename = "3d__var_1_n00006288.h5"
+file = BatsrusHDF5Uniform(filename)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)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.
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.
plot(bd, "p", linewidth=2, color="tab:red", linestyle="--", linewidth=2)scatter(bd, "p")For 2D outputs, we can select the following functions:
contourcontourfpcolormeshplot_surfaceplot_tricontourplot_tricontourfplot_trisurftripcolorwith either quiver or streamplot. By default the linear colorscale is applied. If you want to switch to logarithmic, set argument colorscale=:log.
contour(bd, "p")contourf(bd, "p")
+contourf(bd, "p", levels, plotrange=[-10,10,-Inf,Inf], plotinterval=0.1)plot_surface(bd, "p")plot_trisurf(bd, "p")tricontourf(bd, "p")streamplot(bd, "bx;bz")
+streamplot(bd, "bx;bz", density=2.0, color="k", plotinterval=1.0, plotrange=[-10,10,-Inf,Inf])quiver(bd, "ux;uy", stride=50)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.
To get the index of a certain quantity, e.g. electron number density
ρe_= findfirst(x->x=="rhoS0", bd.head.wnames)wmin, wmax = get_range(bd, var)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.
Settings
This document was generated with Documenter.jl version 1.7.0 on Tuesday 12 November 2024. Using Julia version 1.11.1.
","category":"page"},{"location":"examples/","page":"Examples","title":"Examples","text":"Line Animation","category":"page"},{"location":"examples/","page":"Examples","title":"Examples","text":" \n ","category":"page"},{"location":"examples/","page":"Examples","title":"Examples","text":"1D line animation using pyplot","category":"page"},{"location":"examples/","page":"Examples","title":"Examples","text":" \n | \n
","category":"page"},{"location":"examples/","page":"Examples","title":"Examples","text":"Subplot Line Animation","category":"page"},{"location":"examples/","page":"Examples","title":"Examples","text":" \n ","category":"page"},{"location":"examples/","page":"Examples","title":"Examples","text":"Subplot 1D line animation using pyplot","category":"page"},{"location":"examples/","page":"Examples","title":"Examples","text":" \n | \n
","category":"page"},{"location":"examples/","page":"Examples","title":"Examples","text":"Contour Animation","category":"page"},{"location":"examples/","page":"Examples","title":"Examples","text":" \n ","category":"page"},{"location":"examples/","page":"Examples","title":"Examples","text":"2D contour animation using pyplot","category":"page"},{"location":"examples/","page":"Examples","title":"Examples","text":" \n | \n