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 19.
For a more complete list see Commonly Used Variables.
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 22 |
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.
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 21 describes the variables that are assigned in nekasgn
.
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 23 |
utrans |
convective coefficient | vtrans(ix,iy,iz,ie,ifield) |
See Table 23 |
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 22.
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.
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.
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
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 24 along with their definitions and the corresponding entry in the cbc
array, where \(\mathbf{\hat e}\) denotes a unit vector.
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 |
pa |
pressure | \(p\) | o , on |
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 25 along with their definitions and the corresponding entry in the cbc
array.
These correspond to standard Dirichlet, Neumann, and Robin boundary conditions.
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_c\) | \(\lambda\nabla T\cdot\mathbf{\hat e_n}=h_c(T-T_{\infty})\) | c |
tinf |
ambient temperature, \(T_{\infty}\) | ||
hrad |
heat transfer coefficient, \(h_{rad}\) | \(\lambda\nabla T\cdot\mathbf{\hat e_n}=h_{rad}(T^4-T^4_{\infty})\) | r |
tinf |
ambient temperature, \(T_{\infty}\) |
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. |
---|
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 26 along with their definitions and the corresponding entry in the cbc
array.
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}\) |
Note
Temperature and passive scalar boundary conditions are contextual and will set the boundary condition for the active solution field. See Table 22.
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.
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.