Step 6 Computing stencils

Step 6 - Computing stencils

In this tutorial we will learn how to gather a stencil consisting of data from the current element and its face neighbors with special emphasis on exchanging data over the ghost elements beforehand.

You will find the code to this example in the tutorials/general/step6* files and it creates the executables tutorials/general/t8_step6_stencil.

In the last tutorials we learned how to create a forest, adapt it, pre-allocate element data arrays and store custom data fields in VTU files. In this tutorial we will start by performing all these operations in one go as shown in step 5. Then, when we have our forest and built a data array, we gather data for the local elements of our process. Next, we exchange the data values of the ghost elements and compute a stencil from face neighbors in each element. Finally, the output of several custom data fields is written to .vtu files.

In terms of interacting with t8code, this step mainly builds on the functionality already detailed in the steps 3 to 5. Thus, we focus on the new concept of accessing neighboring element data and doing computations with it by reference of a small finite difference computation.

Computing Finite Differences from Face Neighbor Stencils

For many scientific and industrial applications the access to neighboring data of an element is a crucial feature every mesh library has to provide. The collection of the element's data and its neighbors is called stencil. Stencils can be seen as a compact data structure which allows for convenient and performant finite difference computations yielding numerical approximations of gradients, curvatures, curls, etc. of a data field. An exemplary visualization of a 3x3 cross stencil is shown in the figure below. The uniform grid is divided into three parallelization zones with the already introduced ghost layers at their interfaces.

The data field in our specific case is the physical quantitiy "height" pinned at the midpoint of each mesh element. As an example, we compute the schlieren and a rough messure for the curvature which are the defined as the norms of the central difference formulas shown in the figure above. The actual implementation in the demo is as follows:

const double xslope_m = 0.5/(dx[0] + dx[1])*(stencil[1][1] - stencil[0][1]);
const double xslope_p = 0.5/(dx[1] + dx[2])*(stencil[2][1] - stencil[1][1]);

const double yslope_m = 0.5/(dy[0] + dy[1])*(stencil[1][1] - stencil[1][0]);
const double yslope_p = 0.5/(dy[1] + dy[2])*(stencil[1][2] - stencil[1][1]);

const double xslope = 0.5*(xslope_m + xslope_p);
const double yslope = 0.5*(yslope_m + yslope_p);

const double xcurve = (xslope_p - xslope_m)/0.25/(dx[0] + 2.0*dx[1] + dx[2]);
const double ycurve = (yslope_p - yslope_m)/0.25/(dy[0] + 2.0*dy[1] + dy[2]);

element_data[current_index].schlieren = sqrt(xslope*xslope + yslope*yslope);
element_data[current_index].curvature = sqrt(xcurve*xcurve + ycurve*ycurve);

The above code also accounts for adapted meshes with stencil computations crossing different refinement levels. The gathering of the stencil data is described in the next section.

Accessing Element Face Neighbors in t8code & Gathering the data into a Stencil

t8code provides an API function to access neighboring elements at a given face. The excerpt from t8_forest.h reads:

/** Compute the leaf face neighbors of a forest.
 * \param [in]    forest  The forest. Must have a valid ghost layer.
 * \param [in]    ltreeid A local tree id.
 * \param [in]    leaf    A leaf in tree \a ltreeid of \a forest.
 * \param [out]   neighbor_leafs Unallocated on input. On output the neighbor
 *                        leafs are stored here.
 * \param [in]    face    The index of the face across which the face neighbors
 *                        are searched.
 * \param [out]   dual_face On output the face id's of the neighboring elements' faces.
 * \param [out]   num_neighbors On output the number of neighbor leafs.
 * \param [out]   pelement_indices Unallocated on input. On output the element indices
 *                        of the neighbor leafs are stored here.
 *                        0, 1, ... num_local_el - 1 for local leafs and
 *                        num_local_el , ... , num_local_el + num_ghosts - 1 for ghosts.
 * \param [out]   pneigh_scheme On output the eclass scheme of the neighbor elements.
 * \param [in]    forest_is_balanced True if we know that \a forest is balanced, false
 *                        otherwise.
 * \note If there are no face neighbors, then *neighbor_leafs = NULL, num_neighbors = 0,
 * and *pelement_indices = NULL on output.
 * \note Currently \a forest must be balanced.
 * \note \a forest must be committed before calling this function.
 */
void                t8_forest_leaf_face_neighbors (t8_forest_t forest,
                                                   t8_locidx_t ltreeid,
                                                   const t8_element_t *leaf,
                                                   t8_element_t
                                                   **pneighbor_leafs[],
                                                   int face,
                                                   int *dual_faces[],
                                                   int *num_neighbors,
                                                   t8_locidx_t
                                                   **pelement_indices,
                                                   t8_eclass_scheme_c
                                                   **pneigh_scheme,
                                                   int forest_is_balanced);

A common usage pattern of above routine can be found in the function body of t8_step6_compute_stencil in the demo source file tutorials/general/t8_step6_stencil.cxx. There at each element face, which is given by face argument above, we can decide where to put the data in the stencil. This is done via a switch case construct:

/* `neighids` is a return argument from above `t8_forest_leaf_face_neighbors` call. */
double height = element_data[neighids[0]].height;
switch (iface) {
  case 0: // NORTH
    stencil[0][1] = height;
    break;
  case 1: // SOUT
    stencil[2][1] = height;
    break
 case 2: // WEST
    stencil[1][0] = height;
    break;
 case 3: // EAST
    stencil[1][2] = height;
    break;
}

Once the stencil is filled the finite difference computation from the previous section can be applied. Above code is just an excerpt with details left out. Please refer to the demo code in order to get the full picture.

Running the Demo

If t8code was compiled with MPI support the demo can be exectuted with for example three processes.

$ mpirun -n 3 ${INSTALL_DIR_OF_T8CODE}/bin/t8_step6_stencil

This results in three VTU files and a meta PVTU file.

$ ls
t8_step6_stencil.pvtu  t8_step6_stencil_0000.vtu  t8_step6_stencil_0001.vtu  t8_step6_stencil_0002.vtu

Open t8_step6_stencil.pvtu in Paraview and visualize the height (left figure) and schlieren (right figure) data fields. Schlieren plots are the norm of the gradient of a data field (in our case the height). For smoother pictures the inclined reader is encouraged to increase the initial refinement level, recompile and rerun the demo.