Flow Over a Pebble

.re2 Mesh

NekRS uses a mesh in a custom .re2 format.

If you have a mesh in Exodus or Gmsh format, you can convert that mesh into .re2 format using tools that ship with Nek5000. To get these features, you will need to clone Nek5000 somewhere on your system and then build the exo2nek tool (for converting Exodus meshes) or the gmsh2nek tool (for converting Gmsh meshes).

git clone https://github.com/Nek5000/Nek5000.git
cd Nek5000/tools
./maketools exo2nek

Running the above will place a binary named exo2nek at Nek5000/bin/exo2nek. We recommend adding this to your path.

Now that you have exo2nek, we are ready to convert your mesh (pebble.exo) into .re2 format. This mesh is Hex8. Run the following:

exo2nek

Follow the on-screen prompts. For this case, we have:

• 1 fluid exo file, which is named pebble.exo (you only need to provide pebble as the name)

• 0 solid exo files (this is nonzero when NekRS is solving for both fluid and solid domains; in this example, we apply a surface heat flux, but are not actually solving for temperature in the solid.)

• 0 periodic surface pairs (this is nonzero when applying periodic boundary conditions)

• output file name is pebble.re2 (you only need to provide pebble as the output name)

This will create a mesh named pebble.re2. NekRS is a spectral element code, which means that the solution in each element is represented as an -th order Lagrange polynomial (in each direction). An illustration for a 5th-order NekRS solution is shown in Figure 2 for a 2-D element. Each red dot is a node (GLL quadrature). When you create a mesh for NekRS, you do not see these GLL points in your starting Exodus/Gmsh mesh. Instead, they will be created when you launch NekRS (in other words, you do not need to create unique meshes if you want to run NekRS at different polynomial orders).

Figure 2: Nodal positions for a 5-th order polynomial solution.

.par File

The .par file is used to set up the high-level settings for the case. This file consists of blocks (in square brackets) and parameters.

[OCCA]
  backend = CPU

[GENERAL]
  polynomialOrder = 3
  numSteps = 4000
  dt = targetCFL=2.0 + initial=1e-4
  writeInterval = 500

[PRESSURE]
  residualTol = 1e-06

[VELOCITY]
  boundaryTypeMap = v, o, W, sym
  residualTol = 1e-07
  density = 1000.0
  viscosity = 1e-3

[TEMPERATURE]
  boundaryTypeMap = t, o, f, I
  residualTol = 1e-07
  rhoCp = 4186000
  conductivity = 0.6

[CASEDATA]
  Re = 50
  D = 0.03
  inlet_T = 0

$NEKRS_HOME/bin/nekrs --help par

.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 3 functions (they can be empty):

• UDF_LoadKernels will load any user-defined physics kernels

• 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

#include "udf.hpp"

static dfloat inlet_v;
static dfloat inlet_T;

void UDF_LoadKernels(occa::properties & kernelInfo)
{
  // get parameters from .par file
  dfloat D, Re, mu, rho;
  platform->par->extract("casedata", "D", D);
  platform->par->extract("casedata", "Re", Re);
  platform->par->extract("casedata", "inlet_T", inlet_T);
  platform->options.getArgs("VISCOSITY", mu);
  platform->options.getArgs("DENSITY", rho);

  // send parameter to device
  inlet_v = Re * mu / D / rho;
  kernelInfo["defines/p_inlet_v"] = inlet_v;
  kernelInfo["defines/p_inlet_T"] = inlet_T;
}

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

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

    nrs->P[n] = 0.0; // pressure

    nrs->cds->S[n + 0 * nrs->fieldOffset] = inlet_T; // temperature
  }
}

void UDF_ExecuteStep(nrs_t *nrs, 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). NekRS stores the various solution components on the host as

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

• nrs->P is pressure

• nrs->cds->S are the passive scalars (packed sequentially one after another, each of length nrs->cds->fieldOffset[i])

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 elements in the mesh, and in each element for all the GLL points.

Next, we fetch the parameters we desired from the .par file in UDF_LoadKernels. We fetch the values we provided in the [CASEDATA] block, as well as the viscosity and density we gave in the [VELOCITY] block. 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 and Table 3 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 no-slip boundary condition - it just sets velocity 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->u is the -velocity

• bc->v is the -velocity

• bc->w is the -velocity

• bc->s is the scalar (temperature) solution at the present quadrature point

• bc->flux is the flux (of temperature) at the present quadrature point

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

• bc->scalarId is the passive scalar number (in this example temperature is the 0th scalar and we have no additional scalars)

void velocityDirichletConditions(bcData *bc)
{
  bc->u = 0.0; // x-velocity
  bc->v = 0.0; // y-velocity
  bc->w = p_inlet_v; // z-velocity
}

void scalarDirichletConditions(bcData *bc)
{
  bc->s = p_inlet_T; // temperature
}

void scalarNeumannConditions(bcData *bc)
{
  // power divided by sphere surface area
  bc->flux = 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. Here, we will show you how to use this file for changing the sideset IDs in our mesh. Our input mesh, pebble.exo, did not order the sidesets beginning from 1; instead, the sidesets were numbered as 2, 3, 4, and 5. So, we would like to modify those sideset numbers so that they are 1, 2, 3, 4.

We perform this modification in the usrdat2() function (the other functions must be present, but can be empty - they would be used to perform actions at other "insertion" points in NekRS).

c-----------------------------------------------------------------------
      subroutine useric(i,j,k,eg)
      include 'SIZE'
      include 'TOTAL'
      include 'NEKUSE'
      integer i,j,k,e,eg

      return
      end
c-----------------------------------------------------------------------
      subroutine userchk
      include 'SIZE'
      include 'TOTAL'

      return
      end
c-----------------------------------------------------------------------
      subroutine usrdat
      include 'SIZE'
      include 'TOTAL'

      return
      end
c-----------------------------------------------------------------------
      subroutine usrdat2()
      include 'SIZE'
      include 'TOTAL'

      do iel=1,nelt
      do ifc=1,6
      if (boundaryID(ifc,iel).eq.2) then ! rename sideset 2 to 1 (inlet)
          boundaryID(ifc,iel) = 1
       else if (boundaryID(ifc,iel).eq.3) then ! rename sideset 3 to 2 (outlet)
          boundaryID(ifc,iel) = 2
       else if (boundaryID(ifc,iel).eq.4) then ! rename sideset 4 to 3 (pebble surface)
          boundaryID(ifc,iel) = 3
       else if (boundaryID(ifc,iel).eq.5) then ! rename sideset 5 to 4 (side walls)
          boundaryID(ifc,iel) = 4
       else ! NekRS wants all non-sideset faces to be boundary 0
          boundaryID(ifc,iel) = 0
      endif
      enddo
      enddo

      return
      end
c-----------------------------------------------------------------------
      subroutine usrdat3
      include 'SIZE'
      include 'TOTAL'

      return
      end
c-----------------------------------------------------------------------

The various terms in this function are:

• iel is a loop variable name used to indicate a loop over elements

• ifc is a loop variable name used to indicate a loop over the (six) faces on an element

• boundaryID(ifc, iel) is the sideset ID for velocity

Execution and Postprocessing

When you compile Cardinal, you automatically also get a standalone NekRS executable. You can run the NekRS case with

$NEKRS_HOME/bin/nrsmpi pebble 4

where pebble is the case name and 4 is the number of MPI ranks you would like to use. Note that NekRS does not use any shared memory parallelism (e.g. OpenMP). Running this case will produce NekRS output files named pebble0.f<n>, where <n> is a 5-digit integer indicating the output file number (i.e. they will be sequentially-ordered based on when you told NekRS to write output files).

These files (referred to as "field" files) need to be accompanied by a small "configuration" file in order to load into Paraview. To generate that file, we need to use the visnek script like so

cardinal/scripts/visnek pebble

This will create a file named pebble.nek5000, which has very simple contents:

 filetemplate: pebble%01d.f%05d
 firsttimestep: 1
 numtimesteps: 8

This file tells Paraview how many time steps there are to load and the naming pattern for those files. To open the NekRS output files, you then need to open the pebble.nek5000 file.

First, let's take a look at the mesh to get a better sense for how NekRS stores data in its output files. Figure 3 shows a zoomed-in view near the pebble surface. As you can see, this mesh displays the GLL quadrature points. Outlined in red is a single element in our starting Exodus mesh (pebble.exo); the lines inside each element are connecting the nodes.

Figure 3: NekRS spectral mesh, for polynomial order 5

Note that the pebble surface is "faceted" in the same way as our starting Exodus mesh. This is the case because we used a first-order Hex8 mesh.

Figure 4 shows the velocity, pressure, and temperature predicted by NekRS.

Figure 4: NekRS predicted velocity, pressure, and temperature for laminar flow