- Generated on Fri Nov 6 2026 03:35:31 for Neko by
1.9.8
|
Neko 1.99.9
A portable framework for high-order spectral element flow simulations
|
The configuration of the fluid solver and the flow problem. Contains multiple subobjects for various parts of the setup.
As per the governing equations, Neko requires the value of the density and dynamic viscosity to define the flow problem. These can be provided as rho and mu in the case file.
Alternatively, one may opt to provide the Reynolds number, Re, which corresponds to a non-dimensional formulation of the Navier-Stokes equations. This formulation can effectively be obtained by setting \( \rho = 1 \) and \(
\mu = 1/Re \). This is exactly what Neko does under the hood, when Re is provided in the case file.
Note that if both Re and any of the dimensional material properties are provided, the simulation will issue an error.
As an alternative to providing material properties in the case file, it is possible to do that in a special routine in the user file. This is demonstrated in the rayleigh_benard_cylinder example. Ultimately, both rho and mu have to be set in the subroutine. Additionally, this allows to change the material properties in time. Yet another options is to directly manipulate the case file programmatically in the user_startup routine and inject the material properties there. This is demonstrated in the rayleigh_benard example.
When material properties are constant or only vary in time, one can use the simplified form of the viscous stress tensor in the governing equations. However, when there are spatial variations, it is necessary to use the general (full) form. The variation may come, for example, due to a turbulence model, the modifications in the above-mentioned user routine. The general form of the stress tensor requires solving the 3 equations for the velocity components in a coupled manner, which requires an appropriate linear solver. By default, Neko will use the simplified form of the tensor, and the full one must be selected by the user by setting full_stress_formulation to true.
Spectral vanishing viscosity (SVV) selectively adds dissipation to the high-frequency content of the solution. It was first proposed by Tadmor (1989) and later introduced to the spectral element method (SEM) community by Kirby and Sherwin (2006). Neko currently supports only the one-sided formulation, which applies the high-pass operator to the trial-function gradient in physical space (i.e., to gradients in the x, y, and z directions). The full-stress formulation is also supported for velocity and requires coupled_cg (or fused_coupled_cg on CUDA/HIP). For fluid solves, SVV is available with the implicit pnpn scheme.
For the fluid equations, add svv to the fluid object. For a scalar, add the same object directly to that scalar's configuration:
The optional formulation entry defaults to one-sided; it is shown above to make the operator choice explicit. The required kernel object selects the modal transfer function through its type; currently, only power is supported. For the power kernel, the required power_coefficient controls the modal transfer function; larger values confine the added dissipation to modes nearer the polynomial cut-off. The direction selects the reference-element directions in which the modal filter is applied and defaults to rst. The SVV viscosity nu is multiplied by density internally and may be either a constant value or a registered field. A field configuration uses field_name:
For a field-valued viscosity, the optional time_variable entry controls whether the field is refreshed at every time step and defaults to true.
The SVV operator, including its full-stress variant, is implemented for CPU, CUDA, and HIP backends. It is not currently available with the SX, XSMM, OpenCL, or Metal backends.
Neko currently provides several LES models via the les_model simulation component. The simcomp computes the viscosity and stores in the field registry under the name selected in the case file. For more details, see the documentation of the simcomp. To enable LES in the fluid solver, one simply has to add the nut_field keyword to the configuration and set it to the name of the field generated by the les_model simcomp. This will automatically change the governing equations to feature the full viscous stress tensor, as required for a variable viscosity field.
Note that the full viscous stress tensor requires the equations for the 3 velocity components to be solved in a coupled manner. Therefore, coupled_bicgstab, coupled_cg, or fused_coupled_cg should be used for velocity. The coupled_bicgstab solver is currently available on CPUs only.
This feature is enabled by setting the schwarz_iterations keyword inside the fluid group to an integer larger than zero. In this case, each fluid timestep solves for velocity and pressure multiple times.
The total number of passes is 1 + schwarz_iterations. This feature is often needed to increase the stability of solutions when using overset_interface boundary conditions, i.e., when there are multiple coupled simulations running in tandem.
t and tstep variables are only updated across real timesteps. Therefore, if your user conditions depend on these variables, they remain valid.The optional boundary_conditions keyword can be used to specify boundary conditions. The reason for it being optional, is that periodic boundary conditions are built into the definition of the mesh, so for a periodic box nothings needs to be added to the case file. The TGV example is such a case, for instance. The value of the keyword is an array of JSON objects, each specifying a single boundary condition.
In Neko we usually refer to boundaries as "zones", which in this case are face zones, i.e. a collection of element faces. Which zones the boundary condition is applied to is controlled by the zone_indices keyword, which takes an array of integers. It is up to the user whether to apply a single condition to multiple zones or specify several conditions applied to one zone each. For example, if you have two zones, which should be no-slip walls, you can either create two no_slip conditions, one for each zone, or just create a single condition and apply it to both.
The indices your boundaries have is determined by the mesh. To check them, you can use the mesh_checker utility with the optional --write_zone_indices argument. This will output a zone_indices0.f00000 file that you can inspect in Paraview, and the boundaries will be marked by their index value.
Recall that periodic conditions are built into the mesh, since they are topological in nature. This means that you must not specify any conditions for the corresponding zones. For example, in the turb_pipe example, which is a periodic pipe simulation, two periodic zones comprise the boundary conditions in the streamwise direction. Only one condition, corresponding to zone index 3 (the wall) is the specified in the case file.
It is possible to assign specific names to the boundary conditions through the name keyword. Boundary conditions can then be retireved in the code by using the name or the zone_index where it is applied.
The default name of the boundary conditions is given by the <variable>_bc_<zone_index> pattern. i.e., the pressure boundary condition that applies in zone index 5 can be retrieved by the pressure_bc_5 name.
The conditions to apply are specified by the type keyword inside each of the JSON objects. Each type sets both the velocity and the pressure condition on the zone, and the table below lists what each type corresponds to. An empty entry means the natural condition of the scheme for that field: a homogeneous Neumann condition for the velocity, or a Neumann condition for the pressure with the normal gradient given by the momentum equation. These combinations cover the standard cases. It is also possible to set the velocity and the pressure condition separately, see below.
| Boundary Condition | Velocity | Pressure | Description |
|---|---|---|---|
| symmetry | symmetry | A symmetry plane. Must be axis-aligned. | |
| velocity_value | dirichlet | A Dirichlet condition for velocity. | |
| expression_velocity | expression | A Dirichlet condition for velocity given by mathematical expressions. | |
| expression_pressure | expression | A Dirichlet condition for pressure given by a mathematical expression. | |
| no_slip | no_slip | A no-slip wall. Either stationary or moving. | |
| outflow | zero_dirichlet | A pressure outlet. | |
| normal_outflow | non_normal | zero_dirichlet | An Neumann condition for the surface-normal component of velocity combined with a Dirichlet for the surface-parallel components. |
| outflow+user | user_dirichlet | Same as outflow but with user-specified pressure. | |
| normal_outflow+user | non_normal | user_dirichlet | Same as normal_outflow but with user-specified pressure. |
| outflow+dong | dong | A pressure outlet with the Dong condition applied. | |
| normal_outflow+dong | non_normal | dong | The normal_outflow with the Dong condition applied. |
| shear_stress | shear_stress | Prescribed wall shear stress. | |
| wall_model | wall_model | Shear stress condition based on a wall model for large-eddy simulation. | |
| blasius_profile | blasius_profile | A Blasius velocity profile. | |
| user_velocity | user_dirichlet | The field_dirichlet_vector_t user-defined Dirichlet condition for velocity. | |
| user_pressure | user_dirichlet | The field_dirichlet_t user-defined Dirichlet condition for pressure. | |
| overset_interface | overset_interface | overset_interface (if couple_pressure) | A Dirichlet condition that prescribes values from another Neko simulation running concurrently. |
A more detailed description of each boundary condition is provided below.
symmetry. A symmetry plane. Sets the surface-normal velocity to 0 and applies a homogenous Neumann condition to the surface-parallel components. Requires no additional keywords. velocity_value. A Dirichlet condition for velocity. Suitable for velocity inlets, moving walls, certain freestream conditions, etc. The value is prescribed by the value keyword, that should be an array of 3 reals.
expression_velocity. A Dirichlet condition for velocity, where the value keyword holds an array of three strings, one mathematical expression per component. See expressions for the syntax. This covers inlet profiles and time varying inlets without a user file.
The expressions are evaluated in the boundary points only. An expression that uses neither t nor dt is evaluated once, at setup, and then costs nothing more than velocity_value at runtime; one that uses either of them is re-evaluated once per timestep.
expression_pressure. The same, for pressure, where value is a single string.
no_slip. A standard no-slip wall, which sets velocity to zero relative to the wall. For moving walls, setting the optional argument "moving": true is required. This also requires setting up the ALE module separately. further details can be found in ALE user guide. For stationary walls, no additional keyword is needed. outflow. A standard pressure outlet condition. Requires no additional keywords. normal_outflow. The condition lets the flow escape through the boundary by setting a homogeneous Neumann condition for the surface-normal velocity component, but fixes the values of the surface-parallel components. The latter values are prescribed in the boundary condition's JSON. There are two alternative ways to do the tangential components. One is to simply provide a value vector. This vector is defined in global Cartesian coordinates and internally a projection onto the local basis is performed. Afterwards, the local tangential values are enforced. This works both in the axis-aligned and in the fully general mixed-boundary implementation.
The second option is to provide the values through a field file using the file_name keyword. This is currently supported only by the axis-aligned implementation. The file must contain a three-component vector field in global Cartesian coordinates. Internally, only the tangential components are enforced on the boundary. The file_name and mesh_file_name values use Neko's field-sample naming convention, for example field0.f00000. The mesh_file_name is an optional field sample containing the source mesh coordinates; it can be omitted when the source field contains its own coordinates. Optional interpolation settings are the same as for other field-imported data:
interpolate. Logical flag controlling whether interpolation is used.mesh_file_name. Optional sampled field containing the source mesh coordinates when the source field is defined on a different mesh.interpolation.tolerance. Tolerance for the interpolation search.interpolation.padding. Padding used in the interpolation search.Exactly one of value or file_name must be provided. If value is used, it must be an array of three reals.
outflow+user. Same as outflow, but with user-specified pressure. The pressure is specified via the same interface as user_pressure, see the relevant section for more information.normal_outflow+user. Same as normal_outflow, but with user-specified pressure. The pressure profile is specified via the same interface as user_pressure, see the relevant section for more information. The tangential velocity values are prescribed exactly as for normal_outflow, i.e. using either value or file_name.
outflow+dong. Same as outflow, but additionally applies the Dong boundary condition on the pressure. This is a way to prevent backflow and therefore promotes numerical stability without having to create a mesh "sponge" at the outlet.normal_outflow+dong. Same as normal_outflow, but additionally applies the Dong boundary condition for the pressure to prevent backflow. The tangential velocity values are prescribed exactly as for normal_outflow, i.e. using either value or file_name.
shear_stress. Non-penetration condition combined with a set shear stress vector. Only works with axis-aligned boundaries. The stress value is specified by the value keyword, which should be an array of 3 reals. It is the responsibility of the user to set the vector in the direction parallel to the boundary. wall_model. A shear stress condition where the values are computed by a wall model. Meant to be used for wall-modelled large-eddy simulation. Only works with axis-aligned boundaries. The model is selected using the model keyword. Additional configuration depends on the model selected.spaldingmodel requires specifying kappa and B, which are the log-law constants. This model is suitable for smooth walls.cai_sagaut_model_ii model is a smooth-wall explicit algebraic wall model defined as Model-II of Cai and Sagaut (DOI: 10.1063/5.0048563). It uses the same kappa and B parameters as spalding, and also accepts optional blending parameters p and s, which default to the paper calibration values 1.138 and 217.8. The main advantage of this model is that it is explicit.rough_log_law model requires specifying kappa and B, which are the log-law constants, and z0, which is the characteristic roughness height.most model is a version of the rough_log_law adapted for flows with temperature stratification, such as atmospheric boundary layer (ABL) flows. The model uses Monin-Obukhov stability theory (MOST) to account for the local temperature gradient. More details and required keywords are given below.The richardson model is similar to the most model, but it assesses the stability dependence based on the Richardson number instead of the Obukhov length. More details and required keywords are given below.
All wall models specify their sampling strategy with a sampling object. Its type is either gll, for sampling at an off-wall GLL node, or distance, for sampling at a physical wall-normal distance using global interpolation. The value entry specifies the GLL index or distance, respectively. It may be a scalar or an array. Current wall models require exactly one sampling point per wall node, so use the scalar form. GLL indices start at 1 and may not exceed the polynomial order plus 1. Distances must be positive. If sampling is omitted, the legacy h_index keyword remains available for GLL sampling.
The optional output_h entry in sampling controls whether Neko writes a diagnostic field containing the resolved wall-normal sampling distance. It defaults to true. The output uses the base name wall_model_h_<bc_name>, where <bc_name> is the boundary-condition name. The field is zero away from the wall boundary and contains the sampling distance at wall nodes.
To set sampling values separately for every wall node, set value to "user" and provide n_samples. It specifies the number of samples per wall node and must currently be 1. Neko then calls the corresponding user sampling routine once during setup; see user wall sampling. The user routine receives the boundary-condition name, so descriptive names are helpful when there are multiple wall-modelled boundaries.
A 3D field with the name tau will be registered in the field registry. At the boundary it will store the magnitude of the predicted stress. This can be used to post-process the predictions.
user_velocity, a Dirichlet boundary for more complex velocity profiles. This boundary condition uses a more advanced user interface. user_pressure, a boundary for specified non-uniform pressure profiles, similar in essence to user_velocity. overset_interface, a Dirichlet boundary condition that retrieves values from another simulation with an overlapping domain. All zone indices that belong to the same overset interface should normally be included in one boundary object, avoiding repeated calls to the interpolation routines.
The meshes must overlap by at least one element. Since the donor values come from another concurrent simulation, Neko must be executed in multiple-program-multiple-data (MPMD) mode. The simulations otherwise remain independent. They must use the same fixed timestep; variable timesteps are not currently supported for overset coupling.
The available inputs are:
| Name | Description | Admissible values | Default |
|---|---|---|---|
zone_indices | Face zones comprising the overset interface. | Integer array. | Required |
name | Name assigned to the boundary condition. | Non-empty string. | Generated from the first zone index. |
couple_pressure | Also interpolate and impose pressure from the overlapping simulation. | true or false. | false |
order | Order of the IEXT temporal extrapolation applied on the first interface update of each timestep. | Integer from 1 to 3. | 1 |
relaxation | Under-relaxation factor for subsequent Schwarz corrections within the same timestep. | Real in \((0,1]\). | 1.0 |
interpolation.tolerance | Tolerance used by the global point search. | Positive real. | NEKO_EPS*1e3 |
interpolation.padding | Padding used by the global point search. | Positive real. | 1e-2 |
log | Log the interface RMSE between the current receiver trace and newly interpolated donor data. | true or false. | false |
Pressure coupling is supported and uses the scalar overset implementation. Setting couple_pressure to true applies the same interpolation, IEXT order, relaxation factor, and logging controls to pressure.
In the first interface update of every physical timestep, the IEXT prediction is applied without relaxation. On subsequent Schwarz iterations at the same timestep, the new donor value \(\hat{g}^{k+1}\) is blended with the previously applied interface value \(g^k\) according to
\[ g^{k+1} = (1 - \omega)g^k + \omega\hat{g}^{k+1}, \qquad 0 < \omega \leq 1, \]
where relaxation is \(\omega\). Thus, relaxation does not filter data between physical timesteps and does not affect a case with zero Schwarz iterations. A value of 1.0 recovers the original unrelaxed behavior. Lower values can stabilize an oscillatory Schwarz iteration, but may slow an already contractive iteration. The IEXT order should generally be consistent with the time-integration scheme; higher-order extrapolation may require more Schwarz iterations for stability.
The most model is based on Monin-Obukhov similarity theory (Monin and Obukhov, 1954) and adds a correction to the rough log law according to
\begin{eqnarray*} \frac{\partial{V}}{\partial z} &=& \frac{u_*}{\kappa z}\phi_m\left(\frac{z}{L}\right), \\ \frac{\partial{\theta}}{\partial z} &=& \frac{\overline{(w'\theta')}}{u_* \kappa z}\phi_h\left(\frac{z}{L}\right), \end{eqnarray*}
where \(V\) is the horizontal wind speed (given that \(z\) is the wall-normal direction) and \(\theta\) is the potential temperature.
The formulations of the correction functions \(\phi_m\) and \(\phi_h\) are taken from Dyer 1974 for the convective regime, and from Holstlag and De Bruin 1988 for the stable regime.
The keywords for this model are:
kappa: The von Kàrmàn constant, defaults to 0.4 (as is the standard in the ABL literature).Pr: The turbulent Prandtl number, defaults to 1.0.z0: The characteristic roughness length for momentum.z0h: The characteristic roughness length for heat. If a negative value is given, the roughness length for heat is computed using the formula of Zilitinkevich 1995, with the provided value acting as the constant \(-A_0\) in the Zilitinkevich formula. Defaults to be the same as z0.type_of_temp_bc: Accepted values are the lowercase strings neumann or dirichlet. If neumann, the provided value of bottom_bc_flux_or_temp is used directly as the surface heat flux in the computation of the wall stress. If dirichlet, the value of bottom_bc_flux_or_temp is interpreted as a surface temperature, which is then used to compute a heat flux using the MOST relationship.bottom_bc_flux_or_temp: Value of the surface heat flux if type_of_temp_bc is neumann, or value of the surface temperature if type_of_temp_bc is dirichlet.scalar_field: The name of the scalar field to be used as the potential temperature in the equations.time_dependent_temp_bc: Boolean. If false the value of bottom_bc_flux_or_temp will be kept constant throughout the simulation. If true, the wall model will look for bc_value in neko_const_registry and assign that value at each time step. The value of bc_value can then be updated in the user file, for example in user_check. neumann or dirichlet value for the scalar field to compute the surface shear stress, but it does not set the boundary condition for the scalar. The same boundary condition should be set separately for the scalar (see Boundary conditions).Dyer, A. J. (1974). A review of flux-profile relationships. Boundary-Layer Meteorology, 7(3), 363–372. https://doi.org/10.1007/BF00240838
Holtslag, A. A. M., & De Bruin, H. A. R. (1988). Applied Modeling of the Nighttime Surface Energy Balance over Land. Journal of Applied Meteorology, 27(6), 689–704. https://doi.org/10.1175/1520-0450(1988)027%253C0689:AMOTNS%253E2.0.CO;2
Monin, A. S., & Obukhov, A. M. (1954). Basic laws of turbulent mixing in the surface layer of the atmosphere. Tr Akad Nauk SSSR Geofiz Inst, 24(151), 163–187.
Zilitinkevich, S. S., 1995: Non-local turbulent transport: Pollution dispersion aspects of coherent structure of convective flows. Air Pollution III, H. Power, N. Moussiopoulos, and C. A. Brebbia, Eds., Vol. 1, Air Pollution Theory and Simulation, Computational Mechanics Publications, 53–60.
This Richardson-number based wall model is conceptually similar to the more well-known MOST-based wall model, but it computes the effect of the temperature stratification based on the bulk Richardson number instead of the Obukhov length.
In the convective regime, the surface shear stress, \(\tau\), and surface heat flux, \(\overline{u'\theta'}\) are computed using the formulations of Louis 1979:
\begin{eqnarray*} \tau &=& a^2 V^2 F_m\left(\frac{z}{z_0}, \mathrm{Ri}_b\right), \\ \overline{u'\theta'} &=& \frac{a^2}{R}\, V\, \Delta\theta \, F_h\left(\frac{z}{z_{0h}}, \mathrm{Ri}_b\right). \end{eqnarray*}
Here, \(V\) is the horizontal wind speed (given that \(z\) is the wall-normal direction); \(\theta\) is the potential temperature; \(\mathrm{Ri}_b\) is the bulk Richardson number; \(z_0\) and \(z_{0h}\) are the roughness lengths for momentum and heat, respectively; \(F_m\) and \(F_h\) are stability functions as defined in Louis 1979; and \(a\) and \(R\) are constants, also as defined in Louis 1979.
In the stable regime, the surface shear stress and surface heat flux are computed based on Mauritsen et al. 2007:
\begin{eqnarray*} \tau &=& \frac{V^2}{\left[\ln\left(\dfrac{z}{z_0}\right)\right]^2} \,\frac{f_{\tau}(\mathrm{Ri}_b)}{f_{\tau}(0)} \left(\frac{\ell}{z}\right)^2, \\ \overline{u'\theta'} &=& \frac{\Delta\theta}{\ln\left(\dfrac{z}{z_{0h}}\right)} \,\frac{f_{\theta}(\mathrm{Ri}_b)}{\left|f_{\theta}(0)\right|} \left(\frac{\ell}{z}\right) \frac{u_*}{\mathrm{Pr}}. \end{eqnarray*}
Here, \(V, \theta, \mathrm{Ri}_b, z_0\), and \(z_{0h}\) are the same as above; \(f_{\tau}\) and \(f_{\theta}\) are defined in Mauritsen et al. 2007; \(l\) is a lengthscale (we use \(l = \kappa z\), where \(\kappa=0.4\) is the von Kàrmàn constant); \(u_*\) is the friction velocity, and \(\mathrm{Pr}\) is the turbulent Prandtl number.
The keywords for this wall model are the same as for the MOST model, and a time-varying temperature boundary condition can be applied in the same way as described for the MOST model.
Louis, J.-F. (1979). A parametric model of vertical eddy fluxes in the atmosphere. Boundary-Layer Meteorology, 17(2), 187–202. https://doi.org/10.1007/BF00117978.
Mauritsen, T., Svensson, G., Zilitinkevich, S. S., Esau, I., Enger, L., & Grisogono, B. (2007). A Total Turbulent Energy Closure Model for Neutrally and Stably Stratified Atmospheric Boundary Layers. Journal of the Atmospheric Sciences, 64(11), 4113–4126. https://doi.org/10.1175/2007JAS2294.1.
Instead of a joint type, a boundary condition object can hold a velocity object and a pressure object, each with its own type. Both must be present in this form. This allows combinations that are not covered by the joint types, for example a non_normal velocity condition with a pressure given by an expression.
The natural conditions of the scheme are selected with the neumann type. For the velocity it is a Neumann condition on the viscous traction, which is homogeneous, i.e. no constraint, unless a constant flux is given. For the pressure it is the Neumann condition of the Pn/Pn scheme, with the normal gradient given by the momentum equation, as on velocity Dirichlet boundaries. There the flux is set by the scheme, so the pressure neumann type takes no keywords. Keywords at the top level of the object, such as zone_indices and name, apply to both conditions. Keywords inside velocity and pressure apply to that condition only and take precedence over the top-level ones. A joint type cannot be combined with a velocity or pressure object in the same entry.
The combinations defined by the joint types are the ones that are tested. When combining conditions separately, it is up to the user to make sure that the combination is well posed. In particular, the pressure is only de-meaned when no zone has a Dirichlet pressure condition.
The available velocity condition types are listed below. The keywords they accept are the same as for the corresponding joint types.
| Velocity condition | Description |
|---|---|
| dirichlet | A Dirichlet condition with a constant value. Same as velocity_value. |
| expression | A Dirichlet condition given by expressions. Same as expression_velocity. |
| no_slip | A no-slip wall, see no_slip. |
| symmetry | A symmetry plane, see symmetry. |
| non_normal | A homogeneous Neumann condition for the surface-normal component combined with a Dirichlet for the surface-parallel components, see normal_outflow. |
| neumann | A Neumann condition on the viscous traction. Without keywords it is homogeneous, i.e. no constraint. The optional flux, an array of 3 reals, prescribes a constant traction per unit area, added to the right-hand side of the momentum equation. |
| shear_stress | Prescribed wall shear stress, see shear_stress. |
| wall_model | Shear stress from a wall model, see wall_model. |
| blasius_profile | A Blasius velocity profile, see blasius_profile. |
| user_dirichlet | The user-defined Dirichlet condition. Same as user_velocity. |
| overset_interface | Values from another Neko simulation, see overset_interface. |
The available pressure condition types are the following.
| Pressure condition | Description |
|---|---|
| dirichlet | A Dirichlet condition with a constant value. |
| zero_dirichlet | A zero pressure. The pressure part of outflow. |
| expression | A Dirichlet condition given by an expression. Same as expression_pressure. |
| dong | The Dong outflow condition. The pressure part of outflow+dong. |
| neumann | The Neumann condition of the Pn/Pn scheme, with the normal gradient given by the momentum equation. Takes no keywords. |
| user_dirichlet | The user-defined Dirichlet condition. Same as user_pressure. |
| overset_interface | Values from another Neko simulation, see overset_interface. |
Custom velocity and pressure condition types can be added from user modules with the register_fluid_pnpn_velocity_bc and register_fluid_pnpn_pressure_bc routines. See fluid_pnpn_bc_template.f90 in the examples/programming/user_type_templates directory for a template.
The object initial_condition is used to provide initial conditions. It is mandatory. Note that this currently pertains to both the fluid, but also scalars. The means of prescribing the values are controlled via the type keyword:
user, the values are set inside the compiled user file. as explained in the user defined initial condition section of the user file documentation.uniform, the value is a constant vector, looked up under the value keyword.expression, the value keyword holds an array of three strings, one mathematical expression per velocity component. See expressions for the syntax. For example uniform they can describe an arbitrary profile without a user file.blasius, a Blasius profile is prescribed. The boundary cannot be tilted with respect to the coordinate axes. It requires the following parameters:delta, the thickness of the boundary layer.freestream_velocity, the velocity value in the free stream.approximation, the numerical approximation to the Blasius profile.linear, linear approximation.quadratic, quadratic approximation.cubic, cubic approximation.quartic, quartic approximation.sin, sine function approximation.tanh, hyperbolic tangent approximation of Savaş (2012). In this case delta is the 99% thickness.point_zone, the values are set to a constant base value, supplied under the base_value keyword, and then assigned a zone value inside a point zone. The point zone is specified by the zone_name keyword, and should be defined in the case.point_zones object. See more about point zones.field, where the initial condition is retrieved from a field file. The following keywords can be used:
| Name | Description | Admissible values | Default value |
|---|---|---|---|
file_name | Name of the field file to use (e.g. myfield0.f00034). | Strings ending with f***** | - |
interpolate | Whether to interpolate the velocity and pressure fields from the field file onto the current mesh. | true or false | false |
mesh_file_name | If interpolation is enabled, the name of the field file that contains the mesh coordinates. | Strings ending with f***** | file_name |
interpolation.tolerance | Tolerance for the point search. | Positive real. | NEKO_EPS*1e3 |
interpolation.padding | Padding for the point search. | Positive real. | 1e-2 |
fld files that were written in double precision. To check if your fld file was written in double precision, run the command: #std 4 ... indicates single precision, whereas #std 8 ... indicates double precision. Neko writes single precision fld files by default. To write your files in double precision, set case.output_precision to "double"."interpolate" is set to true, even if the field file matches with the current simulation.The source_terms object should be used to specify the source terms in the momentum equation. The object is not mandatory, by default no forcing term is present. Each source term, is itself a JSON object, so source_terms is just an array of them. Note that with respect to the governing equations, the source terms define \( f^u \), meaning that the values are then multiplied by the density.
For each source, the type keyword defines the kind of forcing that will be introduced. Furthermore, the start_time and end_time keywords can be used to set a time frame for when the source term is active. Note, however, that these keywords have no effect on the user-defined source terms, but their execution can, of course, be directly controlled in the user code. By default, all source terms are active during the entire simulation.
The following types are currently implemented.
constant, constant forcing. Strength defined by the values array with 3 reals corresponding to the 3 components of the forcing.boussinesq, a source term introducing buoyancy based on the Boussinesq approximation, \( \rho \beta (T - T_{ref}) \cdot \mathbf{g} \). Here, \(
\rho \) is density, \( \beta \) the thermal expansion coefficient, \(
\mathbf{g} \) the gravity vector, and \( T_{ref} \) a reference value of the scalar, typically temperature.
Reads the following entries:
scalar_field: The name of the scalar that drives the source term, defaults to "s".reference_value: The reference value of the scalar.g: The gravity vector.beta: The thermal expansion coefficient, defaults to the inverse of ref_value.coriolis, a source term introducing a Coriolis force, defined as \( -2
\Omega \times (u - U_g) \). Here, \( \Omega \) is the rotation vector and \( u \) is the velocity vector, and \( U_g \) is the geostrophic wind. Several ways of setting \( \Omega \) are provided via the following keywords.
rotation_vector: Array with 3 values. Directly assigns \( \Omega \) to the provided vector.omega and phi: Both scalars. The latitude phi should be provided in degrees. Sets \( \Omega = [0, \omega \cos \phi, \omega \sin \phi ] \). Common notation when modelling the atmosphere. This assumes that the \( z
\) axis is normal to the ground.f: Scalar, referred to as the Coriolis parameter, \( f = 2 \omega \sin
\phi \). Sets \( \Omega = [0, 0, 0.5f ] \). This assumes both that \( z
\) axis is normal to the ground and that the ground-normal component of the Coriolis force is negligible.The geostrophic wind is set to 0 for all components by default. Other values are set via the geostrophic_wind keyword.
centrifugal, a source term introducing a centrifugal force, defined as \( -
\Omega \times (\Omega \times r) \). Here, \( \Omega \) is the rotation vector and \( r \) is the position relative to the reference point, which is any point lying on the rotation axis. To define forcing one has to provide \( \Omega \) and the reference point. This is provided via the following keywords.rotation_vector: Array with 3 values. Directly assigns \( \Omega \) to the provided vector.reference_point: Array with 3 values. Deifines any point on the rotaion axis.coriolis and centrifugal source terms in a consistent way.user, the values are set inside the compiled user file, using the corresponding user file subroutine.brinkman, Brinkman permeability forcing inside a pre-defined region.gradient_jump_penalty, perform gradient_jump_penalisation.sponge, adds a sponge term based on a reference velocity field, which is applied in a user-specified region of the domain.field, uses fields in the neko_registry as values of the source term. The fields are selected with the field_names keyword.hpfrt, adds a high-pass filter relaxation term to the momentum equation. This source term damps the highest elementwise Legendre modes of the velocity field. It is configured with filter_weight, the damping strength, and filter_modes, the number of highest modes affected by the filter. See High-pass filter relaxation source term for the definition of the filter and the source term.The same source-term object can be used for scalars. In this case it is added to the scalar source_terms array and acts on that scalar field.
translation Adds a forcing term corresponding to a domain translating with constant velocity. The term comes from the ALE formulation and is written as \( (\mathbf{w}\cdot\nabla)\mathbf{u} \), or, for spatially constant mesh velocity \( \nabla\cdot\mathbf{w} = 0 \), equivalently in conservative form as \( \nabla\cdot(\mathbf{u}\otimes\mathbf{w}) \). The keyword to be provided in the case file is domain_velocity, which expects an array with 3 values.
Note that in this case, we still solve for the absolute velocity, not the relative one. It is simply advection that is affected. Useful to perform simulations on domains that are moving on a periodic direction.
The Brinkman source term introduces regions of resistance in the fluid domain. The volume force \( f_i \) applied in the selected regions are proportional to the fluid velocity component \( u_i \).
\begin{eqnarray*} f_i(x) &=& - B(x) u_i(x), \\ B(x) &=& \kappa_0 + (\kappa_1 - \kappa_0) \xi(x) \frac{q + 1}{q + \xi(x)}, \end{eqnarray*}
where, \( x \) is the current location in the domain, \( \xi: x \mapsto [0,1] \) represent an indicator function for the resistance where \( \xi(x) = 0 \) is a free flow. \( \kappa_i \) describes the limits for the force application at \( \xi(x)=0 \) and \( \xi(x)=1 \). A penalty parameter \( q \) help us to reduce numerical problems.
The indicator function will be defined based on the object type. The following types are currently implemented.
boundary_mesh, the indicator function for a boundary mesh is computed in two steps. First, the signed distance function is computed for the boundary mesh. Then, the indicator function is computed using the distance transform function specified in the case file. This is currently not very well optimized, it will scale by O(log(M)*N), where M is the number of triangles in the boundary mesh and N is the number of grid points in the simulation. To avoid recomputing the distance field for multiple simulations with the same boundary mesh and numerical discretization, the distance field can be cached to a file. This is controlled by the cache keyword. If set to true, the distance field will be saved to a file specified by the cache_file keyword. If the file already exists, it will be loaded instead of recomputing it. The distance field is stored in the Nek5000 .fld file format.point_zone, the indicator function is defined as 1 inside the point zone and 0 outside.field, the indicator function is directly taken from a field in the neko_registry. The field name should be specified in the name keyword.file, the indicator function is directly taken from a field file. The file name should be specified in the file_name keyword. This requires the user to specify the field name in the file, which should be done by setting the field_name keyword.Each object are added to a common indicator field by means of a point-wise max operator. This means that the indicator field will be the union of all the regions defined by the objects.
To assist correct placement and scaling of objects from external sources, the meshes can be transformed using the mesh_transform object. The object can be used to apply a transformation to the boundary mesh. The following types are currently implemented.
none, no transformation is applied.bounding_box, the boundary mesh is transformed to fit inside a box defined by box_min and box_max. The box is defined by two vectors of 3 reals each. The keep_aspect_ratio keyword can be used to keep the aspect ratio of the boundary mesh.After the indicator field is computed, it is filtered using a filter type specified in the case file. The filter is used to smooth the indicator field before computing the Brinkman force. The following types are currently implemented.
none, no filtering is applied.The filtering can be defined for each object separately. Additionally, the filter can be specified for the entire source term, in which case it will be applied to the final indicator field, after all sources have been added.
Additional keywords are available to modify the Brinkman force term.
| Name | Description | Admissible values | Default value |
|---|---|---|---|
brinkman.limits | Brinkman factor at free-flow ( \( \kappa_0 \)) and solid domain ( \( \kappa_1 \)). | Vector of 2 reals. | - |
brinkman.penalty | Penalty parameter \( q \) when estimating Brinkman factor. | Real | \( 1.0 \) |
objects | Array of JSON objects, defining the objects to be immersed. | Each object must specify a type | - |
distance_transform.type | How to map from distance field to indicator field. | step, smooth_step | - |
distance_transform.value | Values used to define the distance transform, such as cut-off distance for the step function. | Real | - |
filter.type | Type of filtering applied to the indicator field either globally or for the current object. | none | none |
mesh_transform.type | Apply a transformation to the boundary mesh. | bounding_box, none | none |
mesh_transform.box_min | Lower left front corner of the box to fit inside. | Vector of 3 reals | - |
mesh_transform.box_max | Upper right back corner of the box to fit inside. | Vector of 3 reals | - |
mesh_transform.keep_aspect_ratio | Keep the aspect ratio of the boundary mesh. | true or false | true |
Example of a Brinkman source term where a boundary mesh and a point zone are combined to define the resistance in the fluid domain. The indicator field for the boundary mesh is computed using a step function with a cut-off distance of \( 0.1 \). The indicator field for the point zone is not filtered.
The optional gradient_jump_penalty object can be used to perform gradient jump penalty as an continuous interior penalty option. The penalty term is performed on the weak form equation of quantity \( T \) (could either be velocity or scalar) as a right hand side term
\( - < \tau |u \cdot n| h^2_{\Omega ^e} G(T) \phi_{t1} \phi_{t2} \frac{\partial \phi_{n}}{\partial n}>\),
where \( <> \) refers to the integral over all facets of the element, \( \tau \) is the penalty parameter, \( |u \cdot n| \) is the absolute velocity flux over the facet, \( h^2_{\Omega ^e} \) is the mesh size, \( G(T) \) is the gradient jump over the facet, \( \phi_{t1} \phi_{t2} \) are the polynomial on the tangential direction of the facet, and finally \( \frac{\partial \phi_{n}}{\partial n} \) is the gradient of the normal polynomial on the facet.
Here in our Neko context where hexahedral mesh is adopted, \( h^2_{\Omega ^e} \) is measured by the average distance from the vertices of the facet to the facet on the opposite side. And the distance of a vertex to another facet is defined by the average distance from the vertex to the plane constituted by 3 vertices from the other facet.
The penalty parameter \( \tau \) could be expressed as the form \( \tau = a * (P + 1) ^ {-b}\), for \( P > 1 \) where \( P \) is the polynomial order while \( a \) and \( b \) are user-defined parameters. The configuration uses the following parameters:
tau, the penalty parameter that can be only used for \( P = 1 \), default to be 0.02.scaling_factor, the scaling parameter \( a \) for \( P > 1 \), default to be 0.8.scaling_exponent, the scaling parameter \( b \) for \( P > 1 \), default to be 4.0.The sponge source term adds a term to each of the momentum equations of the form
\( \mathbf{\lambda} f(\mathbf{x}) ( \mathbf{u}^{bf} - \mathbf{u}) \)
where:
Amplitudes are specified using the amplitudes keyword with an array of 3 reals. Any of those values can be set to 0 to suppress the forcing in that particular direction. For example [1.0, 1.0, 0.0] will multiply the fringe field by 1, 1, and 0 in the x, y and z directions respectively, effectively removing the forcing in the z direction.
The reference velocity field, or baseflow can be set from four methods:
constant, applies constant values according to the value keyword:
field, where the velocity fields are retrieved from an fld file. Uses the same parameters as the field initial condition. field initial condition apply here.user, where the velocity field is set according to what is defined in the user file. Useful for setting velocity fields manually. In this case, the base flow fields must be created and added to the neko_registry (see fortran code snippet below). no-op, where the velocity field is retrieved directly from the registry. Use when multiple sponge objects should use the same baseflow velocity. Finally, the fringe field must be filled by the user. This is done in the user file by adding the field to the neko_registry, by default under the name "sponge_fringe"; the name can be changed with fringe_registry_name. The same applies to the baseflow fields when "method": "user": they are looked up under the prefix "sponge_bf", i.e. sponge_bf_u, sponge_bf_v and sponge_bf_w, which can be changed with baseflow_registry_prefix. The fields must be in the registry before the first time step. Register them in initialize, not in initial_conditions, since the latter is not called when restarting from a checkpoint.
constant or field sponge registers its own baseflow and needs a distinct baseflow_registry_prefix; use no-op with the same prefix to share one instead. A no-op sponge must come after the sponge whose baseflow it reuses, and cannot reuse a user baseflow. Sponges that share a fringe_registry_name act in the same region, so give each one its own fringe to apply them in different regions.initializeIn order to visualize your baseflow and fringe field, you may set dump_fields to true. An fld file will be written to disk as spng_fields.fld(note, not in output_directory) with the fringe field stored as pressure. You may change the name of the field file by setting dump_file_name (must have the extension fld).
The parameters for the sponge source term are summarized in the table below:
| Name | Description | Admissible values | Default value |
|---|---|---|---|
amplitudes | Sponge forcing strength in each Cartesian direction | Array of 3 reals | - |
baseflow.method | Method to define the reference (baseflow) velocity | "constant", "field", "user", "no-op" | - |
baseflow.value | Velocity vector for constant baseflow | Array of 3 reals | - |
baseflow.file_name | File containing baseflow velocity field | String | - |
baseflow.mesh_file_name | Mesh file corresponding to the baseflow field | String | - |
baseflow.interpolate | Whether to interpolate field values to current mesh | Boolean | false |
baseflow.tolerance | Tolerance for interpolation convergence | Real | - |
fringe_registry_name | Name of the fringe mask field in neko_registry | String | "sponge_fringe" |
baseflow_registry_prefix | Prefix of the base flow fields in neko_registry | String | "sponge_bf" |
dump_fields | If true, dumps the fringe and baseflow fields for visualization | Boolean | false |
dump_file_name | Name of the fld file in which to dump the base flow and fringe fields | String ending with fld | spng_fields.fld |
Neko supports the simulation of moving walls through the Arbitrary Lagrangian-Eulerian (ALE) framework. The current implementation allows for an arbitrary number of individually moving or deformable walls, collectively referred to as bodies.
The "ale" block in case file is part of the "fluid" object, and has the following high-level structure:
To run an ALE simulation, the framework must be set up as follows:
"boundary_condition" block, any moving wall must be set to "no_slip" with "moving": true."ale" block, the keyword "enabled" must be set to true."ale" block is present, the "enabled" keyword is mandatory."no_slip" boundary marked as "moving": true must be registered as an ALE body under the "bodies" object."no_slip" boundary with "moving": true.In the following, the main blocks of "ale" object are explained.
To smoothly move the mesh during an ALE simulation, we use a global smooth blending function, \( \phi_{total} \). This function is found by solving the following variable stiffness Laplace equation:
\begin{eqnarray*} \nabla \cdot (h(\mathbf{x}) \nabla \phi_{total}) = 0, \end{eqnarray*}
where \( h(\mathbf{x}) \) is the spatial mesh stiffness. Due to the linearity of the Laplace operator, the total blending function can be decomposed into the sum of individual base shapes \( \phi_i \) associated with each registered body \( i \):
\begin{eqnarray*} \phi_{total}(\mathbf{x}) = \sum_{i} \phi_i(\mathbf{x}). \end{eqnarray*}
Consequently, the problem reduces to solving a separate Laplace equation for each individual body:
\begin{eqnarray*} \nabla \cdot (h(\mathbf{x}) \nabla \phi_i) = 0. \end{eqnarray*}
Currently, the base shapes \( \phi_i \) are computed only once during initialization, and these same blending functions are used throughout the entire simulation.
During the simulation, the mesh velocity \( \mathbf{w}_{mesh}(\mathbf{x}, t) \) at any given grid point \( \mathbf{x} \) is calculated by multiplying the prescribed velocity \( \mathbf{v}_i(\mathbf{x}, t) \) of each body by its local base shape value \( \phi_i(\mathbf{x}) \), and summing the contributions across all registered bodies:
\begin{eqnarray*} \mathbf{w}_{mesh}(\mathbf{x}, t) = \sum_{i} \phi_i(\mathbf{x}) \mathbf{v}_i(\mathbf{x}, t) \end{eqnarray*}
Details regarding the configuration of the mesh stiffness \( h(\mathbf{x}) \) can be found in the Mesh Stiffness section.
Within the "solver" block, the parameters of the linear solver used to solve the Laplace equation are set. This block accepts the following keywords:
| Name | Description | Admissible values | Default value |
|---|---|---|---|
type | Type of linear solver for the Laplace equation | "cg", "gmres" | "cg" |
preconditioner.type | Type of preconditioner to use | "jacobi", "hsmg", "phmg" | "jacobi" |
absolute_tolerance | Absolute tolerance for solver convergence | Positive real | 1.0e-10 |
max_iterations | Maximum number of linear solver iterations | Positive integer | 1000 |
space_size | Krylov space size of gmres before it restarts | Positive integer | 30 |
monitor | Monitor residuals in the linear solver | true or false | false |
output_base_shape | Enables output of the base shape field \( \phi \) | true or false | true |
output_stiffness | Enables output of the computed mesh stiffness field \( h(\mathbf{x}) \) | true or false | false |
import_base_shape | Whether to import \( \phi \) fields from file | true or false | false |
If the output flags are enabled, Neko will generate .fld files during the initialization phase. These files are highly useful for verifying that the mesh deformation fields and stiffness regions are configured correctly before running the simulation:
phi_<body_name>0.f00000: Generated if "output_base_shape": true. Contains the computed base shape \( \phi_i \) for a specific body.phi_total0.f00000: Generated if "output_base_shape": true and there is more than one body registered. Contains the sum of all base shapes ( \( \phi_{total} = \sum \phi_i \)).stiffness0.f00000: Generated if "output_stiffness": true. Contains the global spatial mesh stiffness field \( h(\mathbf{x}) \).absolute_tolerance is set appropriately.user_ale_base_shapes user subroutine. In this case, the internal Laplace solver is bypassed entirely, even if the custom subroutine is only used for one of the ALE bodies. It is thus up to the user to ensure the validity of the base shape. Setting "output_base_shape": true will still write your custom user shapes to .fld files, allowing you to easily visualize and debug your custom implementations. More details about implementing this user subroutine can be found here.import_base_shape is set to true, the Laplace solve will be skipped entirely, and the \( \phi \) fields are instead loaded from previously computed .fld files. When this feature is enabled, the file to be loaded for each body must be specified using the base_shape_import_file keyword within the "bodies" block.One of the available features in the Neko ALE module is the Mesh Preview. This allows users to investigate mesh quality over time without running an expensive full fluid simulation. The preview applies the exact same prescribed motion configured by the user, ensuring the visualized mesh motion perfectly mirrors what will happen during the actual simulation. Enabling this feature generates .fld files named mesh_preview0.f*****, which contain the deformed mesh and the mass matrix.
The "mesh_preview" block accepts the following keywords:
| Name | Description | Admissible values | Default value |
|---|---|---|---|
enabled | Toggles the mesh preview feature on or off | true or false | false |
start_time | Start time for the mesh preview simulation | Positive real | 0.0 |
end_time | End time for the mesh preview simulation | Positive real | - |
output_freq | Number of timesteps between each generated output file | Positive integer | - |
dt | Constant time step size used for the mesh preview | Positive real | - |
"mesh_preview" feature is strictly a pre-processing step. Once the preview completes, whether successfully or due to a failure, the Neko run will terminate and output a corresponding success or failure message. To proceed with the actual fluid simulation, the user must set "mesh_preview.enabled: false" in order for the actual simulation to run. Additionally, if the solver detects an inverted mesh element during the preview phase, it will save the exact time step at which the Jacobian becomes negative to assist with debugging."mesh_preview" feature uses the same time integration order as defined in "case.numerics.time_order". However, the start_time, end_time, constant dt, and output_freq parameters must be explicitly defined within the "mesh_preview" block itself.enabled, Neko will always save the mesh in the .fld files at every output step.The "bodies" block defines an array of objects, where each object represents an individually controlled ALE body. Each body must map to corresponding physical boundary zone using "bodies.zone_indices".
Each individual body object accepts the following general keywords and base kinematics:
| Name | Description | Admissible values | Default value |
|---|---|---|---|
name | The name identifier of the body | String | "body_<body_ID>" |
zone_indices | The physical boundary zones associated with this body | Array of positive integers | - |
oscillation | Sub-object defining the translational oscillation kinematics | JSON object | - |
rotation | Sub-object defining the rotational kinematics applied to the body | JSON object | - |
pivot | Sub-object defining the center point for rotational kinematics | JSON object | - |
stiff_geom | Sub-object defining the mesh stiffness region | JSON object | - |
base_shape_import_file | Name of the \( \phi \) field file to import | String ending with f00000 | - |
"bodies" array, not based on their "zone_indices".no_slip zone IDs are assigned to "zone_indices" of a single ALE body, the code will treat all those boundaries as a unified rigid body.base_shape_import_file keyword is only mandatory if the solver option "ale.solver.import_base_shape" is set to true. In this case, this keyword should be provided for every single registered body.The "oscillation" sub-object defines the harmonic translational motion of the body. It takes the following mandatory keywords:
| Name | Description | Admissible values | Default value |
|---|---|---|---|
oscillation.amplitude | Amplitude of translational oscillation in \( [x, y, z] \) | Array of 3 reals | - |
oscillation.frequency | Frequency of translational oscillation in \( [x, y, z] \) | Array of 3 reals | - |
"oscillation" block is included in the case file, both "amplitude" and "frequency" become mandatory.When translational oscillation is configured, the displacement \( x_i(t) \) and velocity \( v_i(t) \) of the body follow a simple harmonic motion for each active directional component \( i \in \{x, y, z\} \):
\begin{eqnarray*} x_i(t) &=& A_i \sin(2\pi f_i t), \\ v_i(t) &=& 2\pi f_i A_i \cos(2\pi f_i t), \end{eqnarray*}
where \( A_i \) is the oscillation.amplitude and \( f_i \) is the oscillation.frequency for that specific direction.
If the body undergoes rotational motion, the "rotation" sub-object can be configured. Depending on the rotation.type, different parameters become applicable:
| Name | Description | Admissible values | Default value |
|---|---|---|---|
rotation.type | The type of rotational kinematics applied | "harmonic", "ramp", "smooth_step" | - |
rotation.amplitude_deg | Rotational amplitude in degrees (only for harmonic) | Array of 3 reals | - |
rotation.frequency | Rotational frequency (only for harmonic) | Array of 3 reals | - |
rotation.ramp_omega0 | Target angular velocity (only for ramp) | Array of 3 reals | - |
rotation.ramp_t0 | Time constant for the ramp (only for ramp) | Array of 3 reals | - |
rotation.axis | Axis of rotation (only for smooth_step) | 1 (x), 2 (y), 3 (z) | 3 |
rotation.target_angle_deg | Target rotation angle in degrees (only for smooth_step) | Real | - |
rotation.step_control_times | Control times \( [t_0, t_1, t_2, t_3] \) (only for smooth_step) | Array of 4 reals | - |
"rotation" block is included in the case file, a valid "pivot" block to specify the center of rotation must be defined. The "pivot" object is explained here. Additionally, the specific parameters corresponding to the chosen rotation.type become mandatory.The angular velocity vector \( \mathbf{\omega}(t) \) and angular position vector \( \mathbf{\theta}(t) \) are computed based on the selected rotation.type as follows:
1. Harmonic Rotation ("harmonic") Applies a simple harmonic oscillation to the angular velocity and position for each axis \( i \in \{x, y, z\} \):
\begin{eqnarray*} \omega_i(t) &=& A_{rad,i} (2\pi f_i) \cos(2\pi f_i t), \\ \theta_i(t) &=& A_{rad,i} \sin(2\pi f_i t), \end{eqnarray*}
where \( A_{rad,i} \) is the amplitude_deg converted to radians, and \( f_i \) is the frequency.
2. Ramp Rotation ("ramp") Gradually ramps up the angular velocity to a target value for each axis \( i \in \{x, y, z\} \):
\begin{eqnarray*} \omega_i(t) &=& \Omega_{0,i} \left( 1 - \exp\left(-4.6 \frac{t}{t_{0,i}}\right) \right), \\ \theta_i(t) &=& \Omega_{0,i} \left[ t - \frac{t_{0,i}}{4.6} \left( 1 - \exp\left(-4.6 \frac{t}{t_{0,i}}\right) \right) \right], \end{eqnarray*}
where \( \Omega_{0,i} \) is the target angular velocity (ramp_omega0) and \( t_{0,i}\) is the time constant (ramp_t0). The factor of 4.6 ensures the velocity reaches approximately 99% of its target at \( t = t_{0,i} \). At times beyond the ramp parameter ( \( t \gt t_{0,i} \)), the body achieves a steady rotation rate \( \Omega_{0,i} \).
3. Smooth Step Rotation ("smooth_step") Applies a smooth rotation around a single specified axis using a derivative step function, \( \text{dstep}(\tau) \), to reach a target angle and return to zero. The motion along the chosen axis is divided into four phases defined by step_control_times \( [t_0, t_1, t_2, t_3] \), prescribing both the angular velocity \( \omega(t) \) and angular position \( \theta(t) \):
\begin{eqnarray*} \omega(t) &=& \frac{\theta_{rad}}{t_1 - t_0} \text{dstep}\left(\frac{t - t_0}{t_1 - t_0}\right), \\ \theta(t) &=& \theta_{rad} \, S\left(\frac{t - t_0}{t_1 - t_0}\right) \end{eqnarray*}
\begin{eqnarray*} \omega(t) &=& 0, \\ \theta(t) &=& \theta_{rad} \end{eqnarray*}
\begin{eqnarray*} \omega(t) &=& -\frac{\theta_{rad}}{t_3 - t_2} \text{dstep}\left(\frac{t - t_2}{t_3 - t_2}\right), \\ \theta(t) &=& \theta_{rad} \left[ 1 - S\left(\frac{t - t_2}{t_3 - t_2}\right) \right] \end{eqnarray*}
\begin{eqnarray*} \omega(t) &= 0, \\ \theta(t) &= 0, \end{eqnarray*}
where \( \theta_{rad} \) is thetarget_angle_deg converted to radians.The base smooth step function \( S(\tau) \) and its analytical derivative \( \text{dstep}(\tau) \) are defined for \( \tau \in (0, 1) \) as:
\begin{eqnarray*} S(\tau) &=& \frac{1}{1 + \exp(g(\tau))}, \quad \text{where} \quad g(\tau) = \frac{1}{\tau - 1} + \frac{1}{\tau}, \\ \text{dstep}(\tau) &=& -S(\tau)(1 - S(\tau))g'(\tau), \quad \text{where} \quad g'(\tau) = -\frac{1}{(\tau - 1)^2} - \frac{1}{\tau^2}. \end{eqnarray*}
For bounds where \( \tau \le 0 \), \( S(\tau) = 0 \) and \( \text{dstep}(\tau) = 0 \). For bounds where \( \tau \ge 1 \), \( S(\tau) = 1 \) and \( \text{dstep}(\tau) = 0 \).
user_ale_rigid_kinematics or user_ale_mesh_velocity subroutines in the user file (see here for more information). The example "Double_ocyl_cylinder" shows an example for this usage. More details regarding the supported ALE interfaces and subroutines can be found here.The "pivot" sub-object defines the center point around which the body rotates.
| Name | Description | Admissible values | Default value |
|---|---|---|---|
pivot.type | Type of pivot definition | "relative" | "relative" |
pivot.value | The spatial coordinates of the rotation center \( [x, y, z] \) | Array of 3 reals | - |
"pivot") moves rigidly with the body. This means that if a body undergoes both translational oscillation and rotation, "pivot.value" defines the initial position of the rotation center. Throughout the simulation, the location of the pivot point is numerically updated using the translational velocity of the body, even if a custom rigid motion is applied using a user_ale_rigid_kinematics subroutine (see here for more information).The "stiff_geom" sub-object defines a local geometric region where the mesh is kept highly rigid to maintain the original mesh quality in regions of interest, e.g., within the boundary layer.
The global mesh stiffness field \( h(\mathbf{x}) \) is constructed by taking a base stiffness of \( 1.0 \) and adding the maximum local stiffness contribution from all registered bodies:
\begin{eqnarray*} h(\mathbf{x}) = 1.0 + \max_{b \in \text{bodies}} (\text{Stiffness}_b(\mathbf{x})). \end{eqnarray*}
"stiff_geom" definition is mandatory for every registered body.| Name | Description | Admissible values | Default value |
|---|---|---|---|
stiff_geom.type | The shape of the stiffness region | "cylinder", "sphere", "cheap_dist" | - |
stiff_geom.decay_profile | How the stiffness decays away from the body | "gaussian", "tanh" | - |
stiff_geom.gain | The gain multiplier for the stiffness field | Positive real | - |
stiff_geom.cutoff_coef | Controls the steepness of the spatial decay | Positive real | 9.0 (gaussian), 3.5 (tanh) |
stiff_geom.center | Center coordinates \( [x, y, z] \) (only for cylinder, sphere) | Array of 3 reals | - |
stiff_geom.radius | Radius of the geometry (only for cylinder, sphere) | Positive real | - |
stiff_geom.stiff_dist | Distance to maintain stiffness (only for cheap_dist) | Positive real | - |
The local stiffness contribution from a body, \( \text{Stiffness}_b(\mathbf{x}) \), is evaluated based on the chosen decay_profile. These profiles depend on a raw distance \( r \) (detailed in the next section) and a characteristic decay length \( d \), which is defined as:
"radius" if "type" is "cylinder" or "sphere"."stiff_dist" if "type" is "cheap_dist".Based on the decay_profile, the stiffness is calculated as follows:
1. Gaussian Profile ("gaussian") Applies an exponential decay based on the squared normalized distance:
\begin{eqnarray*} \text{Stiffness}_b(\mathbf{x}) = \text{gain} \cdot \exp\left( -\left(\frac{r}{d}\right)^2 \cdot \text{cutoff_coef} \right). \end{eqnarray*}
2. Tanh Profile ("tanh") Applies a smooth hyperbolic tangent transition:
\begin{eqnarray*} \text{Stiffness}_b(\mathbf{x}) = \text{gain} \cdot \left[ 1 - \tanh\left( \frac{r}{d} \cdot \text{cutoff_coef} \right) \right]. \end{eqnarray*}
For a given coordinate \( \mathbf{x} = (x, y, z) \), the raw distance \( r \) used in the stiffness formulas above is calculated based on the selected stiff_geom.type and the body's stiffness center \( C = (c_x, c_y, c_z) \):
"sphere"):\begin{eqnarray*} r = \sqrt{(x - c_x)^2 + (y - c_y)^2 + (z - c_z)^2}. \end{eqnarray*}
"cylinder"): Calculates the distance to the Z-axis passing through the center \( (c_x, c_y) \).\begin{eqnarray*} r = \sqrt{(x - c_x)^2 + (y - c_y)^2}. \end{eqnarray*}
"cheap_dist"): \( r \) is assigned from a precomputed pseudo distance field based on the boundary zone_indices.radius (from the center) or stiff_dist (from the boundary), the mesh stiffness is at its highest. If the gain parameter is set large enough, the mesh within this region moves rigidly with the body, preserving its original element quality without deformation. Users are encouraged to check the ocyl_cylinder3D, ocyl_ellipse3D, and Double_ocyl_cylinder examples to get a better idea of how these parameters are configured in practice.gain (e.g., 1.0e6) is recommended if the mesh immediately surrounding the body is intended to be fully rigid (i.e., \( \phi_i \approx 1 \)). Conversely, if two moving objects are in close proximity, the gain should be kept low enough to ensure the mesh in the gap region remains soft and deformable. Visualizing the generated phi_total0.f00000 file, and checking the mesh quality using mesh_preview are highly recommended.Neko supports checkpointing and restarting ALE simulations with both .chkp and HDF5 checkpoint files. Select HDF5 by setting checkpoint_format to "hdf5"; this requires a build with HDF5 support. No other parameters are needed beyond the usual checkpoint configuration.
Restart Capabilities:
.chkp and HDF5 checkpoints support restarting from the same mesh with a different polynomial order. The mass matrix at the restart time is then used for the lagged mass matrices required by the BDF2 and BDF3 time-integration schemes. The user must decide whether the resulting initial transient error is acceptable."ale.enabled": false, cannot be used as restart_file for an ALE simulation. However, its output field can be loaded as an initial_condition for a subsequent ALE simulation. Saving that field in double precision is recommended.The mandatory velocity_solver and pressure_solver objects are used to configure the solvers for the momentum and pressure-Poisson equation. The following keywords are used, with the corresponding options.
type, solver type.cg, a conjugate gradient solver.pipecg, a pipelined conjugate gradient solver.bicgstab, a bi-conjugate gradient stabilized solver.coupled_bicgstab, a coupled bi-conjugate gradient stabilized solver for CPU backends. It can be used for velocity when viscosity varies in space.cacg, a communication-avoiding conjugate gradient solver.coupled_cg, a coupled conjugate gradient solver. Must be used for velocity when viscosity varies in space.gmres, a GMRES solver. Typically used for pressure.fused_cg, a conjugate gradient solver optimised for accelerators usingfused_coupled_cg, a coupled conjugate gradient solver optimised for accelerators using kernel fusion. Must be used for velocity when viscosity varies in space and device backened is used. using kernel fusion.preconditioner.type, preconditioner type.jacobi, a Jacobi preconditioner. Typically used for velocity.hsmg, a hybrid-Schwarz multigrid preconditioner. Typically used for pressure.phmg, a hybrid ph multigrid preconditioner. Typically used for pressure.ident, an identity matrix (no preconditioner).absolute_tolerance, tolerance criterion for convergence.max_iterations, maximum number of iterations before giving up (default 1000).space_size, Krylov space size of gmres before it restarts (default 30); not used by the other solvers. GMRES stores 2 * space_size + 2 fields, so a smaller value saves memory. A solve that converges within space_size iterations is unaffected by it, so a good choice is the iteration count seen in steady operation.projection_space_size, size of the vector space used for accelerating the solution procedure. If 0, then the projection space is not used. More important for the pressure equation.projection_hold_steps, steps for which the simulation does not use projection after starting or time step changes. E.g. if 5, then the projection space will start to update at the 6th time step and the space will be utilized at the 7th time step.projection_reorthogonalize_basis, logical flag to update and re-orthogonalize the projection basis at each time. This option works only for pressure projection.monitor, monitoring of residuals. If set to true, the residuals will be printed for each iteration.In addition to the above settings, the solvers can be configured with strict convergence criteria. This is done by setting the case.fluid.strict_convergence keyword to true. This will force the solver to converge to the specified tolerance within the specified number of iterations. If the solver does not converge, the simulation will be terminated. This can in some situations cause issues if the initial condition is far from a valid solution. Therefore a user can allow an initial stabilization phase by setting the case.fluid.allow_stabilization keyword to true. In this case, the strict convergence will be ignored untill all components of the velocity field converge within the specified tolerance. After this initial stabilization phase, strict convergence will be enforced for the rest of the simulation.
The multilevel preconditioners, hsmg and phmg, come with an additional set of parameters related to the solution of the coarse grid problem. These parameters are specified as a coarse_grid block under preconditioner.
For hsmg, the following keywords are used:
| Name | Description | Admissible values | Default value |
|---|---|---|---|
coarse_grid.solver | Type of linear solver for the coarse grid, any of the Krylov solvers or TreeAMG tamg | A solver type | cg |
coarse_grid.preconditioner | Type of the preconditioner to use (only valid for a Krylov based solver) | A preconditioner type | jacobi |
coarse_grid.iterations | Number of linear solver iterations (only valid for a Krylov based solver) | An integer | 10 |
coarse_grid.monitor | Monitor residuals in the coarse grid (only valid for a Krylov based solver) | true or false | false |
coarse_grid.levels | Number of AMG levels to construct (only valid for solver type tamg) | An integer | 3 |
coarse_grid.iterations | Number of AMG iterations (only valid for solver type tamg) | An integer | 1 |
coarse_grid.cheby_degree | Degree of the Chebyshev based AMG smoother | An integer | 4 |
For phmg, the following keywords are used:
| Name | Description | Admissible values | Default value |
|---|---|---|---|
pcoarsening_schedule | P-multigrid coarsening schedule (polynomial order, high to low) | Array of integers | [3, 1] |
smoother_iterations | Number of smoother iterations in the p-multigrid parts | An integer | 3 |
smoother_cheby_acc | Type of Chebyshev acceleration | none, jacobi or schwarz | jacobi |
coarse_grid.levels | Number of AMG levels to construct (only valid for solver type tamg) | An integer | 3 |
coarse_grid.iterations | Number of linear solver iterations for coarse grid solver | An integer | 1 |
coarse_grid.cheby_degree | Degree of the Chebyshev based AMG smoother | An integer | 4 |
The coarse levels of phmg are built once, and are otherwise left at the geometry they were built from. If the mesh changes during the simulation, they can be updated with an update block under preconditioner.
| Name | Description | Admissible values | Default value |
|---|---|---|---|
update.enabled | Update the coarse levels when the mesh changes | true or false | false |
update.eigs.enabled | Re-estimate the Chebyshev eigenvalues on an update | true or false | true |
update.eigs.frequency | Re-estimate the eigenvalues every N-th update | An integer | 20 |
update.eigs.warm_start | Start the estimation from the previous eigenvector | true or false | true |
update.eigs.warm_start_iterations | Power iterations used for a warm started estimation | An integer | 20 |
The geometry and the smoother are updated on every mesh change. The eigs keywords only control the eigenvalue re-estimation, which is the expensive part.
The smoother accelerator is only rebuilt for smoother_cheby_acc set to jacobi, where the Jacobi diagonal is recomputed from the new geometry. With schwarz the local solves are left at the initial geometry. The coarse geometry and the eigenvalues are updated either way.
The PHMG update (as for now) also assumes that the mesh connectivity does not change.
The optional flow_rate_force object can be used to force a particular flow rate through the domain. Useful for channel and pipe flows. The configuration uses the following parameters:
direction, the direction of the flow, defined as 0, 1, or 2, corresponding to x, y or z, respectively.value, the desired flow rate.use_averaged_flow, whether value specifies the domain-averaged (bulk) velocity or the volume flow rate.log, whether to print the flow-rate forcing log message each time the forcing is adjusted. Defaults to true.All the parameters are summarized in the table below. This includes all the subobjects discussed above, as well as keyword parameters that can be described concisely directly in the table.
| Name | Description | Admissible values | Default value |
|---|---|---|---|
scheme | The fluid solve type. | pnpn | - |
name | The name associated to the fluid solver. | String | fluid |
Re | The Reynolds number. | Positive real | - |
rho | The density of the fluid. | Positive real | - |
mu | The dynamic viscosity of the fluid. | Positive real | - |
nut_field | The name of the turbulent viscosity field. | String | - |
output_control | Defines the interpretation of output_value to define the frequency of writing checkpoint files. | nsamples, simulationtime, tsteps, never | - |
output_value | The frequency of sampling in terms of output_control. | Positive real or integer | - |
output_format | The file format of field data. | nek5000, adios2, or vtkhdf | nek5000 |
output_mesh_in_all_files | Indicates if the mesh should be written in every output fld file. | true or false | false |
output_filename | The output filename. | String | field |
output_subdivide | Whether to subdivide spectral elements into linear sub-cells for VTKHDF output. | true or false | false |
inflow_condition.type | Velocity inflow condition type. | user, uniform, blasius | - |
inflow_condition.value | Value of the inflow velocity. | Vector of 3 reals | - |
initial_condition.type | Initial condition type. | user, uniform, expression, blasius, point_zone, field | - |
initial_condition.value | Value of the velocity initial condition. | Vector of 3 reals, or of 3 strings if "type" = "expression" | - |
initial_condition.file_name | If "type"="field", the name of the field file to use (e.g. myfield0.f00034). | String ending with f***** | - |
initial_condition.interpolate | If "type"="field", whether to enable interpolation on the current mesh. | true or false | false |
initial_condition.mesh_file_name | If "type"="field" and interpolation is enabled, the name of the field file that contains the mesh coordinates. | Strings ending with f***** | file_name |
initial_condition.interpolation.tolerance | If "type"="field" and interpolation is enabled, the tolerance for the point search. | Positive real. | NEKO_EPS*1e3 |
initial_condition.interpolation.padding | If "type"="field" and interpolation is enabled, the padding for the point search. | Positive real. | 1e-2 |
blasius.delta | Boundary layer thickness in the Blasius profile. | Positive real | - |
blasius.freestream_velocity | Free-stream velocity in the Blasius profile. | Vector of 3 reals | - |
blasius.approximation | Numerical approximation of the Blasius profile. | linear, quadratic, cubic, quartic, sin, tanh | - |
shear_stress.value | The shear stress vector value for sh boundaries | Vector of 3 reals | [0, 0, 0] |
wall_modelling.type | The wall model type for wm boundaries. See documentation for additional config parameters. | rough_log_law, spalding | - |
source_terms | Array of JSON objects, defining additional source terms. | See list of source terms above | - |
boundary_types | Boundary types/conditions labels. | Array of strings | - |
velocity_solver.type | Linear solver for the momentum equation. | cg, pipecg, bicgstab, coupled_bicgstab, coupled_cg, cacg, gmres | - |
velocity_solver.preconditioner.type | Linear solver preconditioner for the momentum equation. | ident, hsmg, jacobi | - |
velocity_solver.absolute_tolerance | Linear solver convergence criterion for the momentum equation. | Positive real | - |
velocity_solver.max_iterations | Linear solver max iteration count for the momentum equation. | Positive integer | 1000 |
velocity_solver.space_size | Krylov space size of gmres for the momentum equation. | Positive integer | 30 |
velocity_solver.projection_space_size | Projection space size for the momentum equation. | Positive integer | 0 |
velocity_solver.projection_hold_steps | Holding steps of the projection for the momentum equation. | Positive integer | 5 |
velocity_solver.monitor | Monitor residuals in the linear solver for the momentum equation. | true or false | false |
pressure_solver.type | Linear solver for the pressure equation. | cg, pipecg, bicgstab, cacg, gmres | - |
pressure_solver.preconditioner.type | Linear solver preconditioner for the pressure equation. | ident, hsmg, jacobi | - |
pressure_solver.absolute_tolerance | Linear solver convergence criterion for the pressure equation. | Positive real | - |
pressure_solver.max_iterations | Linear solver max iteration count for the pressure equation. | Positive integer | 1000 |
pressure_solver.space_size | Krylov space size of gmres for the pressure equation. | Positive integer | 30 |
pressure_solver.projection_space_size | Projection space size for the pressure equation. | Positive integer | 0 |
pressure_solver.projection_hold_steps | Holding steps of the projection for the pressure equation. | Positive integer | 5 |
pressure_solver.projection_reorthogonalize_basis | Whether to enable pressure projection basis reorthogonalization. | true or false | false |
pressure_solver.monitor | Monitor residuals in the linear solver for the pressure equation. | true or false | false |
flow_rate_force.direction | Direction of the forced flow. | 0, 1, 2 | - |
flow_rate_force.value | Bulk velocity or volumetric flow rate. | Positive real | - |
flow_rate_force.use_averaged_flow | Whether bulk velocity or volumetric flow rate is given by the value parameter. | true or false | - |
flow_rate_force.log | Whether to print the flow-rate forcing log message during the volume-flow adjustment. | true or false | true |
freeze | Whether to fix the velocity field at initial conditions. | true or false | false |
strict_convergence | Whether to enforce strict convergence in the linear solvers. | true or false | false |
allow_stabilization | Whether to allow an initial stabilization phase before enforcing strict convergence. | true or false | false |
advection | Whether to compute the advection term. | true or false | true |
full_stress_formulation | Whether to use the full form of the visous stress tensor term. | true or false | false |
1.9.8