Analyzing Wombat Data
We are using Julia to post-process Wombat simulations and produce publication ready plots. Julia offers native parallelism and near C speed in loops with the convenience of Python and IDL. It also interfaces very well with Wombat's Fortran 2008 (e.g. same array order and indexing) and exposes key HPC concepts like vectorization (dot syntax) and parallelization. However, it is structured in many ways like Python (modules, encapsulation), which makes for a very productive to work environment and a great playground to learn programming. Differently to Python, loops are fast and can be parallelized without much work.
Julia has reached version 1.0 end of 2018, which broke some of our code. Until everything is fixed, you can use 0.7 which shows deprecation warnings.
Start Julia like this to optimize for your machine, in this case Intel haswell architecture with 4 cores and -O3 flag:
julia -C haswell -O 3 -p 4
Available architectures are shown with
julia -C help
Our Julia library including serial snapshot reading & plotting examples are part of the wombat source tree and can be found in the "lib/julia/" directory.
You first want to add the Wombat julia library to your search path in your ~/.julia/config/startup.jl
@everywhere push!(LOAD_PATH, "/Users/jdonnert/Dev/src/git/wombat/lib/julia")
Reading Wombat Data
In Julia the next command reads the whole grid density into an array 'rho' and the simulation properties into a structure 'sim':
using Plots
using Wombat
gr() # using GR backend, fastest.
rho, sim = ReadSnap("KH2d0-010-000", "RHO")
vx, sim = ReadSnap(sim, "VX")
vy, sim = ReadSnap(sim, "VY")
vz, sim = ReadSnap(sim, "VZ")
pos = MakeZonePos(sim)
pos[:,1] -= sim.boxsize[1]
pos[:,2] -= sim.boxsize[2]
pos[:,3] -= sim.boxsize[3]
v = sqrt.(vx.^2 + vy.^2 + vz.^2)
r = sqrt.(pos[:,1].^2 + pos[:,2].^2 + pos[:,3].^2)
plot(r,v) # radial velocity profile centered on box
ReadSnap also takes the "sim" structure as an input to avoid processing the manifest every time you're reading data from the same simulation. "sim" carries all info necessary to read the patch data and can become large (10MB). It carries all snapshot information, except the actual data.
We plan to implement parallel analysis using the same patch based approach to the data as Wombat does internally. N Julia tasks read M patches in parallel, resolve boundaries if necessary, compute on every patch and reduce the result (e.g. projection & map making). This should make for an efficient post processing pipeline.
Dot Syntax
Notice in this example the dot . syntax in the example above. The dot tells Julia to operate element wise on an array. This seems dumb, if you're used to Python or IDL, until you realize that Julia is as fast as C because of this. If you're operating on pre-defined arrays you can replace all dots with a macro:
N = 128
x = rand(N,N,N)
y = rand(N,N,N)
r = similar(x)
@. r = exp(x^2 + y^2)^5.9875
This expression will SIMD vectorize and approach C with gcc compiler (not Fortran though).
You can introspect functions and check SIMD vectorization :
function test(x::Array{Float64,3}, y::Array{Float64,3},r::Array{Float64,3})
@. r = exp(x^2 + y^2)^5.9875
return r
end
code_native(test, (Array{Float64,3},Array{Float64,3},Array{Float64,3}))
The last command will show the assembler output of the function test(). There should be vector instructions starting with "v" operating on xmm or ymm registers. Try doing this with Python !
Plots for Analysis and Animations
We recommend the Plots package with the GR backend to make quick plots for analysis and animations/movies. GR is the fastest backend and is being continuously improved. It now plots millions of points in less than a second inlined into iTerm.
using Plots
import ColorBrewer
# init Plots with GR
Plots.gr(color_palette=ColorBrewer.palette("Set1", 9),legend=false,
grid=true, markerstrokewidth=1, markersize=1, linewidth=1, size=(512,512/sqrt(2)))
Then plot something:
Plots.plot(rand(100))
The following code will produce an mp4 movie of an image sequence generated from Wombat data with Kovesi's colormaps:
function KH()
rootdir = "/Users/jdonnert/Volumes/iCloud/Projects/2017/WENO-Wombat/"
datadir = rootdir*"/data/2D_KH/"
clibrary(:colorcet)
col = :bgyw
anim = @animate for i = 1:50
fname = datadir*"compr-"*itoa(i,4)*"-000" # itoa() available in wombat's julia lib
data, sim = ReadSnap(fname, "RHO")
Plots.heatmap(data, size=size(img), show=true, aspect_ratio=1, size=psize,
colorbar=:best, clims=(0.4,1.2), c=col)
end
end
mp4(anim, "test.mp4")
Plots is a very powerful generic package. Refer to the Plots documentation here for more examples and color palette definition.
What if it doesn't work ?
If you get errors like this:
qt.qpa.screen: QXcbConnection: Could not connect to display
Could not connect to any X display.
connect: Connection refused
GKS: can't connect to GKS socket application
Did you start 'gksqt'?
GKS: Open failed in routine OPEN_WS
GKS: GKS not in proper state. GKS must be either in the state WSOP or WSAC in routine ACTIVATE_WS
then you are on a remote machine and GR cannot use the X-forwarding. You may try:
ssh -YC USERNAME@MACHINE
julia
julia> ENV["GRDIR"]=""
julia> Pkg.build("GR")
julia> quit()
julia
julia> using Plots; gr(); histogram(randn(10000))
Another workaround is to let xterm handle the forwarding. To do this open an xterm on the remote machine:
ssh -X USERNAME@MACHINE
xterm
and then start julia inside this xterm. Alternatively, you can tell GR's GKS window system that it is "headless" with before starting julia:
export GKSwstype=100
julia
It will then not attempt to output to a window, and you will have to save your plots into a png file:
savefig( FILENAME )
A Word on Colors
The visual representation of data is a non-trivial undertaking. Colors in particular can bias the impression the viewer has on the data. Additionally, there are issues regarding black & white conversion (lines still distinguishable ?), reproduction on a printer and color blind viewers.
We recommend, colormaps either from the ColorBrewer package (lines in plots) or the Perceptual Colormaps (images). These maps are well behave and ensure good readability of results. They are also very pretty. See colorbrewer2.org and this paper.
The Color Brewer packages are available in julialang via
Pkg.add("ColorBrewer")
import ColorBrewer
and can be made available in the Plots package and the gr backend via:
Plots.gr(color_palette=ColorBrewer.palette("Set1", 9),legend=false, grid=true)
and in PGF via our PGFTools (in wombat/lib/julia/misc) package:
using PGFTools
InitPGF()
which has to be called only once per julia session.
The perceptual Color maps are available in Plots via
julia> using Plots
julia> gr()
julia> clibrary(:colorcet)
julia> plot(rand(128,128), aspect_ratio=1, c=:colorwheel, clims=[0.1,0.9], colorbar=:best)
The list of available colors is here. :colorcet is the perceptually uniform palette recommended above.
Plots for Publications
Highest quality plots for papers can be made with PGFPlots. It generates native LaTex code from the plot commands. This stunning in papers, which are all LaTeX, but the syntax takes some time to learn and is not very intuitive.
A PGF example:
import ColorBrewer
import PGFPlots
import PGFTools # Julius' PGF tools that allow heat map and contour plots.
using LaTeXStrings
PGFTools.InitPGF() # init PGF settings
x = range(0, stop=2*pi, length=256)
y = sin.(x) .* sin.(4*x)
y2 = cos.(0.5*x) .* cos.(8*x)
l1 = PGFPlots.Plots.Linear(x,y, style="Set11,solid", mark="x", legendentry="A sin function") # a line
l2 = PGFPlots.Plots.Linear(x,y2, style="Set12,dashed", mark="none", legendentry="Fun !") # a line
ax = PGFPlots.Axis([l1,l2], ylabel=L"$f(x) = \sin(x)$", xlabel="x",xmin=0, # the axis
xmax=2*pi, ymin=-1.1, ymax=1, ymode="normal", xmode="normal" )
PGFPlots.save("/Users/jdonnert/Desktop/test.pdf", ax) # save it
Consult the PGFPlots documentation or directly the PGF manual. For more examples, you can get the plotting scripts for one of our papers. The MHD-Weno paper contains lots of examples. The Julia source code for the plots is available from the website in the Publications section under the "Source" button.