Fields

The Grid class in Schnek provides a basic method of storing data on a numerical grid. In many simulations of physical systems, the grid will represent a domain in space and the values stored in the grid represent samples or cell averages of physical quantities. This means that each value in the grid will be associated with a coordinate in the physical domain. For a large class of simulations the grid points will be evenly spaced in the physical domain. For these regular grids the coordinates of each grid point can be calculated using the grid size, the extent of the physical domain and some optional stagger of the grid. The stagger specifies for each coordinate dimension whether the value stored in the grid lives in the cell center or on the cell boundary. In the latter case, the coordinate of the grid point is shifted by half a grid cell in the positive direction.

Schnek provides the Field class to make it easier to deal with regular grids on physical domains. The Field class extends the Grid class with some extra functionality. Creating a Field is similar to creating a Grid except that the constructor takes a few more arguments.

Field<double, 3>::IndexType lo(0,0,0), hi(9,9,9);
Field<double, 3>::RangeLimit physLo(-0.5, -0.5, -0.5), physHi(0.5, 0.5, 0.5);
Field<double, 3>::RangeType physRange(physLo, physHi);
Field<double, 3>::Stagger stagger(false, true, false);

Field<double, 3> field(lo, hi, physRange, stagger, 2);

The lo and hi parameters specify the lower and upper indices of the grid’s inner domain. This is thought to be the region which corresponds to the physical domain. The physical range is specified by the second argument, physRange. This is a Rectangular Range of double type that stores the minimum and maximum of the physical coordinates represented by the field. In the example above, the physical coordinate of the grid point at index lo = (0,0,0) is (-0.5, -0.5, -0.5), whereas the physical coordinate of the grid point at index hi = (10,10,10) is (0.5, 0.5, 0.5).

The third argument specifies the staggering of the grid in each direction. The stagger array is an array of boolean. Each boolean indicates whether the grid is staggered in the corresponding direction. In the example the grid is staggered in the y-direction but not in the x or z-direction.

The last argument allows to set the number of ghost cells of the field. Ghost cells are often needed in order to implement boundary conditions. The are also needed in order to implement numerical schemes on distributed systems. In the example, the field has two ghost cells. This means that the grid is extended by two cells in each direction.

The following code illustrates the behaviour of the ghost cells.

lo = field.getLo();
hi = field.getHi();

std::cout << "Outer Lo Index (" << lo[0] << ", " << lo[1] << ", " << lo[2] << ")" << std::endl;
std::cout << "Outer Hi Index (" << hi[0] << ", " << hi[1] << ", " << hi[2] << ")" << std::endl;

This piece of code will write the following output.

Outer Lo Index (-2, -2, -2)
Outer Hi Index (11, 11, 11)

Here the lower and upper index bound of the grid has been shifted by two cells with regard to the bounds passed to the constructor. In order to obtain the bounds of the inner computational domain, excluding the ghost cells, two new methods have been added to the Field class, getInnerLo() and getInnerHi().

lo = field.getInnerLo();
hi = field.getInnerHi();

std::cout << "Inner Lo Index (" << lo[0] << ", " << lo[1] << ", " << lo[2] << ")" << std::endl;
std::cout << "Inner Hi Index (" << hi[0] << ", " << hi[1] << ", " << hi[2] << ")" << std::endl;

The code snippet above will produce the following output.

Inner Lo Index (0, 0, 0)
Inner Hi Index (9, 9, 9)

Here the lower and upper index bounds of the interior region are returned. These are the same as the lo and hi parameters passed to the constructor.

In order to convert between grid indices and physical coordinates the Field class provides the methods indexToPosition() and positionToIndex(). indexToPosition() takes two arguments, the direction and a grid index. It will return the physical coordinate of the grid cell, taking any grid stagger into account. The following code uses indexToPosition() to find the coordinate of the 5th grid cell in the x, y, and z directions.

std::cout << "x_5 = " << field.indexToPosition(0, 5) << std::endl;
std::cout << "y_5 = " << field.indexToPosition(1, 5) << std::endl;
std::cout << "z_5 = " << field.indexToPosition(2, 5) << std::endl;

The output produced by these lines is

x_5 = 0
y_5 = 0.05
z_5 = 0

Although all three indices are the same, the result of indexToPosition() is different for the y coordinate. The reason for this difference is that the grid has been set up to be staggered in the y direction but not in the x or z direction.

The inverse operation is also possible. To obtain the grid index corresponding to a physical position one can use positionToIndex().

std::cout << "Index at x=0.2: " << field.positionToIndex(0, 0.2) << std::endl;
std::cout << "Index at y=0.2: " << field.positionToIndex(1, 0.2) << std::endl;
std::cout << "Index at z=0.2: " << field.positionToIndex(2, 0.2) << std::endl;

The positionToIndex() method returns the integer part of the grid position rescaled to the grid indices. Any non-integer part is always rounded down to the next lower integer. The result of the code above is the following output.

Index at x=0.2: 7
Index at y=0.2: 6
Index at z=0.2: 7

Again the difference between the y-index and the x and z-indices is due to the different staggering on the grid in the different directions. In many cases it is important to calculate both, an integer grid position and a floating point offset, relative to the grid coordinate. This is frequently needed when interpolating values on a grid to arbitrary positions between the grid cells. For this purpose, there is an overloaded method of positionToIndex() that does not return a value nut instead takes two additional references as arguments. This is illustrated in the following example.

int index;
double offset;

field.positionToIndex(0, 0.2, index, offset);
std::cout << "Index (offset) at x=0.2: " << index << " (" << offset << ")" << std::endl;
field.positionToIndex(1, 0.2, index, offset);
std::cout << "Index (offset) at y=0.2: " << index << " (" << offset << ")" << std::endl;
field.positionToIndex(2, 0.2, index, offset);
std::cout << "Index (offset) at z=0.2: " << index << " (" << offset << ")" << std::endl;

Instead of just creating a single integer for each position, this example will return the following pairs of grid index and grid offset.

Index (offset) at x=0.2: 7 (0)
Index (offset) at y=0.2: 6 (0.5)
Index (offset) at z=0.2: 7 (0)

The code for this tutorial can be found here.

Pages

Fields