Several tools are available for post-processing and visualization of data produced by the XSHELLS code:
xsplot.py: a python module for plotting 1D or 2D data.
xsplot: command-line interface to xsplot.py.
xspp: command-line tool for extracting slices, profiles, spectra and more from the field files.
pyxshells.py: a python module that can load and handle the binary XSHELLS field files.
The python module xsplot is provided to load and display data from XSHELLS. It can be used interactively or within scripts. Such Python scripts using matplotlib and xsplot are located in the matplotlib dierectory, and can be called from command line.
xsplot can also be used directly from command line and it will guess the type of file and display it accordingly.
The python modules should be installed by calling make install-py, or explicitly python setup.py install --user from the python subdirectory. In addition, it is convenient to copy the small script python/xsplot in your path, so that you can invoke xsplot from any directory.
Time-dependent quantities, such as the energies and other custom diagnostics which are stored in the energy.job file (see section 2.2.1) are also easy to plot with the xsplot tool.
Example
By default, xsplot displays kinetic and magnetic energy as a function of time:
xsplot energy.jobA list of all available quantities is also shown. To plot other quantities (see section 2.2.1), for instance kinetic energy and viscous dissipation, use: xsplot energy.job -c Eu,D_nu |
The load_diags function in the xsplot module can be used to retrieve conveniently any diagnostic. It returns a dictionary of time series for easy plotting.
Example
To plot the magnetic to kinetic energy ratio as a function of time:
> from pylab import * > import xsplot > d = xsplot.load_diags(’energy.job’) > t = d[’t’] > plot(t, d[’Eb’]/d[’Eu’]) # plot magnetic to kinetic energy ratio |
Fields are stored in binary files (see 1.6), using a custom format. They can be handled after the simulation by the xspp command line program. Alternatively, the pyxshells python module can also read, write and handle these files (see section 3.2.5).
Compile the program by typing make xspp. Invoking it without arguments (by running ./xspp) will print a help screen including the commands and their syntax.
Example
The following will display information about the file fieldU.job (resolution, precision, time of the snapshot, …):
./xspp fieldU.jobTo compute the energy and maximum value of the curl of the field: ./xspp fieldU.job curl nrj maxTo extract the field values along a line spanning the -axis from to , and also display total energy of field: ./xspp fieldU.job line -1,0,0 0.8,0,0 nrjAdd two fields and save the result to a new file (the first file will set the resolution for the result): ./xspp fieldT_0004.job + fieldT0.job save fieldT_total_0004.jobExtract only a given range of spherical harmonic coefficients (2 to 31) and computes the corresponding energy: ./xspp fieldB.job llim 2:31 nrj |
Note that xspp is not parallelized using MPI, so that for very big cases you might run out of memory (although it can operate out-of-core – without actually loading the whole file in memory – in some cases). As a workaround you can always reduce the spherical harmonic truncation while reading your big files with the llim option (see example avobe).
One of the most common usage for xspp is to extract two-dimensional slices of the 3D data stored in spectral representation in the field files. Four types of 2D slices are available:
Meridian cuts (a plane containing the -axis), with the merid command;
Equatorial cuts (the plane ), with the equat command;
Surface data (on a sphere of given radius ), with the surf command;
Disc cuts (an arbitrary plane), with the disc command;
When these commands are given to xspp, NumPy files corresponding to the required cuts are written to the current directory. These NumPy files (*.npy) can then be loaded and displayed using Python with NumPy and matplotlib (see next section).
Example
A meridian cut at :
./xspp fieldU.job meridAn equatorial cut, and a meridian cut at degrees, of the vorticity (curl of U) ./xspp fieldU.job curl equat merid 45Extract the field at the spherical surface closest to , using only the symmetric components. ./xspp fieldU.job sym 0 surf 0.9Make a cut at , using 200 azimuthal points, with field truncated at harmonic degree 60: ./xspp fieldU.job llim 0:60 disc 200 0,0,0.7 |
The python module xsplot is provided to load and display cuts produced by xspp. It can be used interactively or within scripts. Such Python scripts using matplotlib and xsplot are located in the matplotlib dierectory, and can be called from command line. xsplot can also be used directly from command line and will guess the type of cut of your file and display it accordingly.
The python module should be installed by calling make install-py, or explicitly python setup.py install --user from the python subdirectory. In addition, it is convenient to copy the small script python/xsplot in your path, so that you can invoke xsplot from any directory.
Example
Produce a meridian and an equatorial cut with xspp:
./xspp fieldU.job merid equatFrom command prompt, quickly load and plot all components of the field in this meridional slice, as well as in the equatorial plane. xsplot o_merid_0.npy o_equat.npyxsplot has several convenient options. For instance you can plot the first two components only (using -c option), and specify a range for the colormap (using the -z option): xsplot o_equat_0.npy -c 0,1 -z "(-1e-3,2e3)"Alternatively, from an Ipython interpreter (or notebook, or script), load and plot the -component of the field in the meridional and equatorial slices: > import xsplot > a = xsplot.load_slice_npy(’o_merid_0.npy’) > xsplot.plot_slice(a, 2) # plot third (phi) component > d = xsplot.load_slice_npy(’o_equat.npy’) > xsplot.plot_slice(d, 1) # plot second (phi) component |
Several useful scripts for basic and advanced plotting can be found in the matplotlib directory.
From the paraview website: ParaView is an open-source, multi-platform data analysis and visualization application. ParaView users can quickly build visualizations to analyze their data using qualitative and quantitative techniques.
Full spatial fields can be saved to XDMF format, which can be loaded by paraview. Note that the HDF5 library is required for this to work, and must be found by the configure script. If so, Simply run:
make xspp ./xspp fieldB_0004.job hdf5 B_cartesian.h5
The file B_cartesian.h5.xdmf describes the cartesian components of the vector field B on a spherical grid that can be read directly by paraview (if prompted for a loader, select ’XDMF’).
For more complex post-processing, xspp may not be enough. The python module pyxshells allows you to quickly write your own scripts to work directly with the spectral fields stored in the field files output by XSHELLS, cast them to spatial domain, and so on.
type make install-py. You can then import the pyxshells module from any python script. Note that the pyxshells module requires the shtns python module, which can be installed from the shtns directory using:
./configure --enable-python --disable-openmp make python setup.py install --user
Example
Here is an example of what pyxshells can do:
> import pyxshells > f = pyxshells.load_field(’fieldU.job’) > r = f.grid.r # the radial grid > KE = f.energy() # computes the energy > f.sht.set_grid() # prepare the spherical grid > u3D = f.spat_full()ΨΨ# the 3D field in spherical coordinates > ur,ut,up = f.spat_shell(len(f.grid.r)-1) # surface field > ur,ut,up = f.spat_merid() # meridional slice > ur,ut,up = f.spat_equat() # equatorial slice |