SpatialOps
Using stencils with fields

Goal

Illustrate how to use stencil operators within SpatialOps/Nebo

Key Concepts

  1. There are a number of predefined operators/stencils in SpatialOps. The simplest way to obtain these is via the BasicOpTypes struct. Much of this was covered in Inferring operator types
  2. Stencils can be stored and retrieved in the OperatorDatabase which provides simple storage and retrieval of operators by type.
  3. SpatialOps defines many stencils for use on structured, finite-volume meshes. The build_stencils function can be used to easily populate an OperatorDatabase with the predefined operators in SpatialOps.

Example Code

examples/stencils.cpp

#include <spatialops/structured/Grid.h>
using namespace SpatialOps;
// If we are compiling with GPU CUDA support, create fields on the device.
// Otherwise, create them on the host.
#ifdef ENABLE_CUDA
# define LOCATION GPU_INDEX
#else
# define LOCATION CPU_INDEX
#endif
#define PI 3.141592653589793
int main()
{
//----------------------------------------------------------------------------
// Define the domain size and number of points
const DoubleVec length(PI,PI,PI); // a cube of length pi on each side
const IntVec fieldDim( 10, 1, 1 ); // a 1-D problem
//----------------------------------------------------------------------------
// Create fields
const bool bcx=true, bcy=true, bcz=true;
const IntVec bc(bcx,bcy,bcz);
const GhostData nghost(1);
const BoundaryCellInfo sVolBCInfo = BoundaryCellInfo::build<SVolField >( bc, bc );
const BoundaryCellInfo ssxBCInfo = BoundaryCellInfo::build<SSurfXField>( bc, bc );
const MemoryWindow sVolWindow( get_window_with_ghost( fieldDim, nghost, sVolBCInfo) );
const MemoryWindow ssxWindow( get_window_with_ghost( fieldDim, nghost, ssxBCInfo ) );
SVolField x( sVolWindow, sVolBCInfo, nghost, NULL, InternalStorage, LOCATION );
SVolField f( sVolWindow, sVolBCInfo, nghost, NULL, InternalStorage, LOCATION );
SSurfXField dfdx( ssxWindow, ssxBCInfo, nghost, NULL, InternalStorage, LOCATION );
SSurfXField fface( ssxWindow, ssxBCInfo, nghost, NULL, InternalStorage, LOCATION );
//----------------------------------------------------------------------------
// Build a grid. This is a convenient way to set coordinate values that will
// be used below.
const Grid grid( fieldDim, length );
grid.set_coord<XDIR>(x);
//----------------------------------------------------------------------------
// Build the operators (stencils). Here we will use predefined stencils that
// can be built easily.
OperatorDatabase opDB; // holds stencils that can be retrieved easily
build_stencils( grid, opDB ); // builds stencils and places them in opDB
typedef BasicOpTypes<SVolField>::GradX GradX; // x-derivative operator type
typedef BasicOpTypes<SVolField>::InterpC2FX InterpX; // x-interpolant operator type
const GradX& gradx = *opDB.retrieve_operator<GradX >(); // retrieve the GradX operator
const InterpX& interpx = *opDB.retrieve_operator<InterpX>(); // retrieve the InterpX operator
//----------------------------------------------------------------------------
// perform calculations
f <<= sin( x ); // Initialize the function value
dfdx <<= gradx( f ); // gradient of f at the x-faces
fface <<= interpx( f ); // interpolated value of f at the x-faces
return 0;
}

Try this

Modify stencils.cpp to do the following:

  1. Use two interpolants to calculate f at faces and then interpolate back to cell centers:
    f2 <<= interpf2c( interpc2f( f ) );
    To do this:
    • Obtain the interpolant types. Note that interpc2f is already there, and you can use BasicOpTypes
      BasicOpTypes<SVolField>::InterpF2CX
      to get the operator type for face to cell interpolation.
    • Build the f2 field (use dfdx and fface as examples).
  2. Instead of building dfdx and fface as shown in the example, use the SpatialFieldStore as shown in Basics of field creation
    SpatFldPtr<SSurfXField> dfdxPtr = SpatialFieldStore::get<SSurfXField>( x );
    Note that you will have a pointer now, so will need to de-reference it as appropriate when using it in nebo statements.
See also
Inferring operator types