SpatialOps
The 3D Laplacian in action

Using the laplacian operator, $ \frac{\partial^2}{\partial x} + \frac{\partial^2}{\partial y} + \frac{\partial^2}{\partial z} $, as a case-study, this example shows how to create a function template to compute the 3D laplacian on compatible field types.

Goals

  • Illustrate how to use stencil operators
  • Illustrate type inference from fields to operators

Key Concepts

  1. There are a number of predefined operators/stencils in SpatialOps. The simplest way to obtain these is via the BasicOpTypes struct.
  2. Stencils can be stored and retrieved in the OperatorDatabase which provides simple storage and retrieval of operators by type.
  3. Using type inference, highly generic code can be written that is also very robust.
See also
Some natively supported field types
Inferring operator types
Using stencils with fields

Example Code

examples/3D_Laplacian.h

// Here FieldT is a generic field that is assumed to live on a volume.
//
// Valid input field types:
// ------------------------
// SVolField - cell-centered (scalar) volume field
// XVolField - x-staggered volume field (x-velocity/momentum)
// YVolField - y-staggered volume field (y-velocity/momentum)
// ZVolField - z-staggered volume field (z-velocity/momentum)
//
// Inputs:
// -------
// opDB - the OperatorDatabase that holds all operators
// src - the source field (to apply the Laplacian to)
//
// Output:
// -------
// dest - the destination field (result of Laplacian applied to src)
//
template< typename FieldT >
void
calculate_laplacian_3D( const SpatialOps::OperatorDatabase& opDB,
const FieldT& src,
FieldT& dest )
{
using namespace SpatialOps;
//---------------------------------------------------------------------------
// Infer operator types for the 3D laplacian
// If this function is used with a field type that is not a volume field
// then compiler errors result due to invalid type inference for operators.
typedef typename BasicOpTypes< FieldT >::GradX GradX;
typedef typename BasicOpTypes< FieldT >::GradY GradY;
typedef typename BasicOpTypes< FieldT >::GradZ GradZ;
typedef typename BasicOpTypes< FieldT >::DivX DivX;
typedef typename BasicOpTypes< FieldT >::DivY DivY;
typedef typename BasicOpTypes< FieldT >::DivZ DivZ;
//---------------------------------------------------------------------------
// obtain the appropriate operators from the operator database
// The OperatorDatabase is simply a container for operators that
// provides access by type.
GradX& gradX = *opDB.retrieve_operator<GradX>();
GradY& gradY = *opDB.retrieve_operator<GradY>();
GradZ& gradZ = *opDB.retrieve_operator<GradZ>();
DivX& divX = *opDB.retrieve_operator<DivX >();
DivY& divY = *opDB.retrieve_operator<DivY >();
DivZ& divZ = *opDB.retrieve_operator<DivZ >();
//---------------------------------------------------------------------------
// Finally, perform the calculation. This is the easy part!
// Note that this will run on multicore and GPU!
dest <<= divX( gradX( src ) )
+ divY( gradY( src ) )
+ divZ( gradZ( src ) );
}