Bugfix (ProbeVolume): Output IS (with natural ordering) to re-arrange sub-data (post-processing)

[mesnardo]

Bug affects simulations running on multiple MPI processes, when using a volume probe to monitor the solution in a specific window of the computational domain, and when this window spans over multiple domain processes.

When monitoring the solution in a window (volume) that is smaller than the computational domain, we can use the class petibm::misc::ProbeVolume.
The field data in the window can be accessed using the PETSc routine VecGetSubVector.
The routine takes as input parameter an index set that contains index (PETSc ordering) of the data points to grab in the PETSc vector.
When data are spread over multiple processors, the routine VecGetSubVector returns a copy.
PETSc does this by creating a vector scatter context and using it to scatter data from one vector to another (VecScatterBegin and VecScatterEnd).
However, it simply concatenates in the parallel ordering of the vector.
Therefore, the sub vector is not arranged in natural order when output to HDF5 file.

The solution implemented here is to create a second index set, which contains the natural index of the points in the window being monitored, and to save it into the HDF5 file.
At post-processing stage, users can load the index set from HDF5 file and re-arrange the window data in natural ordering.
Here is a small Python example to re-arrange data (in 2D):

import h5py
import numpy


f = h5py.File(str(filepath), 'r')
# Load gridline positions.
x, y = f['mesh']['x'][:], f['mesh']['y'][:]
# Load index set (natural index of positions).
index_set = f['mesh']['IS'][:].flatten()
# Load field values of the volume.
time_str = '0.100000'
values = f['u'][time_str][:]
# Re-arrange values based on index set.
mask = numpy.argsort(index_set)
values = values[mask].reshape((y.size, x.size))
f.close()

Example used to showcase the problematic behavior: 2D flow generated by a cylinder translating at speed 1 in the -x direction. (Simulations ran on 2 processors, Intel(R) Core(TM) i7-3770 CPU @ 3.40GHz, and using 1 K20 GPU to solve the pressure Poisson system) Figure below shows the filled contour of the x-velocity in the domain at time value 3.51. We also add the filled contour of the component from the data monitored with the volume probe (window [-4, 1] x [-1, 1]). Left panel shows the solution with PetIBM-0.5; right panel shows the solution using the branch bugfix-probe where we also output the index set to re-arrange the data in the window; sub-data match data from the computational domain (expected behavior). (Red line represents the value x-velocity, relative to the translational velocity of the cylinder, interpolated at y=0 along an horizontal gridline. Expected behavior: relative velocity near zero in the cylinder, negative just behind the cylinder, approaching 1 far away from the cylinder.)

Change(s) in API: None.

Change(s) in source:

• Add methods writeIS, writeIS_HDF5, writeIS_ASCII, in the class ProbeVolume, to write the Index Set (natural ordering) into file.
• Add new attribute isNatural (class ProbeVolume) that stores the index set for the sub-volume that contains natural index of points to track in the domain.
• Rename attribute is (class ProbeVolume) into isPetsc; it contains the index (PETSc ordering) of the point to track.
• Rename few methods in the class ProbeVolume.

Change(s) in output files:

• Index Set (containing the natural index of the points tracked in the volume) is saved into the output file for the probe. The indices can be used to re-order the sub-vector of data before visualization.

PR ready for review.
Once reviewed, we can update the changelog and bump the version. (-> Done.)