Homec4science
Diffusion tamaas (features/netcdf-mpi) Browse

tamaas/doc/grid.org
d757eb5d4125features/netcdf-mpi

Last Change
Raw File

grid.org
View Options

#+TITLE: Grid
* Description
The =Grid= class is the base container of /Tamaas/. It has two template parameters:
- =T= is the type stored
- =dim= is the dimension of the grid
Common combinations of those template parameters are already instanced in [[../src/grid.cpp][grid.cpp]].
The grid points can be viewed as points in space (physical or Fourier). On each of these points, the =Grid= class allows for storage of an arbitrary number of components. These can be viewed as the scalar, vectorial or higher-order tensorial fields that are stored per-point of the grid.
All data is contiguous in memory, so reinterpretation of the data (such as change in number of components) is possible. It also allows for simple loops for element-wise operations.
* Creation
There are two constructors for =Grid<T, dim>=:
- =Grid()= initializes an empty grid
- =Grid({n}, nb_components)= initializes a grid with specified number of elements, physical dimensions and number of components per point of grid.
* Accessing elements
The =()= operator is implemented to give access to a given component of a given point, or to give access to a given cell in the flattened array, e.g. for a =Grid<Real, 2>= (a surface):
- =grid(0, 3, 2)= will return the 3rd component of the (0, 3) point
- =grid(5)= will return the 6th cell of the flattened array
Calling the =()= with a number of arguments different from (1 or dim+1) will result in an error.
Normally individual points on a grid should not be accessed directly, but through iterators.
* Vectorial operations
The following vectorial operations are defined for both scalar and =Grid= arguments:
- ~+=~
- ~-=~
- ~*=~
- ~/=~
- ~=~
The ~=~ operator calls the =copy()= function, which allows for copies between grids with different types (or the =move()= function).
* Looping with iterators
Because of OpenMP constraints, iterator loops look like this:
#+BEGIN_SRC cpp
auto end_it = primal->end();
#pragma omp parallel for
for (auto it = primal->begin() ; it < end_it ; ++it) {
// do something with *it
}
#+END_SRC
Due to the nature of OpenMP for loops, any additional iterator needs to be offset inside the loop, e.g. :
#+BEGIN_SRC cpp
#pragma omp parallel for firstprivate(t_it, r_it, q_it) reduction(+:num, denum)
for (auto it = primal->begin() ; it < end_it ; ++it) {
// Putting everyone in the correct spot
if (*it > 0) {
// Offsetting other iterators
UInt offset = it - begin_it;
t_it = t_b; t_it += offset;
r_it = r_b; r_it += offset;
q_it = q_b; q_it += offset;
num += *q_it * *t_it;
denum += *r_it * *t_it;
}
}
#+END_SRC
* Looping with multi-component field
The =Grid= class mother of two proxy classes:
- =VectorProxy= which acts as a vector
- =MatrixProxy= which acts as a matrix
These classes /wrap/ on the memory of a grid: they do not own it. They have only one component per point, and define the =()= operator with dim number of arguments (1 for the vector and 2 for the matrix). They allow simple linear algebra:
- matrix-vector multiplication: =VectorProxy::mul()=
- vector dot product: =VectorProxy::dot()=
- more to be implemented
Here's an example of OpenMP loop with proxy assignment
#+BEGIN_SRC cpp
#pragma omp parallel for
for (auto it = freq.begin(dim) ; it < freq.end(dim) ; ++it) {
VectorProxy<T> wavevector(&(*it), dim);
// do stuff with wavevector
}
#+END_SRC
* Miscellaneous
Other functions are defined for =Grid=:
- =dataSize()= returns the size of the flattened array
- =sizes()= returns all sizes
- =getNbComponents()= returns the number of components
- =setNbComponents()= sets the number of components
- =getInternalData()= returns a pointer to the first data
- the =<<= operator is defined so you can print the content of a grid
c4science · Help