Julia is a high-level, high-performance dynamic programming language for technical computing, with syntax that is familiar to users of other technical computing environments. You can obtain Julia from [http://julialang.org], or if you are using Ubuntu/Debian, you can install it with apt-get install julia.

This is the documentation for the Julia module for the Pencil Code. This module contains convenience functions for post-processing data files from the Pencil Code. To use this module, add this to your ~/juliarc.jl so that Julia can find the module.

push!(LOAD_PATH,  ENV["PENCIL_HOME"] * "/julia")

ubuntu$ julia        # Start the Julia interpreter.
...
julia> using Pencil  # Load the module.

NOTE: At the present time, you also need to add the push! line at the top of stand-alone programs.

Plotting and graphics

Julia has several plotting packages. The one I like is PyPlot which uses Python's Matplotlib library. To install PyPlot run the following:

    ubuntu$ sudo apt-get install python-matplotlib

    ubuntu$ julia
    ...
    julia> Pkg.add("PyPlot")

    julia> using PyPlot

    julia> x = linspace(0,2*pi,1000);
    julia> y = sin(3*x + 4*cos(2*x));

    julia> plot(x, y, color="red", linewidth=2.0, linestyle="--")

    julia> title("A sinusoidally modulated sinusoid")

    #
    # Save figure as a PNG:
    #
    julia> savefig("myfigure.png")

    #
    # LaTeX in labels and titles.
    #
    julia> title(L"Plot of $\Gamma_3(x)$")  # L for LaTeX.

    #
    # Colour mesh: plot a 2D grid.
    #
    julia> y = [1:128] * ones(128)';  # Col vector x Row vector

    julia> r2 = (y - 64).^2 + (y' - 64).^2;

    julia> pcolormesh(r2)

    julia> axis([0,128,0,128]) # [xmin, xmax, ymin, ymax]

    julia> savefig("colmesh.png")

    #
    # 3D plotting
    #
    julia> surf(r2)

    julia> mesh(r2)

The PencilPlot module provides a color map called "density" that may be useful when plotting a density variable like rhopmxz. PencilPlot loads PyPlot. Example usage:

    ubuntu$ cd /my/simulation/dir
    ubuntu$ julia
    ...
    julia> using Pencil

    julia> using PencilPlot

    julia> rhopmxz = read_yaver(it=10)["rhopmxz"];

    julia> axis([0, 128, 0, 128])

    #
    # PencilPlot defines the "density" color map.
    #
    julia> pcolormesh(rhopmxz, cmap=ColorMap("density") )

    julia> savefig("rhopmxz.png")

Read time series

Provides:

A function to read the time series from time_series.dat as a Dict(), with one key for each variable in time_series.dat.

Tutorial:

Function	Summary
`read_ts()`	Returns a `Dict()` with keys like `rhopmax`, etc.

  julia> using Pencil

  julia> ts = read_ts();

  julia> typeof(ts)
  Dict{Any,Any} (constructor with 2 methods)

  #
  # List of available keys / variables found in timeseries.dat
  #
  julia> ts["keys"]
  35-element Array{String,1}:
   "it"       
   "t"        
   "dt"       
   "walltime" 
   "nblockmin"
   "nblockmax"
   "nmigmax"  
   "nparmin"  
   "nparmax"  
   "nparbmax" 
   "rhom"     
   "rhomax"   
   "rhomin"   
   "dtdragp"  
   "urms"     
   ...          
   "uxuym"    
   "TTm"      
   "TTmax"    
   "vpxm"     
   "vpx2m"    
   "vpz2m"    
   "xpm"      
   "xp2m"     
   "zpm"      
   "zp2m"     
   "npmax"    
   "rhopm"    
   "rhopmax"  
   "dedragp"  
   "decollp"  

  julia> typeof(ts["rhopmax"])
  Array{Float64,1}

  julia> length(ts["rhopmax"])
  9528

  #
  # Plotting with Python's Matplotlib.
  #
  julia> using PyPlot

  julia> plot( ts["t"] , ts["rhopmax"] )

Optional parameters:

Read particles

Provides:

Option	Description
`datadir="xxx"`	Path to the data directory (default: "data").

Two functions to read particle data from pvar.dat or the PVAR* files. The difference is only in the output: read_pvar is more consistent with most functions in this module in that in returns a Dict with keys like x and vx, while read_particles returns an array of Particle objects which is more cache-friendly for numerical computation.

Tutorial:

Function	Summary
`read_pvar()`	Returns a `Dict()` with keys like `x` and `vx`.
`read_particles()`	Returns an array of `Particle` objects.

  julia> using Pencil

  julia> pars = read_particles();
  INFO: mpvar = 6
  INFO: Read 16384 particles

  #
  # "pars" is an array of 16384 particles.
  #
  julia> typeof(pars)
  Array{Particle,1}

  julia> size(pars)
  (16384,)

  #
  # Each particle has a position (x,y,z) and velocity (u,v,w)
  #
  julia> pars[1].x
  -0.058057133821926316

  julia> pars[1].u
  0.0003632957506802519

  #############

  julia> pvar = read_pvar();
  INFO: mpvar = 6
  INFO: Read 16384 particles

  julia> typeof(pvar)
  Dict{Any,Any} (constructor with 2 methods)

  #
  # The "keys" key has a list of dictionary keys / variables read (from index.pro)
  #
  julia> pvar["keys"]
  6-element Array{Any,1}:
  "y" 
  "vx"
  "vy"
  "x" 
  "z" 
  "vz"

  julia> typeof(pvar["x"])
  Array{Float64,1}

  julia> size(pvar["x"])
  (16384,)

Optional parameters:

TODO:

Read averages

Provides:

Option	Description
`datadir="xxx"`	Path to the data directory (default: "data").
`snapshot=n`	If n > 0, return data from file `PVAR$n`.

Functions to read 2D and 1D averages. These functions return all variable averages as a Dict() unless the var option is specified. They also return all time steps, unless the option it (for "iteration") is specified.

Tutorial:

Function	Summary
`read_xyaver()`	Read `xyaverages.dat` and return the variables as a `Dict()`
`read_xzaver()`	Read `xzaverages.dat` and return the variables as a `Dict()`
`read_yzaver()`	Read `yzaverages.dat` and return the variables as a `Dict()`
`read_xaver()`	Read `xaverages.dat` and return the variables as a `Dict()`
`read_yaver()`	Read `yaverages.dat` and return the variables as a `Dict()`
`read_zaver()`	Read `zaverages.dat` and return the variables as a `Dict()`

  julia> using Pencil

  julia> xyaver = read_xyaver();

  #
  # Self-inspection. List of keys in xyaver.
  #
  julia> collect( keys(xyaver) )
   "uxmz"  
   "ux2mz" 
   "oumz"  
   "uz2mz" 
   "uymz"  
   "keys"  
   "uy2mz" 
   "uxuymz"
   "rhomz" 
   "uzmz"  
   "t"     
   "rhopmz"

  #
  # The key "keys" gives you a list of the available keys / variables read.
  #
  julia> xyaver["keys"]
   "t"     
   "uxmz"  
   "uymz"  
   "uzmz"  
   "ux2mz" 
   "uy2mz" 
   "uz2mz" 
   "uxuymz"
   "oumz"  
   "rhomz" 
   "rhopmz"

  #
  # size() returns (rows,columns). (3416,) is a column vector.
  #
  julia> size( xyaver["t"] ) # (rows,columns)
  (3416,)

  #
  # nit    == number of iterations
  # ngridz == number of cells in z
  # 
  julia> nit, ngridz = size( xyaver["rhomz"] )
  (3416,128)

Optional parameters:

TODO:

Read dimensions

Provides:

Option	Description
`var="rhopmax"`	Read only one variable instead of the entire Dict().
`datadir="xxx"`	Path to the data directory (default: "data").
`varfile="xxx"`	Name of the data file (default: "xyaverages.dat").
`infile="xxx"`	Variables names file (default: "xyaver.in").
`it=12`	If n > 0, return only that iteration (starts at 1).

Functions to parse dim.dat and pdim.dat and return the result as a Dict(). The first file contains general dimensions like nx, mx, nghostx (where mx == nx + 2*nghostx), mvar, nprocx, etc. The second file contains particle dimensions like npar.

Tutorial:

Function	Summary
`read_dim()`	Reads `dim.dat` and returns a `Dict()`.
`read_pdim()`	Reads `pdim.dat` and returns a `Dict()`.

  julia> using Pencil

  #
  # Read data/dim.dat
  #
  julia> dim = read_dim();

  julia> typeof(dim)
  Dict{Any,Any} (constructor with 2 methods)

  julia> dim["keys"]
  33-element Array{Any,1}:
   "mvar"     
   "ny"       
   "m2"       
   "mx"       
   "my"       
   "nghosty"  
   "nygrid"   
   "nghostx"  
   "precision"
   "nprocx"   
   "ipz"      
   "nxgrid"   
   "ipy"      
   ...
   "nprocz"   
   "mxgrid"   
   "mz"       
   "nz"       
   "maux"     
   "ipx"      
   "mzgrid"   

  #
  # Compare data/dim.dat  vs data/proc10/dim.dat
  #

  julia> read_dim()        # Read data/dim.dat

  julia> dim["nprocz"]
  16

  julia> dim["ipz"]
  -1

  julia> read_dim(proc=10) # Read data/proc10/dim.dat

  julia> dim["nprocz"]
  -1

  julia> dim["ipz"]
  10

  #
  # read_pdim() has a similar API.
  #
  julia> pdim = read_pdim();

  julia> typeof(pdim)
  Dict{Any,Any} (constructor with 2 methods)

  julia> pdim["keys"]
  3-element Array{Any,1}:
   "mpvar"     
   "npar"      
   "npar_stalk"

  julia> pdim["npar"]
  16384

Option	Description
`datadir="xxx"`	Path to the data directory (default: "data").
`proc=n`	Read the `dim` or `pdim` file from the `data/proc$n` directory.

Julia module for the Pencil Code

Introduction

Plotting and graphics

Read time series

Provides:

Tutorial:

Optional parameters:

Read particles

Provides:

Tutorial:

Optional parameters:

TODO:

Read averages

Provides:

Tutorial:

Optional parameters:

TODO:

Read dimensions

Provides:

Tutorial:

Optional parameters: