User Routines File (.usr)¶

The user file is a special case file that implements the user interface to Nek5000. What follows is a brief description of the available subroutines and a primer on writing your own user code. An empty user file template can be found in the source directory at Nek5000/core/zero.usr. When setting up a new case, it is recommended to start by copying this file to your run directory when you don’t already have a similar case to use as a template.

To aid in implementing custom physics, the Nek5000 user file is partitioned into subroutines which provide access to specific, commonly used terms in the solver. These subroutines can be broadly divided into 3 categories:

Local routines - called every time step for every field and GLL point
Global routines - called once per time step
Initialization routines - only called once during initialization

All of the routines in the .usr file have access to a common set of global variables and arrays. These are included in the TOTAL common block, included at the beginning of each subroutine. Some of the more commonly used variables are shown in Table 18. For a more complete list see Commonly Used Variables.

Table 18 Commonly used global variables available in `TOTAL`¶
Variable	Description	Note
`pi`	\(\pi\)	`pi=4.0atan(1.0)`
`time`	physical time
`nelv`	number of elements in the velocity mesh on the local MPI rank
`nelt`	number of elements in the temperature mesh on the local MPI rank	see Conjugate Heat Transfer
`nelgv`	total number of elements in the velocity mesh
`nelgt`	total number of elements in the temperature mesh	see Conjugate Heat Transfer
`idsess`	session ID for NekNek	see Overlapping Overset Grids
`ifield`	active solution field	see Table 21

Note

nelt is always greater than or equal to nelv. For non-conjugate heat transfer cases, nelv and nelt are identical. Similarly for nelgv and nelgt.

Also included in the TOTAL block is the field coefficient array. This array contains reference values for fluid properties that are inherited directly from the .par file. For a simulation with constant properties, this array corresponds directly to density, viscosity, etc. as the values from this array are copied directly into the vtrans and vdiff arrays internally by Nek5000 (See uservp for more information on vtrans and vdiff). For simulations with variable properties, this array can be useful as it retains the values assigned in the .par file.

Table 19 The field coefficient array¶
Variable	Entry in `.par` file	Description
`cpfld(1,1)`	`velocity:viscosity`	Reference viscosity
`cpfld(1,2)`	`velocity:density`	Reference density
`cpfld(2,1)`	`temperature:conductivity`	Reference conductivity
`cpfld(2,2)`	`temperature:rhocp`	Reference rho-cp
`cpfld(3,1)`	`scalar01:diffusivity`	Reference diffusivity for scalar 1
`cpfld(3,2)`	`scalar01:density`	Reference density for scalar 1
`cpfld(i+2,1)`	`scalar i:diffusivity`	Reference diffusivity for \(i\)
`cpfld(i+2,2)`	`scalar i:density`	Reference density for scalar \(i\)

Note

The entries in the field coefficient array for velocity and temperature correspond directly to parameters 1, 2, 7, and 8 from the old .rea format. However, no corresponding parameters exist for the passive scalars.

Local Routines¶

The local subroutines uservp, userf, and userq are called at every GLL point, for every field at every time step. The subroutines userbc and useric are similar, but userbc is only called for GLL points on a boundary face with certain boundary conditions and useric is only called during initialization. These subroutines take ix, iy, iz, and eg as arguments, which correspond to the local GLL indexing and global element number. Note that the local element number can be accessed via the gllel array.

NEKUSE¶

The ix, iy, etc. indices can be used to cross-reference the local solution arrays, however, the NEKUSE common block is provided for easier access. This block contains solver variables that may be useful for defining custom models, such as variable properties, or a localized heating rate. Immediately before a local routine is called by Nek5000, the subroutine nekasgn is called. This routine sets many of the commonly used values in NEKUSE, making them available for use in the local subroutines. Table 20 describes the variables that are assigned in nekasgn.

Table 20 Prepopulated `NEKUSE` variables¶
Variable	Description	Solution Array	Note
`x`	x-coordinate	`xm1(ix,iy,iz,ie)`
`y`	y-coordinate	`ym1(ix,iy,iz,ie)`
`z`	z-coordinate	`zm1(ix,iy,iz,ie)`
`r`	r-coordinate	N/A	\(\sqrt{x^2+y^2}\)
`theta`	\(\theta\)-coordinate	N/A	`theta=atan2(y,x)`
`ux`	x-velocity	`vx(ix,iy,iz,ie)`
`uy`	y-velocity	`vy(ix,iy,iz,ie)`
`uz`	z-velocity	`vz(ix,iy,iz,ie)`
`temp`	temperature	`t(ix,iy,iz,ie,1)`
`ps(i)`	passive scalar ‘i’, \(\phi_i\)	`t(ix,iy,iz,ie,i+1)`	\(i=(\) `ifield` \(-2)\)
`pa`	pressure	`pr(ix,iy,iz,ie)`	\(P_N/P_N\) only
`p0`	thermodynamic pressure	`p0th`
`udiff`	diffusion coeffcient	`vdiff(ix,iy,iz,ie,ifield)`	See Table 22
`utrans`	convective coefficient	`vtrans(ix,iy,iz,ie,ifield)`	See Table 22

The active solution field – velocity, temperature, etc. – is indicated by the global variable ifield. The value that corresponds to each field is described in Table 21. ifield is not explicitly passed as an input to any of the user subroutines, rather it is controlled at a higher level directly by Nek5000. It is set by the solver for all of the local routines, but not any of the initialization or global routines.

Table 21 Corresponding solution fields for `ifield`¶
`ifield` value	Field name	Solution variable
1	velocity	\(\mathbf u\)
2	temperature	\(T\)
3	passive scalar 1	\(\phi_1\)
4	passive scalar 2	\(\phi_2\)
\(\ge 5\)	passive scalar \((\) `ifield` \(-2)\)	\(\phi_{ifld-2}\)

uservp¶

This function is only called if the variableProperties key under [PROBLEMTYPE] in the .par file is set to true (see here). It can be used to specify customized or solution dependent material properties. It is called for every GLL-point for every field at every time step. The diffusion and transport coefficients should be set using the variables described in the table below. The diffusion coefficient refers the viscosity, thermal conductivity, or diffusivity for passive scalars. The transport coefficient refers to the coefficient attached to the convective term, typically density for velocity and passive scalars and the product of density and specific heat for temperature.

Table 22 Terms set in `uservp`¶
	Description	Variable in governing equations	solution field(s)
`udiff`	diffusion coefficient	\(\mu\) in the momentum equation	`ifield = 1`
		\(\lambda\) in the energy equation	`ifield = 2`
		\(\Gamma_i\) in the passive scalar equations	`ifield = 3 .. npscal+2`
`utrans`	transport coefficient	\(\rho\) in the momentum equation	`ifield = 1`
		\((\rho c_p)\) in the energy equation	`ifield = 2`
		\(\rho_i\) in the passive scalar equations	`ifield = 3 .. npscal+2`

Warning

The corresponding entries in vdiff and vtrans are overwritten by whatever is assigned to udiff and utrans. Setting vdiff and vtrans directly is not supported.

Example:	The code block below shows how to implement a variable viscosity as a function of temperature, with the density, rho-cp, and thermal conductivity set from the values in the `.par` file using the field coefficient array.

      subroutine uservp(ix,iy,iz,eg) ! set variable properties

      implicit none

      integer ix,iy,iz,eg

      include 'SIZE'
      include 'TOTAL'
      include 'NEKUSE'

      integer e
      real a,b
c     e = gllel(eg)

      a = 0.001 
      b = 10.0

      if (ifield.eq.1) then
        udiff  = a * exp(-b*temp) ! dynamic viscosity
        utrans = cpfld(ifield,2)  ! density
      else
        udiff  = cpfld(ifield,1)  ! conductivity
        utrans = cpfld(ifield,2)  ! rho*cp
      endif

      return
      end

userf¶

This functions sets the source term (which will be subsequently be multiplied by the density) for the momentum equation. It allows the user to effectively add an acceleration term.

Example:	The code block below shows how to implement a body force proportional to the temperature, similar to what would be done for a Boussinesq model for buoyancy.

      subroutine userf(ix,iy,iz,eg) ! set acceleration term
c
c     Note: this is an acceleration term, NOT a force!
c     Thus, ffx will subsequently be multiplied by rho(x,t).
c
      implicit none

      integer ix,iy,iz,eg

      include 'SIZE'
      include 'TOTAL'
      include 'NEKUSE'

      integer e
      real beta
c     e = gllel(eg)

      beta = 1.0

      ffx = 0.0
      ffy = 0.0
      ffz = beta * temp

      return
      end

userq¶

This functions sets the source term for the energy and passive scalar equations. An explicit source term can be set using qvol. In the latest version available from the master branch on github, an implicit source term can be set using avol.

A source term that has the form

\[q'''=\alpha -\beta T\]

can be implemented either entirely explicitly, or semi-implicitly. In general, the implicit term should be used wherever possible as it tends to stabilize the solution. Both approaches are shown below. It is not necessary for \(\alpha\) and \(\beta\) to be constants. They can vary with time, position, or any of the solution variables – including temperature. However, using a solution variable may impose limits on the stability of the solution.

Example:	In the first example, the source term is set entirely explicitly

      subroutine userq(ix,iy,iz,eg) ! set source term

      implicit none

      integer ix,iy,iz,eg

      include 'SIZE'
      include 'TOTAL'
      include 'NEKUSE'

      integer e
      real alpha, beta
c     e = gllel(eg)

      alpha = 1.0
      beta = 10.0

      qvol = alpha - beta*temp !set the explicit source term

      return
      end

Example:	In the second example, the implicit source term is leveraged. Both implementations will yield the same solution for a converged result.

      subroutine userq(ix,iy,iz,eg) ! set source term

      implicit none

      integer ix,iy,iz,eg

      include 'SIZE'
      include 'TOTAL'
      include 'NEKUSE'

      integer e
      real alpha, beta
c     e = gllel(eg)

      alpha = 1.0
      beta = 10.0

      qvol = alpha !set the explicit source term
      avol = beta  !set the implicit source term

      return
      end

userbc¶

This functions sets boundary conditions. It is only called for special boundary condition types – any lowercase value in the cbc array – and only for points on the boundary surface. It includes an additional argument compared to the other Local Routines. The iside variable refers to which side of the element the boundary condition is on. This can be used for accessing the appropriate entry in the boundaryID or cbc arrays. For a more complete list of all supported boundary conditions, see Boundary Conditions.

The available values that can be set for velocity are listed in Table 23 along with their definitions and the corresponding entry in the cbc array, where \(\mathbf{\hat e}\) denotes a unit vector.

Table 23 Velocity boundary conditions set in `userbc`¶
	Description	Definition	`cbc` value
`ux`	x-velocity	\(\mathbf u\cdot\mathbf{\hat e_x}\)	`v`
`uy`	y-velocity	\(\mathbf u\cdot\mathbf{\hat e_y}\)	`v`
`uz`	z-velocity	\(\mathbf u\cdot\mathbf{\hat e_z}\)	`v`
`un`	velocity normal to the boundary face	\(\mathbf u\cdot\mathbf {\hat e_n}\)	`vl`
`u1`	velocity tangent* to the boundary face	\(\mathbf u\cdot\mathbf {\hat e_t}\)	`vl`
`u2`	velocity bitangent* to boundary face	\(\mathbf u\cdot\mathbf {\hat e_b}\)	`vl`
`trx`	traction in the x-direction	\((\boldsymbol{\underline\tau}\cdot\mathbf{\hat e_n})\cdot\mathbf{\hat e_x}\)	`s`, `sh`
`try`	traction in the y-direction	\((\boldsymbol{\underline\tau}\cdot\mathbf{\hat e_n})\cdot\mathbf{\hat e_y}\)	`s`, `sh`
`trz`	traction in the z-direction	\((\boldsymbol{\underline\tau}\cdot\mathbf{\hat e_n})\cdot\mathbf{\hat e_z}\)	`s`, `sh`
`trn`	traction normal to the boundary face	\((\boldsymbol{\underline\tau}\cdot\mathbf{\hat e_n})\cdot\mathbf{\hat e_n}\)	`sl`
`tr1`	traction tangent* to the boundary face	\((\boldsymbol{\underline\tau}\cdot\mathbf{\hat e_n})\cdot\mathbf{\hat e_t}\)	`sl`, `shl`
`tr2`	traction bitangent* to the boundary face	\((\boldsymbol{\underline\tau}\cdot\mathbf{\hat e_n})\cdot\mathbf{\hat e_b}\)	`sl`, `shl`

Warning

*The tangent and bitangent directions are not guaranteed to be consistent between elements in 3D domains.

The available values that can be set for temperature are listed in Table 24 along with their definitions and the corresponding entry in the cbc array. These correspond to standard Dirichlet, Neumann, and Robin boundary conditions.

Table 24 Temperature boundary conditions set in `userbc`¶
	Description	Definition	`cbc` value
`temp`	temperature	\(T\)	`t`
`flux`	heat flux, \(q''\)	\(\lambda\nabla T\cdot\mathbf{\hat e_n}=q''\)	`f`
`hc`	heat transfer coefficient, \(h\)	\(\lambda\nabla T\cdot\mathbf{\hat e_n}=h(T-T_{\infty})\)	`c`
`tinf`	ambient temperature, \(T_{\infty}\)	\(\lambda\nabla T\cdot\mathbf{\hat e_n}=h(T-T_{\infty})\)	`c`

Note

Both heat transfer coefficient and ambient temperature must be specified for a Robin boundary condition.

A few examples are shown next to demonstrate how to set simple boundary conditions.

Example: In the example below, the code sets a parabolic inlet velocity with a constant inlet temperature of zero and a constant wall temperature of one. The temperature field has the same BC of t on both the inlet and the wall, so the velocity BC is accessed to differentiate between the two. Also note that this routine will not be called for ifield=1 for the W boundary, but it will be called for ifield=2 for the t boundary co-located with the W boundary.

      subroutine userbc(ix,iy,iz,iside,eg) ! set up boundary conditions
c
c     NOTE ::: This subroutine MAY NOT be called by every process
c
      implicit none

      integer ix,iy,iz,iside,eg

      include 'SIZE'
      include 'TOTAL'
      include 'NEKUSE'

      integer ie
      character*3 cb3

      ie=gllel(eg) !get local element number
      cb3=cbc(iside,ie,1) !access the velocity boundary condition

      ux   = 0.0
      uy   = 0.0
      uz = 3./2. (1.0-(2.0*y-1.0)**2

      if(cb3.eq.'v  ')
        temp = 0.0 !set inlet temperature to 0
      elseif(cb3.eq.'W  ')
        temp = 1.0 !set wall temperature to 1
      endif
      
      return
      end

Boundary conditions are only applied to boundary faces with the corresponding cbc array value. In the example above, the cbc array does not need to be checked to assign the parabolic velocity profile as it will be ignored for any BC that is not v.

Example:	In this example, the `boundaryID` array is used to differentiate between the inlet and two different walls. The inlet (ID = 1) has a velocity profile and constant temperature value. The walls (IDs 2 and 3 respectively) are set as a positive heat flux on wall 2 and a negative (cooling) heat flux on wall 3. This example corresponds to the example setup shown in usrdat.

      subroutine userbc(ix,iy,iz,iside,eg) ! set up boundary conditions
c
c     NOTE ::: This subroutine MAY NOT be called by every process
c
      implicit none

      integer ix,iy,iz,iside,eg

      include 'SIZE'
      include 'TOTAL'
      include 'NEKUSE'

      integer ie

      ie=gllel(eg)  !get local element number

      ux   = 0.0
      uy   = 0.0
      uz   = 0.0
      temp = 0.0

      if(boundaryID(iside,ie).eq.1)
        uz = 3./2.*(1.0-(2.0*y-1.0)**2)
        temp = 1.0
      elseif(boundaryID(iside,ie).eq.2)
        flux = 1.0
      elseif(boundaryID(iside,ie).eq.3)
        flux = -1.0
      endif

      return
      end

The temperature and passive scalars use the same form of governing equation and are handled practically identically by Nek5000. The boundary conditions available for passive scalars are therefore the same as those available for temperature. These are listed in Table 25 along with their definitions and the corresponding entry in the cbc array.

Table 25 Passive scalar boundary conditions set in `userbc`¶
	Description	Definition	`cbc` value
`temp`	Dirichlet	\(\phi_i\)	`t`
`flux`	Neumann, \(\psi_i\)	\(\Gamma_i\nabla \phi_i\cdot\mathbf{\hat e_n}=\psi_i\)	`f`
`hc`	transfer coefficient, \(\eta\)	\(\Gamma_i\nabla \phi_i\cdot\mathbf{\hat e_n}=\eta(\phi_i-\phi_{i,\infty})\)	`c`
`tinf`	ambient value, \(\phi_{i,\infty}\)		`c`

Note

Temperature and passive scalar boundary conditions are contextual and will set the boundary condition for the active solution field. See Table 21.

Warning

The ps(i) variable array provided by NEKUSE is NOT used to set the passive scalar boundary conditions.

Example:	In this example, the `ifield` variable is used to distinguish between a constant inlet temperature of 0 and a transient concentration for passive scalar 1 which is distributed normally in time.

      subroutine userbc(ix,iy,iz,iside,eg) ! set up boundary conditions
c
c     NOTE ::: This subroutine MAY NOT be called by every process
c
      implicit none

      integer ix,iy,iz,iside,eg

      include 'SIZE'
      include 'TOTAL'
      include 'NEKUSE'

      real time0,sigma

      time0 = 1.0
      sigma = 0.1

      ux   = 0.0
      uy   = 0.0
      uz   = 0.0

      if(ifield.eq.2) then
        temp = 0.0 !set temperature to zero
      elseif(ifield.eq.3) then
        temp = 1./(sigma*sqrt(2.*pi))*exp(-0.5*((time-time0)/sigma)**2) !set scalar 1
      endif

      return
      end

It is not uncommon to use multiple methods of discerning between different boundary conditions. The user may need to implement specific conditions for different fields on different boundaries at different times. This can quickly lead to complex, embedded if-then statements. Use care to ensure you’re setting the correct boundary conditions!

useric¶

This functions sets the initial conditions and behaves similarly to userbc. It is called only during initialization after usrdat3 for every solution field on every GLL point in the domain.

Warning

useric is NOT called for any fields loaded from a restart file.

Global Routines¶

userchk¶

This is a general purpose routine that is executed both during initialization and after every time step. It can be used for a wide variety of purposes, such as monitoring the solution, post-processing results, or implementing entirely new physical solvers.

It is helpful for users to familiarize themselves with the included utility subroutines in Nek5000 as they can make performing complex calculations on the entire dataset much easier. For a list of some of the commonly used subroutines, see here.

Example: Most of the code shown below will only be executed if Nek5000 is run in “post-processing” mode, i.e. numSteps = 0 is specified in the .par file. The solution is rescaled by reference length, velocity, and pressure and written to an output file. The name of the output file will be prepended with the characters dim and it will contain dimensional quantities that can be loaded into Paraview or VisIt for further post processing.

      subroutine userchk()

      implicit none

      include 'SIZE'
      include 'TOTAL'

      integer n1,n2
      real Lref,Uref,rhoref,Pref

      n1=lx1*ly1*lz1*nelv !number of GLL points in the velocity mesh on the local MPI rank
      n2=lx2*ly2*lz2*nelv !number of GLL points in the pressure mesh on the local MPI rank

      Lref   = 0.00915          !hydraulic diameter in meters
      Uref   = 0.768            !reference velocity in m/s
      rhoref = 998.6            !reference density in kg/m3
      Pref   = rhoref*Uref*Uref !reference pressure in Pa

      if(nsteps.eq.0) then ! postprocessing mode

        call cmult(xm1,Lref,n1) !re-dimesionalize the domain
        call cmult(ym1,Lref,n1)
        call cmult(zm1,Lref,n1)
        call cmult(vx,Uref,n1)  !re-dimensionalize the velocity
        call cmult(vy,Uref,n1)
        call cmult(vz,Uref,n1)
        call cmult(pr,Pref,n2)  !re-dimensionalize the pressure

        call prepost(.true.'dim') !write a restart file with dimensional quantities

        call exitt !exit without writing another new restart file

      endif

      return
      end

userqtl¶

This function is only called if the equation key under [PROBLEMTYPE] in the .par file is set to lowMachNS (see here). This function can be used to specify a customized thermal divergence for the low Mach solver. The thermal divergence refers to a non-zero right hand side of the divergence constraint (see Low-Mach Navier-Stokes).

Nek5000 includes a thermal divergence model for a single component ideal gas as denoted by the call to userqtl_scig. To use this model, simply leave userqtl as it is in the template file, zero.usr. The implementation of any other model is left to the user.

Initialization Routines¶

usrdat¶

This function can be used to modify the element vertices and is called before the spectral element mesh (GLL points) has been laid out. It can be used to fill the cbc array based on BoundaryID for 3rd party meshes.

Example:	In the code below, the `cbc` array is filled for a 3rd party mesh. The `boundaryID` array is filled with the Boundary ID values set in Gmsh or the sideset numbers specified in an exodus mesh file. This example corresponds with the setup shown in the second example for userbc.

      subroutine usrdat()   ! This routine to modify element vertices

      implicit none

      include 'SIZE'
      include 'TOTAL'

      integer ie,iside

      do ie=1,nelv !loop over all elements on the local MPI rank
        do iside=1,2*ldim !loop over all sides on the element
          cbc(iside,ie,1)="E  " !first set all faces to internal boundaries
          cbc(iside,ie,2)="E  " !"E  " and "   " both work for this
          if(boundaryID(iside,ie).eq.1) then !boundary 1 is the inlet
            cbc(iside,ie,1) = "v  "
            cbc(iside,ie,2) = "t  "
          elseif(boundaryID(iside,ie).eq.2) then !boundary 2 is a heated wall
            cbc(iside,ie,1) = "W  "
            cbc(iside,ie,2) = "f  "
          elseif(boundaryID(iside,ie).eq.3) then !boundary 3 is another heated wall
            cbc(iside,ie,1) = "W  "
            cbc(iside,ie,2) = "f  "
          elseif(boundaryID(iside,ie).eq.4) then !boundary 4 is the outlet
            cbc(iside,ie,1) = "O  "
            cbc(iside,ie,2) = "I  "
          endif 
        enddo
      enddo

      return
      end

usrdat2¶

This function can be used to modify the spectral element mesh. The geometry information (mass matrix, surface normals, etc.) will be rebuilt after this routine is called. Any changes to the cbc array must be made before or during this call.

Example:	In the code below, the mesh is scaled by a factor of \(1/D_h\), where \(D_h\) is the hydraulic diameter. It is typically convenient to generate a mesh with dimensions, this allows the user to non-dimensionalize the mesh at runtime.

      subroutine usrdat2()   ! This routine to modify mesh coordinates

      implicit none

      include 'SIZE'
      include 'TOTAL'

      integer n
      real Dhyd

      n=lx1*ly1*lz1*nelv !number of GLL points on the local MPI rank

      Dhyd = 0.00915 !hydraulic diameter in meters

      call cmult(xm1,1.0/Dhyd,n)
      call cmult(ym1,1.0/Dhyd,n)
      call cmult(zm1,1.0/Dhyd,n)

      return
      end

usrdat3¶

This function can be used to initialize any additional case/user specific data. The GLL mesh should not be further modified in this routine as the geometric factors are not automatically recomputed.