Container classes (vectors, matrices, ND-arrays, etc.).
NdFlattener is a container to store a deep copy of multidimensional data in a single object that provides a generic []...[] interface. The data and their structure are stored in few contiguous arrays, allowing for an efficient MPI gathering operation.
Which problem does it solve?
In SU2, there is a lot of multidimensional data: data depending on K>1 indices, like a wall property depending on the rank, zone index, and marker index.The data is often stored non-contiguously as a "tree":
- an array of pointers, each one of which points to
- (repeat the previous line another (K-2) times)
- an array of data.
The individual arrays' sizes are usually unrelated to each other. Therefore in addition to this tree of pointers-to-...-to-data, usually there is a tree of pointers-to-...-to-size for each lower level. Instead of arrays, there might be getter functions for the data and the sizes.In order to use multidimensional data behind such a non-generic interface as an argument in a function call, we usually make the callee aware of the interface and expose all the trees: Either as separate pointers to arrays (or functions), which is tedious, or by handing an object in which they naturally live, which exposes much more information.If the callee is an MPI communication operation, this is not possible. Here, a frequent pattern is to copy the data into a single contiguous array, which is then communicated alongside with the structure (lengths of the arrays).The
NdFlattener container is a systematic way to store multidimensional data in a single contiguous array, and the structure in few additional contiguous arrays. It exposes them through a generic []...[] interface and allows for an efficient MPI_Allgatherv-like operation.To demonstrate the usage of
NdFlattener, let us store a property of type su2double that depends on the rank, zone index, and marker index.
by a "recursive function" like this:
auto f_local =
make_pair( nZone, [=](unsigned long iZone){
return make_pair( (Geometry of iZone)->GetnMarker() , [=](unsigned long iMarker){
return YOUR_PROPERTY(iZone, iMarker);
});
});
NdFlattener<2> nd_local(f_local);
It might be safer to explicitly capture what you need in the lambda expressions.You can also construct an NdFlattener without any arguments, and use NdFlattener::initialize(f_local).f_local is
- a pair of an (exclusive) upper bound for the first index and a function that maps the first index to
- a pair of an (exclusive) upper bound for the second index and a function that maps the second index to
- (iterate this if there are more layers)
- the value corresponding to these indices. The template parameter "2" is the number of indices. The data type (default: su2double) and the index type (default: unsigned long) are further optional template parameters.
by "collective communication" like this:
NdFlattener<3> nd_global(Nd_MPI_Environment(), &nd_local);
nd_global's first index is the rank, then the indices of nd_local follow.The struct Nd_MPI_Environment contains the relevant information about the parallel environment. The default values are suitable to communicate su2double data. In order to communicate data of another MPI datatype, use e.g. Nd_MPI_Environment(MPI_UNSIGNED_LONG).You can also construct an NdFlattener without any arguments, and use NdFlattener::initialize(mpi_env, &nd_local).
Refreshing
If only the data but not the indices (i.e. the structure of sublists' lengths) have changed, you can refresh them with the call
nd_local.refresh(f_local); // , or
nd_global(Nd_MPI_Environment(), &nd_local); // respectively.
An NdFlattener constructed or initialized from a "recursive function" or "collective communication" must be refreshed in the same way.We also provide a function NdFlattener::initialize_or_refresh, that initializes if not initialized, and refreshes otherwise.
Look-up
You can access the
NdFlattener via a "[]...[]" interface that resembles multidimensional arrays or pointer-to-pointer-... arrays. If you provide all of the K indices, this returns a reference to the data element. If you provide less than K indices, you obtain an IndexAccumulator object to which you can later pass more indices. With respect to the above example:
std::cout << nd_global[rank][zone][marker]; // e.g. 4.0
auto nd_g_rank = nd_global[rank];
nd_g_rank[zone][marker] *= 2.0;
std::cout << nd_g_rank[zone][marker]; // e.g. 8.0
You may use the IndexAccumulator repeatedly as long as the accessed NdFlattener exists. As NdFlatteners are always deep copies, the above block would not change nd_local or the non-contiguous data from which it was copied.Each index is checked to lie in the correct range dictated by the NdFlattener and the previous indices. You can obtain the minimal overflowing index with size:
auto nd_g_rank = nd_global[rank];
std::cout << nd_g_rank[zone].size();
When all indices except the last one are fixed, the corresponding data is stored contiguously and a pointer to this 1D array can be retrieved by
nd_global[rank][zone].data()
Constness
The reference returned by the last
operator[], and the pointer returned by
data(), are const if
- the accessed NdFlattener object was const, or
- one of the intermediate IndexAccumulator objects was qualified as const.
This reflects the semantic meaning of const:
NdFlattener<2> nd2(...);
template<class T> set(T& t) { t[0] = 1.0; }
set( nd2[0] ); // OK
template<class T> set_const(T const& t) { t[0] = 1.0; }
set_const( nd2[0] ); // does not compile
Unit tests
The interface described here is tested in UnitTests/Common/toolboxes/ndflattener_tests.cpp, which may also serve as an addition example of how to use the
NdFlattener.
If your array has K indices, instantiate this class with template parameter K, which is derived recursively from this class with template parameter (K-1). In each layer, there is an array: Of type Data for K=1, and of type Index for K>1.The data array of K=1 contains the values of the array A in lexicographic ordering: [0]...[0][0][0], [0]...[0][0][1], ..., [0]...[0][0][something], [0]...[0][1][0], [0]...[0][1][1], ..., [0]...[0][1][something], ..., [0]...[1][0][0], [0]...[1][0][1], ..., [0]...[1][0][something], ... ..., [1]...[0][0][0], [1]...[0][0][1], ..., [1]...[0][0][something], ... ... ... Let us call each row in this representation a "last-index list". Note that they might have different lengths, so "something" might stand for different values here. Last-index lists can also be empty.The indices array of K=2 contains the indices of the (K=1)-data array at which a new last-index list starts. If a last-index list is empty, the are repetitive entries in the indices array.The indices array of K=3 contains the indices of the (K=2)-indices array at which the last-but-two index increases by one or drops to zero. If the last-but-two index is increased by more than one, there are repetitive entries.Etc. etc, up to the indices array at layer K.To form such a structure, we typically iterate twice through the pointer-to-pointer-... array: The first time we get to know how much space to reserve in each layer, then we allocate it and during the second iteration we copy the data.
1.9.8