Chapter 4

4.1 Setting up the geometry

Note that in case of any changes in the SIZE file, a recompilation is necessary.

4.1.1 Uniformly Distributed Mesh

Suppose you wish to simulate flow through an axisymmetric pipe, of radius R = 0.5  and length L = 4  . You estimate that you will need 3 elements in radial (y  ) direction, and 5 in the x  direction, as depicted in Fig. 4.1. This would be specified by the following input file (called to genbox:

2                      spatial dimension  
1                      number of fields  
#    comments:   This is the box immediately behind the  
#                refined cylinder in Ugo’s cyl+b.l. run.  
Box 1                         Pipe  
-5 -3                         Nelx  Nely  
0.0   4.0   1.0               x0  x1   ratio  
0.0   0.5   1.0               y0  y1   ratio  
v  ,O  ,A  ,W  ,   ,          BC’s:  (cbx0, cbx1, cby0, cby1, cbz0, cbz1)


Figure 4.1: Axisymmteric pipe mesh

4.1.2 Graded Mesh


Figure 4.2: Axisymmteric pipe mesh, graded

Suppose you wish to have the mesh be graded, that you have increased resolution near the wall. In this case you change ratio in the y  -specification of the element distribution. For example, changing the 3 lines in the above genbox input file from

-5 -3                         Nelx  Nely  
0.0   4.0   1.0               x0  x1   ratio  
0.0   0.5   1.0               y0  y1   ratio


-5 -4                         Nelx  Nely  
0.0   4.0   1.0               x0  x1   ratio  
0.0   0.5   0.7               y0  y1   ratio

yields the mesh shown in Fig. 4.2.

4.1.3 User-Specified Distribution


Figure 4.3: Axisymmteric pipe mesh, user specified

You can also specify your own, precise, distribution of element locations. For example, another graded mesh similar to the one of the preceding example could be built by changing the genbox input file to contain:

-5  4                                               Nelx  Nely  
0.0   4.0   1.0                                     x0  x1   ratio  
0.000    0.250    0.375    0.450    0.500           y0  y1 ... y4

Here, the positive number of elements for the y  direction indicates that genbox is expecting Nely+1 values of y  positions on the y  -element distribution line. This is the genbox default, which explains why it corresponds to Nely >  0. The corresponding mesh is shown in Fig. 4.3.

4.1.4 Mesh Modification in Nek5000

For complex shapes, it is often convenient to modify the mesh direction in the simulation code, Nek5000. This can be done through the usrdat2 routine provided in the .usr file. The routine usrdat2 is called by nek5000 immediately after the geometry, as specified by the .rea file, is established. Thus, one can use the existing geometry to map to a new geometry of interest.

For example, suppose you want the above pipe geometry to have a sinusoidal wall. Let x := (x,y)  denote the old geometry, and x ′ := (x′,y′)  denote the new geometry. For a domain with y ∈ [0,0.5]  , the following function will map the straight pipe geometry to a wavy wall with amplitude A  , wavelength λ  :

y ′(x,y) = y + yA sin(2πx∕λ).
Note that, as y -→ 0  , the perturbation, yA sin(2 πx∕λ)  , goes to zero. So, near the axis, the mesh recovers its original form.

In nek5000, you would specify this through usrdat2 as follows

      subroutine usrdat2  
      include ’SIZE’  
      include ’TOTAL’  
      real lambda  
      ntot = nx1*ny1*nz1*nelt  
      lambda = 3.  
      A      = 0.1  
      do i=1,ntot  
         argx         = 2*pi*xm1(i,1,1,1)/lambda  
         ym1(i,1,1,1) = ym1(i,1,1,1) + ym1(i,1,1,1)*A*sin(argx)  
      param(59) = 1.  ! Force nek5 to recognize element deformation.  

Note that, since nek5000 is modifying the mesh, postx will not recognize the current mesh unless you tell it to, because postx looks to the .rea file for the mesh geometry. The only way for nek5000 to communicate the new mesh to postx is via the .fld file, so you must request that the geometry be dumped to the .fld file. This is done by modifying the OUTPUT SPECIFICATIONS, which are found near the bottom of the .rea file. Specifically, change

   T      VELOCITY  
   T      PRESSURE  


   T      COORDINATES                       <------  CHANGE HERE  
   T      VELOCITY  
   T      PRESSURE  

The result of above changes is shown in Fig. 4.4.


Figure 4.4: Axisymmteric pipe mesh

4.1.5 Cylindrical/Cartesian-transition Annuli

PIC (a) Annuli mesh PIC (b) Annuli mesh

Figure 4.5: Cylinder mesh

An updated version of genb6, known as genb7, is currently under development and designed to simply/automate the construction of cylindrical annuli, including basic transition-to-Cartesian elements. More sophisticated transition treatments may be generated using the GLOBAL REFINE options in prenek or through an upgrade of genb7, as demand warrants. Example 2D and 3D input files are provided in the nek5000/doc files box7.2d and box7.3d. Figure 4.5a shows a 2D example generated using the box7.2d input file, which reads:

2                      spatial dimension  
1                      number of fields  
#    comments  
Y                   cYlinder  
3 -24 1             nelr,nel_theta,nelz  
.5 .3               x0,y0 - center of cylinder  
ccbb                descriptors: c-cyl, o-oct, b-box (1 character + space)  
.5 .55 .7 .8        r0 r1 ... r_nelr  
0  1  1             theta0/2pi theta1/2pi  ratio  
v  ,W  ,E  ,E  ,    bc’s (3 characters + comma)

An example of a mesh is shown in Fig. 4.5a. The mesh has been quad-refined once with oct-refine option of prenek. The 3D counterpart to this mesh could joined to a hemisphere/Cartesian transition built with the spherical mesh option in prenek.

4.2 Extrusion/Mirroring

4.2.1 Building Extruded Meshes with n2to3

In nek5000/tools, there is a code n2to3.f that can be compiled with your local fortran compiler (preferably not g77). By running this code, you can extend two dimensional domains to three dimensional ones with a user-specified number of levels in the z-direction. Such a mesh can then be modified using the mesh modification approach. Assuming you have a valid two-dimensional mesh, n2to3 is straightforward to run. Below is a typical session, upon typing n2to3 the user is prompted at the command line

 Input old (source) file name:  
 Input new (output) file name:  
 input number of levels: (1, 2, 3,... etc.?):  
 input z min:  
 input z max:  
 input gain (0=custom,1=uniform,other=geometric spacing):  
 This is for CEM: yes or no:  
 Enter Z (5) boundary condition (P,v,O):  
 Enter Z (6) boundary condition (v,O):  
 this is cbz: v  O   <---  
      320 elements written to h3e.rea  

In this context CEM stands for computational electromagnetics, a spin-off of the original Nek5000 code.

The domain in which the fluid flow/heat transfer problem is solved consists of two distinct subdomains. The first subdomain is that part of the region occupied by fluid, denoted Ωf  , while the second subdomain is that part of the region occupied by a solid, denoted Ωs  . These two subdomains are depicted in Fig. 1.1. The entire domain is denoted as D = Ω  ∪ Ω
      f    s  . The fluid problem is solved in the domain Ωf  , while the temperature in the energy equation is solved in the entire domain; the passive scalars can be solved in either the fluid or the entire domain.

We denote the entire boundary of Ωf  as ∂Ωf  , that part of the boundary of Ωf  which is not shared by Ωs  as ----
∂ Ωf  , and that part of the boundary of Ωf  which is shared by Ωs  . In addition,      ----
∂ Ωs,∂Ωs  are analogously defined. These distinct portions of the domain boundary are illustrated in Fig.1.1. The restrictions on the domain for Nek5000 are itemized below.

Fully general three-dimensional meshes generated by other softwares packages can be input to PRENEK as import meshes.

4.3 Moving Geometry

If the imposed boundary conditions allow for motion of the boundary during the solution period (for example, moving walls, free-surfaces, melting fronts, fluid layers), then the geometry of the computational domain is automatically considered in Nek5000 as being time-dependent.

For time-dependent geometry problems, a mesh velocity w is defined at each collocation point of the computational domain (mesh) to characterize the deformation of the mesh. In the solution of the mesh velocity, the value of the mesh velocity at the moving boundaries is first computed using appropriate kinematic conditions (for free-surfaces, moving walls and fluid layers) or dynamic conditions (for melting fronts). On all other external boundaries, the normal mesh velocity on the boundary is always set to zero. In the tangential direction, either a zero tangential velocity condition or a zero tangential traction condition is imposed; this selection is automatically performed by Nek5000 based on the fluid and/or thermal boundary conditions specified on the boundary. However, under special circumstances the user may want to override the defaults set by Nek5000, this is described in the PRENEK manual in Section 5.7.1 If the zero tangential mesh velocity is imposed, then the mesh is fixed in space; if the zero traction condition is imposed, then the mesh can slide along the tangential directions on the boundary. The resulting boundary-value-problem for the mesh velocity is solved in Nek5000 using a elastostatic solver, with the Poisson ratio typically set to zero. The new mesh geometry is then computed by integrating the mesh velocity explicitly in time and updating the nodal coordinates of the collocation points.

Note that the number of macro-elements, the order of the macro-elements and the topology of the mesh remain unchanged even though the geometry is time-dependent. The use of an arbitrary-Lagrangian-Eulerian description in Nek5000 ensures that the moving fronts are tracked with the minimum amount of mesh distortion; in addition, the elastostatic mesh solver can handle moderately large mesh distortion. However, it is the responsibility of the user to decide when a mesh would become "too deformed" and thus requires remeshing. The execution of the program will terminate when the mesh becomes unacceptable, that is, a one-to-one mapping between the physical coordinates and the isoparametric local coordinates for any macro-element no longer exists.

4.4 Boundary and initial conditions

4.4.1 Boundary Conditions

The boundary conditions for the governing equations given in the previous section are now described.

The boundary conditions can be imposed in various ways:

The general convention for boundary conditions in the .rea file is

Since there are no supporting tools that will correctly populate the .rea file with the appropriate values, temperature, velocity, and flux boundary conditions are typically lower case and values must be specified in the userbc subroutine in the .usr file.

4.4.2 Fluid Velocity

Two types of boundary conditions are applicable to the fluid velocity : essential (Dirichlet) boundary condition in which the velocity is specified; natural (Neumann) boundary condition in which the traction is specified. For segments that constitute the boundary ∂Ωf  , see Fig. 1.1, one of these two types of boundary conditions must be assigned to each component of the fluid velocity. The fluid boundary condition can be all Dirichlet if all velocity components of u  are specified; or it can be all Neumann if all traction components t = [- pI + μ(∇u + (∇u )T)]⋅n  , where I  is the identity tensor, n  is the unit normal and μ  is the dynamic viscosity, are specified; or it can be mixed Dirichlet/Neumann if Dirichlet and Neumann conditions are selected for different velocity components. Examples for all Dirichlet, all Neumann and mixed Dirichhlet/Neumann boundaries are wall, free-surface and symmetry, respectively. If the nonstress formulation is selected, then traction is not defined on the boundary. In this case, any Neumann boundary condition imposed must be homogeneous; i.e., equal to zero. In addition, mixed Dirichlet/Neumann boundaries must be aligned with one of the Cartesian axes.

For flow geometry which consists of a periodic repetition of a particular geometric unit, the periodic boundary conditions can be imposed, as illustrated in Fig. 1.1.

Identifier Description Parameters No of Parameters

P periodic periodic element and face 2
V Dirichlet velocity u,v,w 3
O outflow - 0
W wall (no slip) - 0
F flux flux 1
SYM symmetry - 0
A axisymmetric boundary - 0
MS moving boundary - 0
ON Outflow, Normal - 0
E Interior boundary Neighbour element ID 2

Table 4.1: Primitive boundary conditions (flow velocity)

Identifier Description

v user defined Dirichlet velocity
t user defined Dirichlet temperature
f user defined flux

Table 4.2: User defined boundary conditions (flow velocity)

The open(outflow) boundary condition ("O") arises as a natural boundary condition from the variational formulation of Navier Stokes. We identify two situations

For a fully-developed flow in such a configuration, one can effect great computational efficiencies by considering the problem in a single geometric unit (here taken to be of length L), and requiring periodicity of the field variables. Nek5000 requires that the pairs of sides (or faces, in the case of a three-dimensional mesh) identified as periodic be identical (i.e., that the geometry be periodic).

For an axisymmetric flow geometry, the axis boundary condition is provided for boundary segments that lie entirely on the axis of symmetry. This is essentially a symmetry (mixed Dirichlet/Neumann) boundary condition in which the normal velocity and the tangential traction are set to zero.

For free-surface boundary segments, the inhomogeneous traction boundary conditions involve both the surface tension coefficient σ  and the mean curvature of the free surface.

Passive scalars and Temperature

The three types of boundary conditions applicable to the temperature are: essential (Dirichlet) boundary condition in which the temperature is specified; natural (Neumann) boundary condition in which the heat flux is specified; and mixed (Robin) boundary condition in which the heat flux is dependent on the temperature on the boundary. For segments that constitute the boundary    ′     ′
∂Ω f ∪ ∂ Ωs  (refer to Fig. 2.1), one of the above three types of boundary conditions must be assigned to the temperature.

The two types of Robin boundary condition for temperature are : convection boundary conditions for which the heat flux into the domain depends on the heat transfer coefficient hc  and the difference between the environmental temperature T
  ∞ and the surface temperature; and radiation boundary conditions for which the heat flux into the domain depends on the Stefan-Boltzmann constant/view-factor product hrad  and the difference between the fourth power of the environmental temperature T∞ and the fourth power of the surface temperature.

Identifier Description Parameters No of Parameters

T Dirichlet temperature/scalar value 1
O outflow - 0
P periodic boundary - 0
I insulated (zero flux) for temperature 0

Table 4.3: Primitive boundary conditions (Temperature and Passive scalars)

Identifier Description

t user defined Dirichlet temperature
c Newton cooling
f user defined flux

Table 4.4: User defined boundary conditions (Temperature and Passive scalars)

The boundary conditions for the passive scalar fields are analogous to those used for the temperature field. Thus, the temperature boundary condition menu will reappear for each passive scalar field so that the user can specify an independent set of boundary conditions for each passive scalar field.

4.4.3 Internal Boundary Conditions

In the spatial discretization, the entire computational domain is subdivided into macro-elements, the boundary segments shared by any two of these macro-elements in Ωf  and Ωs  are denoted as internal boundaries. For fluid flow analysis with a single-fluid system or heat transfer analysis without change-of-phase, internal boundary conditions are irrelevant as the corresponding field variables on these segments are part of the solution. However, for a multi-fluid system and for heat transfer analysis with change-of-phase, special conditions are required at particular internal boundaries, as described in the following.

For a fluid system composes of multiple immiscible fluids, the boundary (and hence the identity) of each fluid must be tracked, and a jump in the normal traction exists at the fluid-fluid interface if the surface tension coefficient is nonzero. For this purpose, the interface between any two fluids of different identity must be defined as a special type of internal boundary, namely, a fluid layer; and the associated surface tension coefficient also needs to be specified.

In a heat transfer analysis with change-of-phase, Nek5000 assumes that both phases exist at the start of the solution, and that all solid-liquid interfaces are specified as special internal boundaries, namely, the melting fronts. If the fluid flow problem is considered, i.e., the energy equation is solved in conjunction with the momentum and continuity equations, then only the common boundary between the fluid and the solid (i.e., all or portion of ∂Ω-′
   f  in Fig. 1.1) can be defined as the melting front. In this case, segments on  --′
∂Ω f  that belong to the dynamic melting/freezing interface need to be specified by the user. Nek5000 always assumes that the density of the two phases are the same (i.e., no Stefan flow); therefore at the melting front, the boundary condition for the fluid velocity is the same as that for a stationary wall, that is, all velocity components are zero. If no fluid flow is considered, i.e., only the energy equation is solved, then any internal boundary can be defined as a melting front. The temperature boundary condition at the melting front corresponds to a Dirichlet condition; that is, the entire segment maintains a constant temperature equal to the user-specified melting temperature Tmelt  throughout the solution. In addition, the volumetric latent heat of fusion ρL  for the two phases, which is also assumed to be constant, should be specified.

4.4.4 Initial Conditions

For time-dependent problems Nek5000 allows the user to choose among the following types of initial conditions for the velocity, temperature and passive scalars:

A tabulated summary of the compatibility of these initial condition options with various other solution strategies/parameters is given in the appendix.

4.5 Mesh Partioning for Parallel Computing

Genmap is spectral graph partitioning tool, similar to e.g. METIS, which partitions the graph associated to the mesh to assure optimal communication time in HPC applications. Let us consider a simple mesh such as the one in Fig. 4.6. The vertices are distributed in a random fashion, which is the way they may be provided by some mesh generator. Let us assume the vertices are here given as

V1 = (- 1,0), V2 = (0,1), V3 = (- 1,2), V4 = (- 1,1), V5 = (0,2), V6 = (0,0), V7 = (1,1), V8 = (1,0)

The geometry is already stored in the .rea file by the point coordinates, and not vertex numbers

Element  1 = [V1 V6 V2 V4]
x     =
 1,...,4  -1. 0. 0. -1.
y1,...,4 =  0. 0. 1. 1.
Element  2 = [V8 V7 V2 V6]
x1,...,4 =  1. 1. 0. 0.
y1,...,4 =  0. 1. 1. 0.
Element  3 = [V5 V3 V4 V2]
x1,...,4 =  0. -1. -1. 0.
y1,...,4 =  2 . 2. 1. 1.


Figure 4.6: Two-dimensional mesh

Let us a regard the mesh in Fig. 4.6 as a graph of N  vertices and M  edges, G (VN ,EM )  . We define the Laplacian matrix associated to a graph G  as L(G )  . We define as the degree of a node Vi  the number of incident edges, e.g. in Fig. 4.6 deg(V2) = 4  and deg (V6) = 3  .

L(G )  =   i = j    deg(Vi)
     ij     i ⁄= j    - 1 if there is and edge (i,j) and 0 otherwise

       | 1|   2   0    0  - 1   0   - 1  0    0|
       || 2|   0   4    0  - 1  - 1  - 1 - 1   0||
       ||                                       ||
       || 3|   0   0    2  - 1  - 1  0    0    0||
L(G) = | 4|  - 1 - 1  - 1  3    0   0    0    0|
       || 5|   0  - 1  - 1  0    2   0    0    0||
       || 6|  - 1 - 1   0   0    0   3    0   - 1||
       || 7|   0  - 1   0   0    0   0    2   - 1||
       ( 8|   0   0    0   0    0   - 1 - 1   2)

Properties of L(G )

The main ides of the spectral bisection algorithm is

compute \(v_2\) eigenvector corresponding to \(\lambda_2(L(G))\)  
for i=1,N  
  if v_2(i) < 0 put vertex \(V_i\) in N_{-}  
  else put vertex \(V_i\) in N_{-}

The eigenvectors and eigenvalues are computed using Lanczos’s algorithm. These steps are repeated recursively on each of the two branches of the graph N - ,N+   . This is possible since according to Fiedler’s theorems the graph N- is connected, N+   connected only if no v2(i) = 0  , and for each subgraph G1   the algebraic connectivities satisfy λ2(L (G1)) ≤ λ2 (L (G))  .

To run the genmap code be sure that the Nek tools are up-to-date and compiled. At command line type: genmap NOTE-If the executables for the tools were not placed in the bin directory(default), include the path to the genmap executable. We give here the output for the .rea file in the Kovasznay example

Input (.rea) file name:  
Input mesh tolerance (default 0.2):  
NOTE: smaller is better, but generous is more forgiving for bad meshes.  
 reading .rea file data ...  
 start locglob_lexico:           8         960        7680  0.10000000000000001  
 locglob:           1           1        7680  
 locglob:           3        1254        7680  
 done locglob_lexico:        1254        1254        7680           8  
 start periodic vtx:         960        1254  
 done periodic vtx  
 start rec_bisect:         960  
 done:    0.1%  
 done:   99.4%  
 done rec_bisect  

The user is prompted for .rea file name and should enter only the prefix of the .rea file. The user is prompted for mesh tolerance value. Typically a value of .05 is sufficient. Increasing or decreasing this value should make very little difference in the mesh generation. However, if given an error from genmap, the tolerance may need to be made slightly more generous.

A successful genmap run will produce a .map file with the proper processor decomposition.

NOTE: For large element counts, it is not uncommon for genmap to be produce a few disconnected sets. These sets are typically under 7 elements large and will not affect optimization of the NEK5000 run. If a disconnected set is produced, genmap will output the following warning to stdout.

 not connected   N0   NEL  Nsets   Nlarge sets

Here, N0 is the number of elements disconnected from the set of NEL elements, Nsets is the counter of disconnected sets found, and Nlarge sets is the number of sets greater than 64 elements in size. Nlarge sets should always be 0. If not, please contact someone on the developer team so we can be sure to have a more optimal partition of your mesh.

Genmap outputs an ordered set of numbers which are organized as follows Line number 1 contains the header nel, nactive, depth, d2, npts, nrank, noutflow

For the Kovasnay flow on an 8 element mesh with periodic boundary conditions we have 8 12 3 8 32 12 0

Next we have the data (one line per element, listed in order of global element number) ==== 6 12 11 6 5

This means that elemnt one (since we are on the first line) belongs to group 6, and this element is given by vertices in unique ordering. The vertices are ordered in symmetric ordering (starting at 1)

3 - 4 | | 1 - 2

To distribute amongst processors, one just takes as many consecutive processors as one wants.

1This manual is old may soon be deprecated