4_data_visualisation

We provide here an overview of the python routines that display time series and process volumes and slices exported by the solver. After trying these few routines on the exemple below, the reader is advised to test his ability to reproduce the three benchmarks (1, 2, 3). Not only does this hands-on exercise will help him become much more familiar to the routines; it will also validate the installation of Coral against three known cases. NB: benchmark3 is perhaps a little computationally intensive for mere workstations.

Time series

---

After your first run, the Timeseries subdirectory should contain a few files, depending on how you set up the coral.timeseries file:

ls Timeseries/
tFull_XYavg_z00023.dat  tFull_XYavg_z00035.dat  time.dat  tw_volAvg.dat  uu_volAvg.dat  vv_volAvg.dat  ww_volAvg.dat

For plotting a family of volume-averaged time series, run the plot_volAvg_timeseries.py python routine with a list of variable as suffix. For instance, for time series of kinetic energy along each axis, type:

python plot_volAvg_timeseries.py uu vv ww

For horizontally-averaged quantities, run the plot_xyAvg_timeseries.py python routine, with a list of pairs of names and indices for z. For instance, for plotting the full temperature field at index 23 and 35, enter:

python plot_xyAvg_timeseries.py tFull 23 tFull 35

Upon examination of these files, you will soon realize that plotting the time series contained in the file my_var.dat as a function of time can be as simple as:

import numpy as np
import matplotlib.pyplot as plt

time  = np.fromfile('time.dat',   dtype=np.float_)
myVar = np.fromfile('my_var.dat', dtype=np.float_)
plt.plot(time, myVar)
plt.show()

---

Volumes

---

Full vs dealiased data

By default Coral exports volume on the full collocation grid. For clarity, we recall that if (Nx,Ny,Nz) modes are requested by the user in the coral.parameters.in file, alias-free quadratic terms are computed on a grid the size (3/2 Nx, 3/2 Ny, 3/2 Nz). Volumes are exported in the ./Volumes subdirectory of a given run. Full-grid volumes can be recognized by their _full.dat suffix.

Considering that the grid-size is larger than the actual number of modes (due to the dealiasing procedure), keeping these "full" files is wasteful. We could interpolate these volumes on a looser grid the size (Nx, Ny, Nz) without loosing information, saving a factor 27/8 in disk space. From your ${PROJECTS} directory, all the "full" volumes of the my_run folder will be interpolated on a smaller grid with the command:

python dealiase_volumes.py ./my_run/

After completion of the call, the folder ${PROJECTS}/my_run/Volumes will contain both full files and dealiased (smaller) files, recognizable with their _twoThirds.dat prefix. Full files may be erased at this point.

By hand

Any file in the folder ${PROJECTS}/my_run/Volumes can be loaded into a numpy array using:

# for full files
a = np.fromfile('./filename_full.dat', dtype=np.float_).reshape(3*NX//2,3*NY//2,3*NZ//2)
# for dealiased files
a = np.fromfile('./filename_twoThirds.dat', dtype=np.float_).reshape(NX,NY,NZ)

Notice that data is stored with Z being the fastest index (contiguous in memory), and X being the slowest.

Finally the collocation grid is equispaced in the X and Y direction, as it should be for Fourier modes. The Gauss-Chebyshev grid is used along Z.

Although it is convenient to be able to read the data yourself, for convenience, you may prefer to use the python class provided in process_data.py. This class, described below, takes care of all the technicalities (determining the resolution, defining grids, etc.) for you, so that your focus can be on exploring the data, as opposed to python book-keeping.

---

Using the existing Python routines

Loading some data

A python class for loading (and processing) volumes exists in process_data.py. Let us get familiar with this class working on an example in an interactive python session (I like ipython; Jupyter is also very popular). Once you're familiar with the functionalities, use them in your own python scripts.

In the following, In[] and Out[] are the iPython interface, and of course should be ignored. Start by loading the class:

In [1]: from process_data import plane_layer_volume

Load some volumes into the class instance--let us name it a. Suppose we are interested in three variables--the first linear variable, the second quadratic variable, the fourth linear variable, as defined in coral.equations--at timestep 200:

In [2]: a = plane_layer_volume( list_var_str=['linear_var01', 'quadra_var02', 'linear_var04'], time_int=200)

The class contains a look-up table, which displays a confirmation of which volumes are currently loaded :

In [3]: a.lut  
Out[3]: ['linear_var01', 'quadra_var02', 'linear_var04']

The volumes that corresponds to variable a.lut[i] (with i an integer) is stored in a.dat[i] in the form of a rank-3 numpy array, e.g.:

In [4]: a.dat[2].shape    
Out[4]: (64, 64, 48)

Differentiating Volumes

Volumes can be differentiated by calling .xdiff, .ydiff or .zdiff methods, and providing their index. For differentiating w.r.t. z quadra_var02, located at index 1 according to the look-up table (recall that python is 0-based indexing...)enter:

In [5]: a.zdiff(1)
# possibly displays an efficiency warning. should be ignored
In [6]: # let us check that a new volume has been created and added to the list:
In [7]: a.lut                                                                                                        
Out[7]: ['linear_var01', 'quadra_var02', 'linear_var04', '[d/dz]quadra_var02']
In [8]: # the last volume has the same shame as the others
In [9]: a.dat[-1].shape                                                                                              
Out[9]: (64, 64, 48)
In [10]: # we may differentiate the same (or another) volume again, say w.r.t. x this time                            
In [11]: a.xdiff(3)              
In [12]: a.lut                  
Out[12]:
['linear_var01',
 'quadra_var02',
 'linear_var04',
 '[d/dz]quadra_var02',
 '[d/dx][d/dz]quadra_var02']

Horizontal averaging, Vertical profiles

Volumes can be horizontally averaged using .horizontally_average(index). This appends a profile (i.e. a rank-1 array that depends on z only) to the list a.profile. This list is referenced by a separate look-up table: a.lut_profiles. Very similarly to volumes, the profile listed as a.lut_profiles[i] (with i an index) is stored as a rank-1 numpy array in a.profiles[i]:

In [24]: a.horizontally_average(1)
In [25]: a.horizontally_average(3)
In [26]: a.lut_profiles          
Out[26]: ['<quadra_var02>_{x,y}', '<[d/dz]quadra_var02>_{x,y}']
In [28]: a.profiles[0].shape            
Out[28]: (48,)

Display profiles

Profiles have values on the Gauss-Chebyshev grid, with is the a.zgrid attribute. For plotting a profile, we merely need to do:

In [30]: import matplotlib.pyplot as plt                                                                              
In [31]: plt.plot(a.profiles[0], a.zgrid) #this results in a "vertical" plot
Out[31]: [<matplotlib.lines.Line2D at 0x7f8f3e5eb0b8>]
In [32]: plt.show()   

Display slices of a volume

Display a slice of a given volume with .plot_slice. The keyword argument varInt refers to the volume position in the look-up table (starting at 0). A second keyword argument must be entered: x=i or y=i or z=i, depending on the kind of slice, where i is the index on the collocation grid. (Reminder: the Gauss Chebyshev grid is upside down: 0 is the top, the endpoint -1 is the bottom.):

In [33]: fig=a.plot_slice(varInt=4, z=10)
In [34]: fig
Out[34]: <Figure size 640x480 with 1 Axes>
In [35]: fig.show()       
In [36]: fig=a.plot_slice(varInt=2, x=0)
In [37]: fig.show()

---

Slices

---

Above, we have learned how to slice a volume. Of course, it is possible to directly export a slice on the disk at runtime. These slices are stored in the folder ./Slices of a given run. The filename format is: varName_kindOfSlice_timeInteger_full.dat where:

• varName is either linear_var?? or quadra_var?? with ?? a (possibly zero-padded) integer that refers to the variable, as defined in coral.equations.
• kindOfSlice informs us on which index (x,y, or z) is kept constant. The value of this index follows, padded by zeros on 5 characters.
• timeInteger refers to the time step of the given output, padded on 8 characters.

In sum, slices files may have names such as linear_var01_z00032_time00000200_full.dat or quadra_var04_z00014_time00000300_full.dat.

Colorplot with colorplot_slice.py

Among the python scripts that are automatically copied by the prepare_directory.sh script during the creation of a new folder for a run, you will find colorplot_slice.py. This script is useful for plotting any data that is in the ./Slices folder. From a terminal, call the python script, and follow the guidance offered by the prompt:

[ben@occigen]> python3 colorplot_slice.py 
Linear (l) or Quadratic (q) variable? [l/q] >> l
Var. number? [integer]                      >> 4
Slice position? [eg x00001, y00512, z00064] >> y00001
Time? [integer]                             >> 7000

Alternatively, you can provide a matplotlib colormap name as command line argument, e.g.:

[ben@occigen]> python3 colorplot_slice.py RdYlBu

Advanced usage

For opening slices by hand, remember the index ordering : Z is fastest than Y which is fastest than X. Also remember that, by contrast with volumes, slices are not dealiased. They are always plotted on the full grid, with the 3/2 factor along both directions. Hence, to load an XY-, a YZ-, or an XZ-slice into a numpy array, use:

XYSlice = np.fromfile('./Slices/varName_z00014_time00000200_full.dat', dtype=np.float_).reshape(3*NX//2, 3*NY//2)
XZSlice = np.fromfile('./Slices/varName_y00010_time00001000_full.dat', dtype=np.float_).reshape(3*NX//2, 3*NZ//2)
YZSlice = np.fromfile('./Slices/varName_x00001_time00001000_full.dat', dtype=np.float_).reshape(3*NY//2, 3*NZ//2)