Homec4science
Diffusion tamaas (features/netcdf-mpi) Browse

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

Last Change
Raw File

grid.org
View Options

#+TITLE: Grid
* Description
The =Grid= class is a 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[], L[], 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.
* 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.
* 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
Looping over points of a grid and setting the proxy object is easy but a bit verbose. Here is an example:
#+BEGIN_SRC cpp
UInt N = grid.getNbPoints();
VectorProxy<Real> proxy(nullptr, nb_components);
for (UInt i = 0 ; i < N ; ++i) {
proxy.setPointer(grid.getInternalData() + i * proxy.dataSize());
// proxy is now a vector wrapping on a point of the grid
}
#+END_SRC
Looping with OpenMP requires the proxy objects to be private to each thread:
#+BEGIN_SRC cpp
UInt N = grid.getNbPoints();
VectorProxy<Real> main_proxy(nullptr, 1);
VectorProxy<Real> aux_proxy(nullptr, 1);
// firstprivate for anything declared outside of for-loop
#pragma omp parallel for firstprivate(main_proxy, aux_proxy)
for (UInt i = 0 ; i < N ; ++i) {
main_proxy.setPointer(grid.getInternalData() + i * main_proxy.dataSize());
aux_proxy.setPointer(other_grid.getInternalData() + i * main_proxy.dataSize());
// do calcultations
}
#+END_SRC
* Miscellaneous
Other functions are defined for =Grid=:
- =dataSize()= returns the size of the flattened array
- =sizes()= returns all sizes
- =lengths()= returns all physical lengths
- =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