Flow Over a Pebble

OCCA block

This block is used to specify whether NekRS is running on CPU or GPU (and which vendor). For our case, we will run this on CPU so we set the backend as CPU.

GENERAL block

The [GENERAL] block describes the time stepping, simulation end control, and polynomial order. A NekRS output file is written every 500 time steps. The numSteps field indicate how many time steps to take. We will use a polynomial order of . We will begin with our first time step as seconds, but allow NekRS to adaptively increase this to get a target Courant-Friedrichs-Lewy (CFL) number of 2.0.

MESH block

The [MESH] block in this example is used to list the ordering of sidesets as they pertain to the ordering of the boundary conditions in the boundaryTypeMap parameters in other blocks. In other words, the first boundary condition given in these lists corresponds to mesh sideset 2, and so on.

FLUID VELOCITY, FLUID PRESSURE, SCALAR TEMPERATURE blocks

Next, the [FLUID VELOCITY] and [FLUID PRESSURE] blocks describe the solution of the pressure Poisson equation and velocity Helmholtz equations. The [SCALAR TEMPERATURE] block describes the solution of the energy equation. NekRS can also solve for an arbitrary number of coupled passive scalars, which would be represented using [SCALAR FOO] blocks, where FOO would be replaced with a string describing the name of the scalar (e.g. [SCALAR K] for the turbulent kinetic energy when using a Reynolds Averaged Navier-Stokes (RANS) model).

In these blocks, residualTol is used to indicate the solver tolerance. Here, you also specify the type of boundary conditions to apply to each sideset (you only specify boundary conditions in the [FLUID VELOCITY] and [SCALAR FOO] blocks, because the pressure and velocity solves are really indicating together a solution to the momentum conservation equation). The boundaryTypeMap is used to specify the mapping of boundary IDs to types of boundary conditions. NekRS uses character strings to refer to the type of boundary condition; full information can be found on the NekRS website.

Some of the common boundary condition strings are shown in Table 2. When you specify these boundary conditions, the boundary conditions requiring user information (e.g. udfDirichlet, udfNeumann, and udfRobin) will trigger the requirement for you to provide functions by the same name in the .oudf file.

When you populate boundaryTypeMap in the input file, you simply list the character string for your desired boundary condition in the same order as the sidesets. If your sidesets are already ordered sequentially, nothing else needs to be done. Otherwise, you need to provide the boundaryIDMap field in the [MESH] block as described earlier.

In each of the equation blocks, different parameter names are used to indicate the properties which are multiplying the time and diffusive terms.

For convenience, NekRS allows you to pull parameters from the .par file elsewhere throughout your case files. In order to streamline the model setup, you can add any number of user-defined blocks by setting the userSections parameter and then providing those blocks. In this example, we illustrate this for a [CASEDATA] block where we define local variables which we arbitrarily name name Re, D, and inlet_T which we'll use in our boundary and initial conditions.

.udf File

The .udf file contains additional C++ code which is typically used for setting initial conditions and postprocessing. This file at a minimum must contain 2 functions (they can be empty):

• UDF_Setup is called on initialization time, and is typically where initial conditions are applied

• UDF_ExecuteStep is called once on each time step, and is typically where postprocessing is applied

static dfloat inlet_v;
static dfloat inlet_T;
static dfloat Re;
static dfloat D;
static dfloat mu;
static dfloat rho;

#ifdef __okl__

#include "pebble.oudf"

#endif

void UDF_Setup0(MPI_Comm comm, setupAide &options)
{
  // get parameters from .par file
  platform->par->extract("casedata", "D", D);
  platform->par->extract("casedata", "Re", Re);
  platform->par->extract("casedata", "inlet_T", inlet_T);
  platform->options.getArgs("FLUID VISCOSITY", mu);
  platform->options.getArgs("FLUID DENSITY", rho);
}

void UDF_LoadKernels(deviceKernelProperties & kernelInfo)
{
  // send parameter to device
  inlet_v = Re * mu / D / rho;
  kernelInfo.define("p_inlet_v") = inlet_v;
  kernelInfo.define("p_inlet_T") = inlet_T;
}

void UDF_Setup()
{
  auto mesh = nrs->meshV;

  // arbitrary arrays on the host where we will fill some initial data
  std::vector<dfloat> U(nrs->fluid->fieldOffsetSum, 0.0);
  std::vector<dfloat> temp(mesh->Nlocal, 0.0);

  for (int n = 0; n < mesh->Nlocal; ++n)
  {
    U[n + 0 * nrs->fieldOffset] = 0.0; // x-velocity
    U[n + 1 * nrs->fieldOffset] = 0.0; // y-velocity
    U[n + 2 * nrs->fieldOffset] = inlet_v; // z-velocity

    temp[n] = inlet_T; // temperature
  }

  // copy host data into the device data
  nrs->fluid->o_U.copyFrom(U.data(), U.size());
  nrs->scalar->o_solution("temperature").copyFrom(temp.data(), temp.size());
}

void UDF_ExecuteStep(double time, int tstep)
{
}

Here, we use the UDF_Setup function to set initial conditions (to the same values we applied as our inlet boundary conditions). Because NekRS is a GPU code, most data structures only exist on the device. Therefore, in order to set initial conditions, we declare some empty vectors of the appropriate lengths for velocity and temperature, fill them with our desired initial conditions, and then copy them into the device.

NekRS stores the various solution components on the device as

• nrs->fluid->o_U is velocity (all three components are packed one after another, with each "slice" of length nrs->fieldOffset)

• nrs->fluid->o_P is pressure

• nrs->scalar->o_solution("foo") is the passive scalar named FOO in the .par file

By default, initial conditions are taken to be zero unless otherwise specified; because we did not specify an initial condition for pressure, this means our initial pressure will be zero in this problem.

Our problem only has one passive scalar (temperature). You will have additional passive scalars if doing RANS modeling (e.g. and would be passive scalars) or if you have additional conservation equations for mass transport, etc. We apply the initial condition by looping over all the degrees of freedom on the present MPI rank.

Our .udf file has two additional (optional) functions which we use to set up other aspects of the simulation. The UDF_Setup0 function is where we extract parameters from the .par file. We fetch the values we provided in the [CASEDATA] block, as well as the viscosity and density we gave in the [FLUID VELOCITY] block.

Finally, the UDF_LoadKernels function is where we can load GPU kernels and send parameters to the device. The kernelInfo lines in particular are sending the inlet_v and inlet_T values from host to device, so that we can later access them from the .oudf file. It is customary when sending a value to device to prepend p_ to the name, to help with clarity.

.oudf File

When you define the boundary condition types in the .par file, you also need to set the values of those boundary conditions (the non-trivial ones) n the .oudf file. The names of these functions are pre-defined by NekRS, and are paired up to the character strings you set earlier. These function names are shown in Table 2 in the .oudf column. Note that a --- indicates that no user-defined kernels are necessary to define that boundary condition (i.e. NekRS already has all the info it needs to apply a zero-Dirichlet boundary condition - it just sets the solution to zero).

For each of these functions, the bcData struct contains all information about the current boundary that is "calling" the function:

• bc->id is the boundary ID

• bc->uxFluid is the -velocity

• bc->uyFluid is the -velocity

• bc->uzFluid is the -velocity

• bc->sScalar is the scalar solution at the present quadrature point; if there are multiple scalars, this should be placed inside an if (isField("scalar foo")) to ensure you only set the value of the FOO scalar.

• bc->fluxScalar is the flux of the scalar at the present quadrature point; if there are multiple scalars, this should be placed inside an if (isField("scalar foo")) to ensure you only the the flux of the FOO scalar.

• bc->idxVol is the index corresponding to the current quadrature point

void udfDirichlet(bcData *bc)
{
  if (isField("fluid velocity"))
  {
    bc->uxFluid = 0.0; // x-velocity
    bc->uyFluid = 0.0; // y-velocity
    bc->uzFluid = p_inlet_v; // z-velocity
  }

  if (isField("scalar temperature"))
    bc->sScalar = p_inlet_T; // temperature
}

void udfNeumann(bcData *bc)
{
  // power divided by sphere surface area
  bc->fluxScalar = 315892;
}

.usr File

NekRS is under rapid development, and many features from Nek5000 are currently being ported over to NekRS. To access those Fortran-based features, you can optionally include a .usr file, which contains funtions for setting initial conditions, modifying the mesh, etc. We do not need this file for this case, so we will defer to a later tutorial to demonstrate it.