BoundaryLayerDynamics.jl v1.0: a modern codebase for atmospheric boundary-layer simulations

We present BoundaryLayerDynamics.jl, a new code for turbulence-resolving simulations of atmospheric boundary-layer flows as well as canonical turbulent flows in channel geometries. The code performs direct numerical simulation as well as large-eddy simulation using a hybrid (pseudo)spectral and finite-difference approach with explicit time advancement. Written in Julia, the code strives to be flexible and adaptable without sacrificing performance, and extensive automated tests aim to ensure that the implementation is and remains correct. We show that the simulation results are in agreement with published results and that the performance is on par with an existing Fortran implementation of the same methods.

1 Introduction

Since Deardorff's early studies (Deardorff, 1969, 1970a, b), numerical simulations of the three-dimensional, unsteady flow field have become an integral part of microscale atmospheric boundary-layer (ABL) research. Direct numerical simulation (DNS) provides an extremely accurate tool to study fundamental properties of turbulent flows and their scaling from low to moderate Reynolds numbers (Moin and Mahesh, 1998). Large-eddy simulation (LES) provides a numerical model for a wide range of ABL flow phenomena at realistic Reynolds numbers while relying on modest, well-supported modeling assumptions (Meneveau and Katz, 2000; Stoll et al., 2020). Together, DNS and LES constitute the backbone of the computational study of turbulent flow dynamics and have contributed many insights to our current understanding.

Many different implementations of these methods are in use for ABL research and continue to be actively developed. Projects such as PALM (Maronga et al., 2020), OpenFOAM (Chen et al., 2014), and the Weather Research and Forecasting (WRF) model (Skamarock et al., 2021) are developed in a community effort as open-source codes with broad applicability. In addition, many research groups have their own codes that have been developed and extended over decades and are passed person to person. Some studies also rely on commercial software such as Ansys Fluent, although not having access to implementation details is problematic for scientific reproducibility.

Codes differ along many dimensions, perhaps most importantly in the physical models that are implemented. The numerical methods used to compute a solution can also lead to important differences, particularly for LES, where the smallest resolved scales of motion are an integral part of the turbulence dynamics (Kravchenko and Moin, 1997). The performance characteristics of a code are also central as most turbulence-resolving simulations continue to be limited by their computational cost. Furthermore, there can be important differences in the effort and model-specific experience required for setting up simulations and making changes and additions to the source code while ensuring the correctness of the results. When developing a new code or selecting an existing model for a simulation, these different qualities have to be weighed against each other and trade-offs are inevitable.

The advent of the Julia programming language (Bezanson et al., 2017) represents a shift in the landscape of possible trade-offs between conflicting goals. Publicly launched in 2012 and stabilized with version 1.0 in 2018, Julia promises to combine the performance of Fortran, C, and C++ with the convenience of Python, MATLAB, and R. Automatic memory management, dynamic typing with type inference, multiple dispatch, and a built-in modern package manager facilitate rapid development of clear, concise code that keeps orthogonal functionality separate. At the same time, Julia's “just-ahead-of-time” compilation model allows code to run with no or minimal computational overhead.

In this paper, we present and discuss BoundaryLayerDynamics.jl (Schmid et al., 2023b), a new code for turbulence-resolving flow simulation optimized for ABL research. The code has been written to provide core functionality for DNS and LES of channel-flow configurations with a focus on making it easy to use, adapt, and extend the code without jeopardizing the correctness of the results. To achieve better trade-offs along these dimensions, the implementation relies on the Julia programming language, on automated testing, and on a modular design.

The suitability of the new code is of course not limited to simulations of ABL turbulence. In fact, none of the current functionality is specific to ABL applications. However, the choice of physical models and numerical methods is guided by the needs of such applications and future developments will similarly prioritize those use cases.

The code solves the incompressible Navier–Stokes equations relying on (pseudo)spectral and finite difference methods for discretization in the horizontal and vertical directions respectively, making use of the Message Passing Interface (MPI) for parallelization in the vertical direction and performing fully explicit time integration. The details of the numerical methods are described in Sect. 2. The implementation is validated via a number of automated tests as described in Sect. 3, where we also present a validation against DNS and LES results computed with different codes. The performance analysis presented in Sect. 4 shows that the computational cost is comparable to a Fortran implementation of the same numerical approach and that parallel performance scales favorably up to the maximum supported number of parallel processes.

BoundaryLayerDynamics.jl is open-source software and available under the MIT License through the official Julia package repository and on GitHub (https://github.com/efpl-columbia/BoundaryLayerDynamics.jl, last access: 22 November 2023), where the public repository of the package is currently hosted. The version described in this article is archived on Zenodo (Schmid et al., 2023b).

2 Governing equations and numerical methods

The choice of governing equations and numerical methods is guided by the goal of studying the turbulent flow dynamics in the atmospheric boundary layer. Other turbulent flows in engineering and in the natural environment are also considered insofar as their requirements do not conflict with those of atmospheric boundary-layer flows.

It is well-established that the Navier–Stokes equations are an extremely accurate physical model for flows with a Knudsen number Kn≪1 and that compressibility effects are minimal for flows with a Mach number $\mathit{Ma}\lesssim \mathrm{1}/\mathrm{3}$ (Panton, 2013). Since these conditions are met for atmospheric boundary-layer flows, the incompressible Navier–Stokes equations,

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial u_i}{\partial t}}+u_j\left({\displaystyle \frac{\partial u_i}{\partial x_j}}-{\displaystyle \frac{\partial u_j}{\partial x_i}}\right)={\displaystyle \frac{\mathrm{1}}{\mathit{Re}}}{\displaystyle \frac{{\partial }^{\mathrm{2}}{u}_i}{\partial x_j\partial x_j}}-{\displaystyle \frac{\partial p}{\partial x_i}}+f_i\end{array}$$

and

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{\partial u_i}{\partial x_i}}=\mathrm{0}\end{array}$$

are used as the mathematical model for this work, given here with the rotational form of the advection term (Orszag, 1971a), for which the total kinematic pressure $p=p^{\mathrm{gauge}}/\mathit{\rho }+\frac{\mathrm{1}}{\mathrm{2}}{u}_i{u}_i$ includes both the static and dynamic pressure. Quantities are non-dimensionalized with a length scale ℒ and a velocity 𝒰, producing the Reynolds number $\mathit{Re}=\mathcal{U}\mathcal{L}/\mathit{\nu }$. The Cartesian coordinates x_i denote the primary horizontal (i=1, usually streamwise), secondary horizontal (i=2, usually cross-stream), and vertical (i=3) directions, while the corresponding velocity components are given as u_i. The term f_i denotes the components of a body force, usually gravity or a constant pressure gradient, and the kinematic viscosity ν and fluid density ρ are assumed scalar constants.

While this formulation can serve as a reasonable representation of a neutrally stratified atmospheric boundary layer, there are many important ABL processes that are not included. Coriolis forces, temperature, and humidity in particular are of central importance for many applications. The current code is meant to provide the core functionality necessary for ABL flow simulations and serve as a foundation for a more comprehensive set of physical and numerical models that can be added over time.

For flows with a moderate to high Reynolds number, current computational capabilities generally do not permit resolving the full range of scales of motion. In this case, the filtered Navier–Stokes equations,

$$\begin{array}{}\text{(3)}& \begin{array}{rl}{\displaystyle \frac{\partial {\stackrel{\mathrm{̃}}{u}}_i}{\partial t}}& +{\stackrel{\mathrm{̃}}{u}}_j\left({\displaystyle \frac{\partial {\stackrel{\mathrm{̃}}{u}}_i}{\partial x_j}}-{\displaystyle \frac{\partial {\stackrel{\mathrm{̃}}{u}}_j}{\partial x_i}}\right)+{\displaystyle \frac{\partial {\mathit{\tau }}_{ij}^{\mathrm{sgs}}}{\partial x_j}}\\ & ={\displaystyle \frac{\mathrm{1}}{\mathit{Re}}}{\displaystyle \frac{{\partial }^{\mathrm{2}}{\stackrel{\mathrm{̃}}{u}}_i}{\partial x_j\partial x_j}}-{\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{p}}{\partial x_i}}+{\stackrel{\mathrm{̃}}{f}}_i\end{array}\end{array}$$

and

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{\partial {\stackrel{\mathrm{̃}}{u}}_i}{\partial x_i}}=\mathrm{0},\end{array}$$

are used as the computational model, where ${\stackrel{\mathrm{̃}}{u}}_i$ represents the spatially filtered velocity field, i.e., ${\stackrel{\mathrm{̃}}{u}}_i=\int G(\mathit{r},\mathit{x})u_i(\mathit{x}-\mathit{r})\mathrm{d}\mathit{r}$ with G defining the filtering operation. The subgrid-scale stress tensor ${\mathit{\tau }}_{ij}^{\mathrm{sgs}}={\mathit{\tau }}_{ij}^{\mathrm{R}}-\frac{\mathrm{1}}{\mathrm{3}}{\mathit{\tau }}_{kk}^{\mathrm{R}}{\mathit{\delta }}_{ij}$ represents the anisotropic component of the residual stress tensor ${\mathit{\tau }}_{ij}^{\mathrm{R}}=\stackrel{\mathrm{̃}}{u_i{u}_j}-{\stackrel{\mathrm{̃}}{u}}_i{\stackrel{\mathrm{̃}}{u}}_j$ and has to be modeled as a function of the resolved velocity ${\stackrel{\mathrm{̃}}{u}}_i$. The modified pressure $\stackrel{\mathrm{̃}}{p}={\stackrel{\mathrm{̃}}{p}}^{\mathrm{gauge}}/\mathit{\rho }+\frac{\mathrm{1}}{\mathrm{2}}{\stackrel{\mathrm{̃}}{u}}_i{\stackrel{\mathrm{̃}}{u}}_i+\frac{\mathrm{1}}{\mathrm{3}}{\mathit{\tau }}_{ii}^{\mathrm{R}}$ now includes contributions from the filtered gauge pressure ${\stackrel{\mathrm{̃}}{p}}^{\mathrm{gauge}}$, the resolved kinetic energy, and the unresolved kinetic energy. The forcing term ${\stackrel{\mathrm{̃}}{f}}_i$ is simply the spatially filtered f_i. In the following, the same notation is used to represent both the unfiltered (DNS) and filtered (LES) equations to simplify the notation.

The past decades have seen several efforts to develop a suitable model for ${\mathit{\tau }}_{ij}^{\mathrm{sgs}}$ (Smagorinsky, 1963; Schumann, 1975; Bardina et al., 1980; Germano et al., 1991; Meneveau et al., 1996; Porté-Agel et al., 2000; Bou-Zeid et al., 2005). The current implementation includes the static Smagorinsky (1963) subgrid-scale model:

$$\begin{array}{}\text{(5)}& {\mathit{\tau }}_{ij}^{\mathrm{sgs}}=-\mathrm{2}{l}_{\mathrm{S}}^{\mathrm{2}}\mathcal{S}{S}_{ij},\end{array}$$

where $S_{ij}=\mathrm{1}/\mathrm{2}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)$ is the resolved strain rate, $\mathcal{S}=\sqrt{\mathrm{2}{S}_{ij}{S}_{ij}}$ is the characteristic or total strain rate, and l_S is the Smagorinsky length scale, taken to be the product of the filter width Δ and a constant Smagorinsky coefficient C_S.

Specifying boundary conditions for turbulent flows remains a challenging research problem, since chaotic velocity fluctuations have to be prescribed in a physically accurate manner. To avoid these difficulties, turbulence-resolving simulations are often run with periodic boundary conditions. While this requires that the problem is formulated such that it can be approximated with a periodic flow field, unphysical border regions are avoided and accurate results can be obtained in the whole domain as long as the domain is large enough to accommodate all relevant scales of motion. Furthermore, the flow field can then be expressed in terms of periodic basis functions. For a horizontally periodic domain of size L₁×L₂, the velocity field can be written as

$$\begin{array}{}\text{(6)}& u_i(x_{\mathrm{1}},x_{\mathrm{2}},x_{\mathrm{3}})=\sum \sum _{\begin{array}{c}-\mathrm{\infty }<{\mathit{\kappa }}_{\mathrm{1}}<\mathrm{\infty }\\ -\mathrm{\infty }<{\mathit{\kappa }}_{\mathrm{2}}<\mathrm{\infty }\end{array}}{\widehat{u}}_i^{{\mathit{\kappa }}_{\mathrm{1}}{\mathit{\kappa }}_{\mathrm{2}}}\left(x_{\mathrm{3}}\right){\mathrm{e}}^{i{\mathit{\kappa }}_{\mathrm{1}}\mathrm{2}\mathit{\pi }{x}_{\mathrm{1}}/L_{\mathrm{1}}}{\mathrm{e}}^{i{\mathit{\kappa }}_{\mathrm{2}}\mathrm{2}\mathit{\pi }{x}_{\mathrm{2}}/L_{\mathrm{2}}}.\end{array}$$

With similar expressions for the pressure p and the forcing f_i, we can rewrite the governing equations for a single mode with wavenumber κ₁ in x₁ direction and κ₂ in x₂ direction,

$$\begin{array}{}\text{(7)}& \begin{array}{rl}{\displaystyle \frac{\partial }{\partial t}}{\widehat{u}}_i^{{\mathit{\kappa }}_{\mathrm{1}}{\mathit{\kappa }}_{\mathrm{2}}}\left(x_{\mathrm{3}}\right)& +\sum \sum _{\begin{array}{c}-\mathrm{\infty }<{\mathit{\kappa }}_{\mathrm{1}}^{\prime }<\mathrm{\infty }\\ -\mathrm{\infty }<{\mathit{\kappa }}_{\mathrm{2}}^{\prime }<\mathrm{\infty }\end{array}}{\widehat{u}}_j^{({\mathit{\kappa }}_{\mathrm{1}}-{\mathit{\kappa }}_{\mathrm{1}}^{\prime })({\mathit{\kappa }}_{\mathrm{2}}-{\mathit{\kappa }}_{\mathrm{2}}^{\prime })}\left(x_{\mathrm{3}}\right)\\ & \left({\widehat{D}}_j^{{\mathit{\kappa }}_{\mathrm{1}}^{\prime }{\mathit{\kappa }}_{\mathrm{2}}^{\prime }}{\widehat{u}}_i^{{\mathit{\kappa }}_{\mathrm{1}}^{\prime }{\mathit{\kappa }}_{\mathrm{2}}^{\prime }}\left(x_{\mathrm{3}}\right)-{\widehat{D}}_i^{{\mathit{\kappa }}_{\mathrm{1}}^{\prime }{\mathit{\kappa }}_{\mathrm{2}}^{\prime }}{\widehat{u}}_j^{{\mathit{\kappa }}_{\mathrm{1}}^{\prime }{\mathit{\kappa }}_{\mathrm{2}}^{\prime }}\left(x_{\mathrm{3}}\right)\right)\\ & +{\widehat{D}}_j^{{\mathit{\kappa }}_{\mathrm{1}}{\mathit{\kappa }}_{\mathrm{2}}}{\widehat{\mathit{\tau }}}_{ij}^{\mathrm{sgs}{\mathit{\kappa }}_{\mathrm{1}}{\mathit{\kappa }}_{\mathrm{2}}}\left(x_{\mathrm{3}}\right)\\ & ={\displaystyle \frac{\mathrm{1}}{\mathit{Re}}}{\left({\widehat{D}}_j^{{\mathit{\kappa }}_{\mathrm{1}}{\mathit{\kappa }}_{\mathrm{2}}}\right)}^{\mathrm{2}}{\widehat{u}}_i^{{\mathit{\kappa }}_{\mathrm{1}}{\mathit{\kappa }}_{\mathrm{2}}}\left(x_{\mathrm{3}}\right)\\ & -{\widehat{D}}_i^{{\mathit{\kappa }}_{\mathrm{1}}{\mathit{\kappa }}_{\mathrm{2}}}{\widehat{p}}^{{\mathit{\kappa }}_{\mathrm{1}}{\mathit{\kappa }}_{\mathrm{2}}}\left(x_{\mathrm{3}}\right)+{\widehat{f}}_i^{{\mathit{\kappa }}_{\mathrm{1}}{\mathit{\kappa }}_{\mathrm{2}}}\left(x_{\mathrm{3}}\right),\end{array}\end{array}$$

where ${\widehat{D}}_{\mathrm{1}}^{{\mathit{\kappa }}_{\mathrm{1}}{\mathit{\kappa }}_{\mathrm{2}}}=\frac{\mathrm{2}\mathit{\pi }i{\mathit{\kappa }}_{\mathrm{1}}}{L_{\mathrm{1}}}$, ${\widehat{D}}_{\mathrm{2}}^{{\mathit{\kappa }}_{\mathrm{1}}{\mathit{\kappa }}_{\mathrm{2}}}=\frac{\mathrm{2}\mathit{\pi }i{\mathit{\kappa }}_{\mathrm{2}}}{L_{\mathrm{2}}}$, and ${\widehat{D}}_{\mathrm{3}}^{{\mathit{\kappa }}_{\mathrm{1}}{\mathit{\kappa }}_{\mathrm{2}}}=\frac{\partial }{\partial x_{\mathrm{3}}}$ are the differential operators, with i denoting the imaginary unit. For direct numerical simulations the subgrid-scale term is omitted. The continuity equation becomes

$$\begin{array}{}\text{(8)}& {\widehat{D}}_i^{{\mathit{\kappa }}_{\mathrm{1}}{\mathit{\kappa }}_{\mathrm{2}}}{\widehat{u}}_i^{{\mathit{\kappa }}_{\mathrm{1}}{\mathit{\kappa }}_{\mathrm{2}}}\left(x_{\mathrm{3}}\right)=\mathrm{0}.\end{array}$$

In the vertical direction, turbulence is not homogeneous for ABL flows and periodic boundary conditions are not applicable. For a channel geometry, boundary conditions have to be specified for u_i(x₃=0) and u_i(x₃=L₃) with the constraint that

$$\begin{array}{}\text{(9)}& \underset{\mathrm{0}}{\overset{L_{\mathrm{2}}}{\int }}\underset{\mathrm{0}}{\overset{L_{\mathrm{1}}}{\int }}{u}_{\mathrm{3}}(x_{\mathrm{3}}=\mathrm{0})-u_{\mathrm{3}}(x_{\mathrm{3}}=L_{\mathrm{3}})\mathrm{d}{x}_{\mathrm{1}}\mathrm{d}{x}_{\mathrm{2}}=\mathrm{0},\end{array}$$

which can be obtained from integrating the continuity equation over the whole domain. A number of engineering flows such as smooth-wall open- and closed-channel flows can be modeled with Dirichlet and Neumann boundary conditions. The complex boundaries of atmospheric flows can require significant modeling effort, and simulations generally have to partially resolve surfaces (e.g., immersed boundary method, terrain-following coordinates) or represent their effect with a wall model for ${\mathit{\tau }}_{i\mathrm{3}}^{\mathrm{sgs}}$, usually formulated for the discretized equations (Piomelli and Balaras, 2002). The current implementation includes an algebraic equilibrium rough-wall model defined similar to Mason and Callen (1986) with

$$\begin{array}{}\text{(10)}& {\mathit{\tau }}_{i\mathrm{3}}^{\mathrm{sgs}}(x_{\mathrm{3}}=\mathrm{0})={\displaystyle \frac{-{\mathit{\kappa }}^{\mathrm{2}}\sqrt{u_{\mathrm{1}}(x_{\mathrm{3}}^{\mathrm{ref}}{)}^{\mathrm{2}}+u_{\mathrm{2}}(x_{\mathrm{3}}^{\mathrm{ref}}{)}^{\mathrm{2}}}}{\log (x_{\mathrm{3}}^{\mathrm{ref}}/z_{\mathrm{0}}{)}^{\mathrm{2}}}}{u}_i\left(x_{\mathrm{3}}^{\mathrm{ref}}\right)\end{array}$$

for $i=\mathrm{1},\mathrm{2}$ and ${\mathit{\tau }}_{\mathrm{33}}^{\mathrm{sgs}}(x_{\mathrm{3}}=\mathrm{0})=\mathrm{0}$, where z₀ is the roughness length, κ≈0.4 is the von Kármán constant, and $x_{\mathrm{3}}^{\mathrm{ref}}>z_{\mathrm{0}}$ is a reference height at which the (resolved) velocity is obtained, usually chosen as the first grid point. To improve the near-wall behavior of the subgrid-scale model, the Smagorinsky length scale is adjusted to $l_{\mathrm{S}}^{-n}={\left(C_{\mathrm{S}}\mathrm{\Delta }\right)}^{-n}+{\left(\mathit{\kappa }{x}_{\mathrm{3}}\right)}^{-n}$ as proposed by Mason and Thomson (1992), with n=2 as the default value.

For the numerical solution of Eqs. (7) and (8), we limit ourselves to N₁×N₂ wavenumbers at N₃ vertical grid points. The wavenumbers are selected symmetrically around κ_i=0, i.e., $\left|{\mathit{\kappa }}_{\mathrm{1}}\right|\le (N_{\mathrm{1}}-\mathrm{1})/\mathrm{2}$ and $\left|{\mathit{\kappa }}_{\mathrm{2}}\right|\le (N_{\mathrm{2}}-\mathrm{1})/\mathrm{2}$. This results in an odd number of wavenumbers in each direction and avoids the need for a special treatment of Nyquist frequencies. Since ${\widehat{\mathit{\varphi }}}^{-{\mathit{\kappa }}_{\mathrm{1}}-{\mathit{\kappa }}_{\mathrm{2}}}={\widehat{\mathit{\varphi }}}^{{\mathit{\kappa }}_{\mathrm{1}}{\mathit{\kappa }}_{\mathrm{2}}\ast }$ for any real-valued ϕ, we only need to explicitly solve for half the modes and can obtain the others through complex conjugation. In vertical direction, equidistant grid points are selected from the interval [0,1], which is then mapped to the domain with a function $x_{\mathrm{3}}:\left[\mathrm{0},\mathrm{1}\right]\to \left[\mathrm{0},L_{\mathrm{3}}\right]$, ζ↦x₃(ζ). This function can be used for grid stretching in the vertical direction; the choice of $x_{\mathrm{3}}:\mathit{\zeta }\mapsto L_{\mathrm{3}}\mathit{\zeta }$ defines a uniform grid. A staggered arrangement of grid points with ζ_C at the center of the N₃ segments and ζ_I at the N₃−1 interfaces between them, i.e.,

$$\begin{array}{}\text{(11)}& {\mathit{\zeta }}_C\in \left\{{\displaystyle \frac{\mathrm{1}/\mathrm{2}}{N_{\mathrm{3}}}},{\displaystyle \frac{\mathrm{3}/\mathrm{2}}{N_{\mathrm{3}}}},\mathrm{\dots },{\displaystyle \frac{N_{\mathrm{3}}-\mathrm{1}/\mathrm{2}}{N_{\mathrm{3}}}}\right\}\end{array}$$

for $u_{\mathrm{1}},u_{\mathrm{2}},p,f_{\mathrm{1}},f_{\mathrm{2}},{\mathit{\tau }}_{ii}^{\mathrm{sgs}},{\mathit{\tau }}_{\mathrm{12}}^{\mathrm{sgs}}$, along with

$$\begin{array}{}\text{(12)}& {\mathit{\zeta }}_I\in \left\{{\displaystyle \frac{\mathrm{1}}{N_{\mathrm{3}}}},{\displaystyle \frac{\mathrm{2}}{N_{\mathrm{3}}}},\mathrm{\dots },{\displaystyle \frac{N_{\mathrm{3}}-\mathrm{1}}{N_{\mathrm{3}}}}\right\}\end{array}$$

for $u_{\mathrm{3}},f_{\mathrm{3}},{\mathit{\tau }}_{\mathrm{13}}^{\mathrm{sgs}},{\mathit{\tau }}_{\mathrm{23}}^{\mathrm{sgs}}$, avoids the need to specify boundary conditions for the pressure field, prevents odd–even decoupling, and results in a smaller effective grid spacing (Ferziger et al., 2020). When running large-eddy simulations, the discretization implicitly defines the spatial filter G and the filter width Δ is taken to be $\mathrm{\Delta }=\sqrt[\mathrm{3}]{{\mathrm{\Delta }}_{\mathrm{1}}{\mathrm{\Delta }}_{\mathrm{2}}{\mathrm{\Delta }}_{\mathrm{3}}}$ (Scotti et al., 1993) with ${\mathrm{\Delta }}_{\mathrm{1}}=L_{\mathrm{1}}/N_{\mathrm{1}}$, ${\mathrm{\Delta }}_{\mathrm{2}}=L_{\mathrm{2}}/N_{\mathrm{2}}$, and ${\mathrm{\Delta }}_{\mathrm{3}}=\mathrm{1}/N_{\mathrm{3}}\mathrm{d}{x}_{\mathrm{3}}/\mathrm{d}\mathit{\zeta }$.

The horizontal derivatives ${\widehat{D}}_{\mathrm{1}}^{{\mathit{\kappa }}_{\mathrm{1}}{\mathit{\kappa }}_{\mathrm{2}}}$ and ${\widehat{D}}_{\mathrm{2}}^{{\mathit{\kappa }}_{\mathrm{1}}{\mathit{\kappa }}_{\mathrm{2}}}$ can be computed exactly. For the vertical derivative ${\widehat{D}}_{\mathrm{3}}^{{\mathit{\kappa }}_{\mathrm{1}}{\mathit{\kappa }}_{\mathrm{2}}}$, we use central second-order finite differences on the staggered ζ nodes (Moin and Verzicco, 2016) as well as the analytical derivative of x₃(ζ), i.e.,

$$\begin{array}{}\text{(13)}& \begin{array}{rl}{\left.{\widehat{D}}_{\mathrm{3}}^{{\mathit{\kappa }}_{\mathrm{1}}{\mathit{\kappa }}_{\mathrm{2}}}{\widehat{\mathit{\varphi }}}^{{\mathit{\kappa }}_{\mathrm{1}}{\mathit{\kappa }}_{\mathrm{2}}}\right|}_{\mathit{\zeta }}& ={\left.{\displaystyle \frac{\mathrm{d}\mathit{\zeta }}{\mathrm{d}{x}_{\mathrm{3}}}}\right|}_{\mathit{\zeta }}{\left.{\displaystyle \frac{\partial {\widehat{\mathit{\varphi }}}^{{\mathit{\kappa }}_{\mathrm{1}}{\mathit{\kappa }}_{\mathrm{2}}}}{\partial \mathit{\zeta }}}\right|}_{\mathit{\zeta }}\\ & ={\left.{\displaystyle \frac{\mathrm{d}\mathit{\zeta }}{\mathrm{d}{x}_{\mathrm{3}}}}\right|}_{\mathit{\zeta }}{\displaystyle \frac{{\widehat{\mathit{\varphi }}}^{{\mathit{\kappa }}_{\mathrm{1}}{\mathit{\kappa }}_{\mathrm{2}}}(\mathit{\zeta }+\mathit{\delta }\mathit{\zeta }/\mathrm{2})-{\widehat{\mathit{\varphi }}}^{{\mathit{\kappa }}_{\mathrm{1}}{\mathit{\kappa }}_{\mathrm{2}}}(\mathit{\zeta }-\mathit{\delta }\mathit{\zeta }/\mathrm{2})}{\mathit{\delta }\mathit{\zeta }}}\\ & +\mathcal{O}\left(\mathit{\delta }{\mathit{\zeta }}^{\mathrm{2}}\right)\end{array}\end{array}$$

for any field ϕ, where $\mathit{\delta }\mathit{\zeta }\equiv \mathrm{1}/N_{\mathrm{3}}$ is the grid spacing in the ζ coordinate. Vertical derivatives are therefore evaluated at the opposite set of grid points to the ones where ϕ is defined, as typical for staggered grids. At the boundary, one-sided second-order stencils are employed. This approximation of the vertical derivatives results in a truncation error of order 𝒪(δζ²).

The non-linear advection term of Eq. (7) requires further approximations. First, some of the terms (e.g., $i=\mathrm{1},j=\mathrm{3}$) are evaluated at the opposite set of vertical grid points than where they are required and have to be interpolated. With the simple interpolation $\widehat{\mathit{\varphi }}\left(\mathit{\zeta }\right)=\mathrm{1}/\mathrm{2}(\widehat{\mathit{\varphi }}\left(\mathit{\zeta }-\mathit{\delta }\mathit{\zeta }/\mathrm{2}\right)+\widehat{\mathit{\varphi }}\left(\mathit{\zeta }+\mathit{\delta }\mathit{\zeta }/\mathrm{2}\right))+\mathcal{O}\left(\mathit{\delta }{\mathit{\zeta }}^{\mathrm{2}}\right)$, the truncation error generally increases but remains of the order 𝒪(δζ²). Furthermore, the double sum can only be computed over the resolved range of wavenumbers, i.e., $\left|{\mathit{\kappa }}_{\mathrm{1}}^{\prime }\right|\le (N_{\mathrm{1}}-\mathrm{1})/\mathrm{2}$ and $\left|{\mathit{\kappa }}_{\mathrm{2}}^{\prime }\right|\le (N_{\mathrm{2}}-\mathrm{1})/\mathrm{2}$, producing another truncation error that decreases exponentially with the number of resolved wavenumbers. The same applies to the non-linear expressions involved in the evaluation of ${\mathit{\tau }}_{ij}^{\mathrm{sgs}}$. This discretization of the advection term in rotational form conserves kinetic energy in the absence of time-integration errors (Mansour et al., 1979) as long as the grid is uniform.

To simplify the computation of non-linear terms and avoid evaluating expensive convolutions, those terms are computed on $N_{\mathrm{1}}^{\mathrm{PD}}\times N_{\mathrm{2}}^{\mathrm{PD}}$ equidistant grid points in the physical domain, relying on the fast Fourier transform (FFT) algorithm for forward and backward transforms (Orszag, 1969, 1971b). In principle, $N_{\mathrm{1}}^{\mathrm{PD}}$ and $N_{\mathrm{2}}^{\mathrm{PD}}$ are parameters that can be chosen independently of N₁ and N₂, but the choice of $N_i^{\mathrm{PD}}\ge \mathrm{1}+\mathrm{3}{\mathit{\kappa }}_i^{\max }$ avoids introducing aliasing errors for a simple product of two variables such as the resolved advection term (Patterson and Orszag, 1971). In this case, the physical-domain evaluation is equivalent to a true spectral Galerkin method computing the convolution of Eq. (7) over all resolved wavenumbers. Contributions from wavenumbers $\left|{\mathit{\kappa }}_i\right|>(N_i-\mathrm{1})/\mathrm{2}$ are discarded upon return to the Fourier domain, and vertical derivatives can be computed before or after the horizontal Fourier transforms as the two operations commute.

For the more complex non-linear expressions introduced by the SGS model, full dealiasing is generally not feasible and physical-domain evaluations incur aliasing errors in addition to the truncation errors. This approach, dubbed the pseudospectral method by Orszag (1971b), can achieve similar accuracy to a Galerkin method (Orszag, 1972). While it is common to set $N_i^{\mathrm{PD}}=N_i$ for pseudospectral approximations and only discard the Nyquist wavenumber, $N_i^{\mathrm{PD}}$ can in principle be chosen freely for more control over aliasing errors. Furthermore, it can be beneficial to choose different values of $N_i^{\mathrm{PD}}$ for each non-linear term, since computing the resolved advection requires only 9 transforms and is known to be sensitive to aliasing errors (Kravchenko and Moin, 1997; Margairaz et al., 2018), while the evaluation of the SGS model requires 15 transforms.

Combining the velocity components into a single vector $\widehat{\mathit{u}}$ of length $N_{\mathrm{1}}\times N_{\mathrm{2}}\times (\mathrm{3}{N}_{\mathrm{3}}-\mathrm{1})$, the spatially discretized momentum equation can be written as

$$\begin{array}{}\text{(14)}& {\displaystyle \frac{\mathrm{d}\widehat{\mathit{u}}}{\mathrm{d}t}}=\mathrm{Adv}\left(\widehat{\mathit{u}}\right)+{\displaystyle \frac{\mathrm{1}}{\mathit{Re}}}\mathbf{\Delta }\widehat{\mathit{u}}+{\displaystyle \frac{\mathrm{1}}{\mathit{Re}}}{\widehat{\mathit{b}}}_{\mathrm{\Delta }}-\mathbf{G}\widehat{\mathit{p}}+\widehat{\mathit{f}}.\end{array}$$

Here, Adv(⋅) is the non-linear advection operator that includes both resolved and subgrid-scale contributions. Δ is the linear operator that computes the Laplacian of each velocity component, with ${\widehat{\mathit{b}}}_{\mathrm{\Delta }}$ denoting the contributions from vertical boundary conditions. G is the linear operator that computes the gradient of the pressure, discretized as a vector $\widehat{\mathit{p}}$ of size $N_{\mathrm{1}}\times N_{\mathrm{2}}\times N_{\mathrm{3}}$. The forcing term $\widehat{\mathit{f}}$, a vector of the same length as $\widehat{\mathit{u}}$, can be constant in time and space (e.g., pressure-driven channel flow), vary in time only (e.g., constant-mass-flux channel flow), vary in space only (e.g., baroclinic flow), or even vary in time and space as a function of $\widehat{\mathit{u}}$, which can be used to model vegetation drag or to represent complex geometry with an immersed boundary method. The discrete continuity equation,

$$\begin{array}{}\text{(15)}& \mathbf{D}\widehat{\mathit{u}}+{\widehat{\mathit{b}}}_{\mathrm{D}}=\mathbf{0},\end{array}$$

contains the linear divergence operator D with the contributions ${\widehat{\mathit{b}}}_{\mathrm{D}}$ from the vertical boundary conditions of u₃.

This hybrid approach, relying on spectral approximations in horizontal direction (pseudospectral for the evaluation of ${\mathit{\tau }}_{ij}^{\mathrm{sgs}}$) and second-order-accurate finite differences in vertical direction, has long been employed for computational studies of turbulent flows in channel geometries (Moin and Kim, 1982; Moeng, 1984; Albertson and Parlange, 1999a, b). It combines the fast convergence and low dissipation of spectral methods (Giacomini and Giometto, 2021) with the ease of parallelization and simple handling of boundary conditions of finite differences. Conversely, handling complex domains and non-periodic boundaries can be problematic, though still possible (Chester et al., 2007; Schmid, 2015; Li et al., 2016).

Following Perot (1993), we obtain the expressions for time integration through a block LU decomposition of the fully discretized equations. This results in expressions in the style of the fractional step method (Chorin, 1968; Temam, 1969), but it avoids the need for boundary conditions for the intermediate velocity and the pressure on a staggered grid and can easily be adapted when new terms are included or different numerical methods are employed.

Adams–Bashforth methods solving ordinary differential equations of the form $\mathrm{d}\mathit{u}/\mathrm{d}t=f\left(\mathit{u}\right),\mathit{u}\left(t_{\mathrm{0}}\right)={\mathit{u}}_{\mathrm{0}}$ can be written as ${\mathit{u}}^{(n+\mathrm{1})}={\mathit{u}}^{\left(n\right)}+\mathrm{\Delta }t{\sum }_{i=\mathrm{0}}^{s-\mathrm{1}}{\mathit{\beta }}_{i}f\left({\mathit{u}}^{(n-i)}\right)$, where β_i are the coefficients of the method, s is the order of accuracy, and superscripts denote the time step (Hairer et al., 1993). With β₀=1 this corresponds to the forward Euler method (s=1), while ${\mathit{\beta }}_{\mathrm{0}}=\mathrm{3}/\mathrm{2},{\mathit{\beta }}_{\mathrm{1}}=-\mathrm{1}/\mathrm{2}$ gives second-order accuracy (s=2). Applied to the momentum Eq. (14), this can be written as

$$\begin{array}{}\text{(16)}& {\widehat{\mathit{u}}}^{(n+\mathrm{1})}={\widehat{\mathit{u}}}^{\left(n\right)}+\mathrm{\Delta }t\sum _{i=\mathrm{0}}^{s-\mathrm{1}}{\mathit{\beta }}_i\left(\mathrm{F}\left({\widehat{\mathit{u}}}^{(n-i)}\right)-\mathbf{G}{\widehat{\mathit{p}}}^{(n-i)}\right),\end{array}$$

where terms are grouped with the definition $\mathrm{F}\left(\widehat{\mathit{u}}\right)\equiv \mathrm{Adv}\left(\widehat{\mathit{u}}\right)+\frac{\mathrm{1}}{\mathit{Re}}\mathbf{\Delta }\widehat{\mathit{u}}+\frac{\mathrm{1}}{\mathit{Re}}{\widehat{\mathit{b}}}_{\mathrm{\Delta }}+\widehat{\mathit{f}}$ to simplify the notation. Together with the continuity Eq. (15), the fully discretized equations become

$$\begin{array}{}\text{(17)}& {\widehat{\mathit{u}}}^{(n+\mathrm{1})}+\mathbf{G}{\widehat{\mathit{\phi }}}^{(n+\mathrm{1})}={\widehat{\mathit{u}}}^{\left(n\right)}+\mathrm{\Delta }t\sum _{i=\mathrm{0}}^{s-\mathrm{1}}{\mathit{\beta }}_i\mathrm{F}\left({\widehat{\mathit{u}}}^{(n-i)}\right)\end{array}$$

and

$$\begin{array}{}\text{(18)}& \mathbf{D}{\widehat{\mathit{u}}}^{(n+\mathrm{1})}=-{\widehat{\mathit{b}}}_{\mathrm{D}}\end{array}$$

if we group the pressure contributions with ${\widehat{\mathit{\phi }}}^{(n+\mathrm{1})}\equiv \mathrm{\Delta }t{\sum }_{i=\mathrm{0}}^{s-\mathrm{1}}{\mathit{\beta }}_i{\widehat{\mathit{p}}}^{(n-i)}$. Similarly, explicit s-stage Runge–Kutta methods can be written as ${\mathit{u}}^{(n,i)}={\sum }_{k=\mathrm{0}}^{i-\mathrm{1}}\left({\mathit{\alpha }}_{ik}{\mathit{u}}^{(n,k)}+\mathrm{\Delta }t{\mathit{\beta }}_{ik}f\left({\mathit{u}}^{(n,k)}\right)\right)$ for $i=\mathrm{1},\mathrm{\dots },s$, with ${\mathit{u}}^{(n,\mathrm{0})}={\mathit{u}}^{\left(n\right)}$ and ${\mathit{u}}^{(n+\mathrm{1})}={\mathit{u}}^{(n,s)}$ (Gottlieb et al., 2009). This is referred to as the Shu–Osher form (Shu and Osher, 1988), where the coefficients α_ik and β_ik are not uniquely determined by the Butcher tableau of the method and can be chosen to minimize storage requirements. In this case, the fully discretized equations can be written as

$$\begin{array}{}\text{(19)}& {\widehat{\mathit{u}}}^{(n,i)}+\mathbf{G}{\widehat{\mathit{\phi }}}^{(n,i)}=\sum _{k=\mathrm{0}}^{i-\mathrm{1}}\left({\mathit{\alpha }}_{ik}{\widehat{\mathit{u}}}^{(n,k)}+\mathrm{\Delta }t{\mathit{\beta }}_{ik}\mathrm{F}\left({\widehat{\mathit{u}}}^{(n,k)},t^{(n,k)}\right)\right)\end{array}$$

and

$$\begin{array}{}\text{(20)}& \mathbf{D}{\widehat{\mathit{u}}}^{(n,i)}=-{\widehat{\mathit{b}}}_{\mathrm{D}}\end{array}$$

with the definition ${\widehat{\mathit{\phi }}}^{(n,i)}\equiv \mathrm{\Delta }t{\sum }_{k=\mathrm{0}}^{i-\mathrm{1}}{\mathit{\beta }}_{ik}{\widehat{\mathit{p}}}^{(n,k)}$.

Each step or stage requires the solution of a system of equations in the form $\widehat{\mathit{u}}+\mathbf{G}\widehat{\mathit{\phi }}=\widehat{\mathit{a}}$ and $\mathbf{D}\widehat{\mathit{u}}=\widehat{\mathit{b}}$. This can be solved with the LU decomposition:

$$\begin{array}{}\text{(21)}& \begin{array}{rl}\left(\begin{array}{cc}\mathbf{I}& \mathbf{G}\\ \mathbf{D}& \mathbf{0}\end{array}\right)\left(\begin{array}{c}\widehat{\mathit{u}}\\ \widehat{\mathit{\phi }}\end{array}\right)& =\left(\begin{array}{cc}\mathbf{I}& \mathbf{0}\\ \mathbf{D}& -\mathbf{DG}\end{array}\right)\left(\begin{array}{cc}\mathbf{I}& \mathbf{G}\\ \mathbf{0}& \mathbf{I}\end{array}\right)\left(\begin{array}{c}\widehat{\mathit{u}}\\ \widehat{\mathit{\phi }}\end{array}\right)\\ & =\left(\begin{array}{cc}\mathbf{I}& \mathbf{0}\\ \mathbf{D}& -\mathbf{DG}\end{array}\right)\left(\begin{array}{c}{\widehat{\mathit{u}}}^{\star }\\ \widehat{\mathit{\phi }}\end{array}\right)\\ & =\left(\begin{array}{c}\widehat{\mathit{a}}\\ \widehat{\mathit{b}}\end{array}\right),\end{array}\end{array}$$

where ${\widehat{\mathit{u}}}^{\star }\equiv \widehat{\mathit{u}}+\mathbf{G}\widehat{\mathit{\phi }}$ has been introduced. The steps to compute the solution are therefore

$$\begin{array}{}\text{(22)}& \begin{array}{rl}{\widehat{\mathit{u}}}^{\star }& =\widehat{\mathit{a}},\\ \mathbf{DG}\widehat{\mathit{\phi }}& =\mathbf{D}{\widehat{\mathit{u}}}^{\star }-\widehat{\mathit{b}},\text{ and }\\ \widehat{\mathit{u}}& ={\widehat{\mathit{u}}}^{\star }-\mathbf{G}\widehat{\mathit{\phi }},\end{array}\end{array}$$

where the second step requires solving a linear system. Since the operator (DG) has no coupling between different wavenumbers, Eq. (22) can be decomposed into N₁×N₂ tridiagonal systems of size N₃. For ${\mathit{\kappa }}_{\mathrm{1}}={\mathit{\kappa }}_{\mathrm{2}}=\mathrm{0}$ the system is singular due to the fact that the governing equations only include the gradient of the pressure and do not place any restrictions on the absolute magnitude of the pressure variable. This system therefore has to be solved iteratively or the equations have to be regularized, e.g., by specifying an arbitrary value for one element of $\widehat{\mathit{\phi }}$. The current implementation relies on the Thomas algorithm (Quarteroni et al., 2007) to solve the tridiagonal systems and includes the forward Euler and second-order Adams–Bashforth methods for time integration as well as the strong stability preserving Runge–Kutta methods SSPRK (2,2) and SSPRK (3,3) (Gottlieb et al., 2009). Adding other explicit methods is straightforward, provided they can be formulated in a similar fashion.

The simulation code is written in the Julia programming language, relying on the Julia bindings to the FFTW library (Frigo and Johnson, 2005) for fast Fourier transforms. For parallelization, the domain is vertically split into up to N₃ blocks that are computed by separate processes exchanging information through the Message Passing Interface (MPI).

Figure 1 Error convergence for transient two-dimensional laminar flows. The left panel shows second-order convergence as the vertical grid resolution is refined for flows set up along the vertical direction and a randomly chosen horizontal direction. The other panels show first-, second-, and third-order convergence as the time resolution is refined for a Taylor–Green vortex set up along the two horizontal directions, in which case the spatial discretization is exact and the order of convergence of the time integration methods is measured. Grid lines show the formal order of convergence for each case.

3 Model validation

The validation efforts presented in this section aim to confirm that the numerical methods are implemented faithfully and that these methods produce physically relevant results. To maintain this confidence as the code is inevitably modified, a focus is placed on automated tests that can be rerun after every change. A set of automated unit tests verifies the expected order of accuracy when computing individual terms of the discretized equations for prescribed velocity fields and when applying the time-integration algorithms to ordinary differential equations. A set of automated integration tests verifies that the solution to canonical transient two-dimensional laminar flows can be simulated with the expected order of accuracy. Finally, fully turbulent flow solutions are computed and compared to published results produced with different codes. These tests are not automated since they require significant computational resources and have no analytical solution to compare against so there is some degree of judgment required to evaluate the quality of the solution.

The automated tests of individual terms make use of the fact that the implemented numerical methods are exact for certain velocity fields. The diffusion term and the pressure solver are exact for a function that is the product of truncated Fourier series along horizontal dimensions and a quadratic polynomial in vertical direction. The advection term is only exact for a linear function in vertical direction due to linear interpolations, although the term is still second-order accurate. By constructing such a function with randomized parameters, each term can be computed numerically as well as analytically, and matching values give a high degree of confidence in the correctness of the implementation. Furthermore, we can verify the order of convergence when computing the terms at different grid resolutions for a velocity field that cannot be handled exactly by the implemented methods. The time integration methods are verified in a similar way by solving ordinary differential equations that have analytical solutions with different time steps. We also verify that the tridiagonal solver is exact for random inputs.

To test the full solver including time integration, the automated tests include a number of laminar flow problems, currently the transient Poiseuille and Couette flows as well as decaying Taylor–Green vortices. Numerical solutions computed at different resolutions are then compared to the analytical solution to ensure that the order of convergence corresponds to formal order of the numerical methods, as shown in Fig. 1. Dimensional parameters such as domain sizes and velocity scales are again chosen randomly since parameters that are zero or unity can mask errors in the solution. For Poiseuille and Couette flows, this includes the horizontal direction of the flow. The two-dimensional Taylor–Green vortices are oriented both in horizontal and vertical planes. For the former, the spatial discretization is exact so the test case verifies the order of convergence of the time integration method.

The above test cases have been verified in single-process (serial) mode as well as in multi-process (parallel) mode. Multi-process tests are run both with a vertical resolution greater than and equal to the number of processes since those configurations sometimes rely on different code paths. Since the tests are automated and run within minutes on consumer hardware, they can be rerun whenever changes are made to the code to ensure that any future version of the code still satisfies all the tested properties.

Turbulence-resolving flow simulations require substantially more computation and are therefore not included in the automated tests that are meant to be run routinely during code development. The validation cases presented below are chosen such that they represent scientifically relevant flow systems while keeping the computational cost moderate. Each case can be simulated in about 2 h using 32 MPI processes on a single compute node of the Intel Skylake generation.

The results of direct numerical simulations are not supposed to depend on the exact method used for modeling the flow, at least for lower-order flow statistics. For relatively low Reynolds numbers, simulations have been run with many different codes and with a wide range of parameters such as domain sizes, aspect ratios, and grid resolutions, so the expected simulation results are well-established and have been validated against wind-tunnel measurements (Kim et al., 1987; del Álamo and Jiménez, 2003; Lee and Moser, 2015).

Figure 2 Direct numerical simulation (DNS) of a turbulent channel flow at Re_τ≈180, validated against data published by Lee and Moser (2015). Mean profiles are shown for the streamwise velocity $u_{\mathrm{1}}^{+}$, the advective transport $u_{\mathrm{1}}^{+}{u}_{\mathrm{3}}^{+}$, the diffusive transport $\partial u_{\mathrm{1}}^{+}/\partial x_{\mathrm{3}}^{+}$, as well as the production 𝒫⁺ and (pseudo)dissipation ε⁺ of turbulent kinetic energy. The last panel shows contours of the premultiplied turbulent kinetic energy spectra $E_{ii}^{+}$ along the streamwise ($k_{\mathrm{1}}^{+}=\mathrm{2}\mathit{\pi }{\mathit{\kappa }}_{\mathrm{1}}/L_{\mathrm{1}}^{+}$) and cross-stream ($k_{\mathrm{2}}^{+}=\mathrm{2}\mathit{\pi }{\mathit{\kappa }}_{\mathrm{2}}/L_{\mathrm{2}}^{+}$) direction. The superscript ⁺ marks values in inner units, i.e., non-dimensionalized with the friction velocity u_τ and the kinematic viscosity ν.

In Fig. 2, we show a comparison of a closed-channel flow at Re_τ≈180 with data published by Lee and Moser (2015). The friction Reynolds number ${\mathit{Re}}_{\mathit{\tau }}=u_{\mathit{\tau }}\mathit{\delta }/\mathit{\nu }$ is based on the half-channel height δ and the friction velocity $u_{\mathit{\tau }}^{\mathrm{2}}=.\mathit{\nu }\frac{\partial u_{\mathrm{1}}}{\partial x_{\mathrm{3}}}{|}_{x_{\mathrm{3}}=\mathrm{0}}$ here. The simulation is run with a bulk Reynolds number of ${\mathit{Re}}_{\mathrm{b}}=U_{\mathrm{b}}\mathit{\delta }/\mathit{\nu }=\mathrm{20}\mathrm{000}/\mathrm{7}$, where the vertically averaged bulk velocity U_b is held constant by a dynamically adjusted pressure forcing. The solution is computed in a domain of size 4πδ in streamwise and 2πδ in cross-stream direction. The velocity field is discretized with 255×191 Fourier modes at 96 vertical grid points that are spaced according to a sinusoidal grid transform $x_{\mathrm{3}}\left(\mathit{\zeta }\right)=\mathit{\delta }+\mathit{\delta }\sin \left((\mathrm{2}\mathit{\zeta }-\mathrm{1})\mathit{\eta }\mathit{\pi }/\mathrm{2}\right)/\sin (\mathit{\eta }\mathit{\pi }/\mathrm{2})$ with η=0.97. The mean statistics computed over ∼17.5 large-eddy turnover times $T_{\mathit{\tau }}=\mathit{\delta }/u_{\mathit{\tau }}$ after a spinup time of ∼3.5 T_τ closely match the results from Lee and Moser (2015).

For large-eddy simulation, validation is not as straightforward since results remain relatively sensitive to differences in the modeling approach and in the grid resolution. To validate the new implementation, we limit ourselves to a comparison with a pre-existing Fortran implementation of the same physical and numerical models (Giometto et al., 2017) and refer to previous publications for validation studies and discussions of limitations of the modeling approach (Porté-Agel et al., 2000; Yue et al., 2007, 2008; Giometto et al., 2016).

Figure 3 Large-eddy simulation (LES) of a turbulent channel flow at Re_τ=10⁸ with an aerodynamically rough wall and a channel height of $h/z_{\mathrm{0}}={\mathrm{10}}^{\mathrm{4}}$, validated against a tried and tested Fortran code with the same numerical approach (Giometto et al., 2017). All values are non-dimensionalized with the friction velocity u_τ and the roughness length z₀. Mean profiles are shown for the streamwise velocity u₁, the resolved transport u₁u₃, the subgrid-scale transport ${\mathit{\tau }}_{\mathrm{13}}^{\mathrm{sgs}}$, as well as the production 𝒫 and (pseudo-)dissipation ε of resolved turbulent kinetic energy. The last panel shows contours of the resolved turbulent kinetic energy spectra E_ii along the streamwise direction, premultiplied with the wavenumber $k_{\mathrm{1}}=\mathrm{2}\mathit{\pi }{\mathit{\kappa }}_{\mathrm{1}}/L_{\mathrm{1}}$.

In Fig. 3 we show the results of this comparison for an open-channel flow at Re_τ=10⁸ driven by a constant body force f₁. The friction Reynolds number ${\mathrm{Re}}_{\mathit{\tau }}=u_{\mathit{\tau }}h/\mathit{\nu }$ is based on the channel height h and the friction velocity $u_{\mathit{\tau }}^{\mathrm{2}}=h{f}_{\mathrm{1}}$ here. The lower surface is characterized by a roughness length z₀ that results in a non-dimensional channel height of $h/z_{\mathrm{0}}={\mathrm{10}}^{\mathrm{4}}$. The solution is computed in a domain of size 2πh in streamwise and $(\mathrm{4}/\mathrm{3})\mathit{\pi }h$ in cross-stream direction. The velocity field is discretized with 63×63 Fourier modes at 64 equidistant vertical grid points. The mean statistics computed over 200 large-eddy turnover times $T_{\mathit{\tau }}=h/u_{\mathit{\tau }}$ after a spinup time of 50 T_τ closely match for the two separate implementations.

Combined, these validation efforts provide ample evidence that the implementation matches the mathematical formulation of the methods and that those methods are capable of accurately simulating flow physics, within the limitations of the physical models. A comprehensive set of easily repeatable validation tests serves both to verify the current implementation and to ensure that future developments do not jeopardize correctness. This should not only facilitate adding new functionality but also help make changes to existing functionality and avoid getting locked into design decisions that might prove suboptimal for future developments.

4 Performance and scaling

Defined in a broad way, performance can be understood as the time required to obtain a solution at the required quality given the available computational resources. We can examine how the time changes as a function of the required quality and the available resources (relative performance) or how fast different methods arrive at a solution for fixed quality and resources (absolute performance). However, it is difficult to measure the overall quality of a turbulent flow simulation in a quantitative way since the system is chaotic and analytic solutions are not available. Furthermore, there is a great diversity of computational resources that vary along important dimensions such as floating point operations per second (FLOPS), memory bandwidth and latency, network bandwidth and latency, and many more. We therefore narrow the scope of the performance analysis to the question of the time required to obtain a solution given the specified simulation parameters and the number of compute nodes, as measured on a fairly typical high-performance computing (HPC) system.

For the implemented explicit time integration schemes, the computational cost of a single evaluation of the right-hand side of $\mathrm{d}\mathit{u}/\mathrm{d}t=f\left(\mathit{u}\right)$ fully characterizes the overall cost of a simulation, which is a simple function of the number of steps and the evaluations per step (i.e., stages of a Runge–Kutta method). The computational cost of a single evaluation of f(u) depends primarily on the number of Fourier modes and vertical grid points and whether a subgrid-scale term is modeled (LES) or not (DNS). The impact of other parameters such as the type of pressure forcing (constant flux vs. constant force), boundary conditions, and grid transformations is imperceptible.

Figure 4 Performance and scaling of individual terms and full time step, as measured on Intel Xeon Ice Lake nodes of the Stampede2 system at the Texas Advanced Computing Center. Individual terms are computed at a resolution of $\mathrm{256}\times \mathrm{256}\times \mathit{\lambda }{N}_{\mathrm{p}}$, where $\mathit{\lambda }\in \mathit{\{}\mathrm{1},\mathrm{2},\mathrm{4},\mathrm{8},\mathrm{16}\mathit{\}}$ is the number of vertical grid points per MPI process and N_p is the total number of MPI processes. Dotted lines indicate weak scaling, dashed lines indicate strong scaling, and the grid lines correspond to perfect scaling. The last panel shows the overall performance for the resolution $\mathrm{256}\times \mathrm{256}\times \mathrm{1280}$ (highlighted in orange on other panels) in comparison to a pre-existing Fortran implementation of the same numerical approach (Giometto et al., 2017). DNS performance is shown with the symbol +, LES performance with the symbol ×.

To assess relative performance, Fig. 4 shows the time required to compute the advection, diffusion, and pressure terms for different numbers of compute nodes and vertical grid points. The figure displays both strong scaling, where the number of processes is varied for a problem of fixed size, and weak scaling, where the problem size is varied in proportion to the computational resources. These results shows that the advection term contributes most to the overall cost, while the pressure term exhibits the most problematic scaling behavior.

Computing the advection term is a global operation in horizontal direction but only involves neighboring nodes in vertical direction. The bulk of the computational work consists of computing discrete Fourier transforms, which are local to each MPI process and scale as 𝒪(Nlog N), where N is the number of modes. It appears that this cost dominates over the cost of communication, resulting in near-perfect strong and weak scaling. When computing subgrid-scale stresses with a static Smagorinsky model, additional transforms are required and the cost increases to almost twice as much without affecting the scaling behavior.

Computing the pressure term has no data dependency between horizontal modes but is a sequential, global process in vertical direction (Thomas algorithm). Horizontal modes can be processed in batches to stagger the sequential passes up and down the domain, where the size of those batches is a tuning parameter that represents a trade-off between maximizing parallelism and minimizing per-batch overhead. The resulting performance shows imperfect weak and strong scaling. Scaling appears to improve when there are more vertical grid points per process, increasing the work-to-communication ratio. While the overall cost appears to remain at most about a quarter of the cost of the advection term, it is possible that the two costs are even closer for some combinations of hardware configurations and simulation parameters, in which case it could be worth optimizing the batch size parameter of the pressure solver.

Computing the diffusion term only involves neighboring vertical grid points and has no global data dependencies. This results in near-perfect weak scaling. Strong scaling is not quite perfect, which is explained by the fact that there is very little work to do for each grid point so the work-to-communication ratio is low. This has no discernible effect on the overall scaling behavior however, as computing the diffusion term is always at least an order of magnitude less work than computing the advection term.

To assess absolute performance, Fig. 4 includes a comparison with a Fortran code that implements the same numerical methods (Giometto et al., 2017). While such a comparison does not answer the question of whether either code is making optimal use of the computational resources, it does respond to the practical question of whether there are any performance trade-offs when substituting the new code for a codebase that has been actively used for turbulence research for over 2 decades. The comparison shows that the overall performance of both implementations is of a similar order of magnitude, with the new Julia implementation showing somewhat better scaling and significantly faster DNS performance.

It appears that the new Julia code has avoided introducing excessive overhead without much effort devoted to performance optimization. That it even surpasses the performance of the Fortran implementation is likely explained by two factors. First, the new code is formulated with the Fourier domain representation at its center, which makes it easier to avoid unnecessary Fourier transforms than in the physical-space formulation of the Fortran implementation. Second, the new code makes different trade-offs between work and communication which appear to be more suitable for modern hardware. Some of these insights will flow back to the Fortran implementation, reducing the performance discrepancy between the two codes.

Overall, the performance characteristics of the new code are as expected. The computational cost is dominated by the Fourier transforms necessary to compute the non-linear term, while the pressure solver shows the least favorable scaling properties, and the overall performance is comparable to a Fortran implementation of an equivalent numerical scheme. The analysis shows that for the implemented numerical scheme and current HPC hardware, performance is optimized by reducing the number and size of Fourier transforms and choosing an efficient implementation of the fast Fourier transform algorithm. Other details matter less as long as the computational cost can be kept significantly below the cost of the non-linear term.

Future improvements are likely to focus on parallelism along horizontal coordinate directions, either through multi-threaded CPU code or through GPGPU computation, allowing the code to scale to larger systems. There is also room for optimization in how communication is handled, which might become important if the work-to-communication is decreased through additional parallelism. However, performance optimizations always have to be weighed against their impact on code simplicity and ease of adaptation. Since the computational cost of flow simulations is a strongly non-linear function of the grid resolution, large performance differences are required for a practical difference in the scientific problems that can be tackled.

On the suitability of Julia for high-performance computation

Conceptually, there are two main differences in the performance characteristics of Julia compared to languages like Fortran, C, and C++. The first difference concerns the handling of types. For a statically typed language like Fortran, the type of every variable is specified and therefore available to the compiler, which can use the information to generate efficient machine code. In Julia, machine code is generated when a function is called for the first time, at which point the types of the function arguments are known. For subsequent function calls with the same argument types, the compiled code is reused while a new copy of the function is compiled for calls with different argument types. If the types of all variables inside a function can be derived from the type of function arguments, the compiler has access to the same information as for a statically typed language and can in principle generate the same or equivalent machine code. The second difference concerns the handling of memory allocations. For a language with (semi-)manual memory management like Fortran, memory is allocated and freed either manually with explicit commands or following deterministic rules. In Julia, memory is automatically allocated whenever required for the operation that is performed and the program is periodically interrupted to examine which of the memory is still in use and which memory can be freed (garbage collection). Vectorized Julia code written in a naive way often allocates large amounts of memory for intermediate results, but writing Julia code that operates in place and avoids such allocations is relatively easy with experience. Therefore, there is no a priori reason that code written in Julia should be significantly slower (or faster) than code written in Fortran, C, or C++, although Julia does make it much easier to accidentally include expensive operations that result in poor performance.

Julia was chosen as implementation language for this code with the goal of improving the ease of development, hoping that the negative effect on performance could be kept minimal. However, the experience so far has shown that the net impact on performance might even be positive. Performance optimizations can be seen as a continuum from low level (e.g., vectorized CPU instructions, optimal use of CPU cache) to high level (e.g., choice of algorithms, speed–accuracy trade-offs). At the lower end of this range, Julia relies on the LLVM compiler framework and a number of pre-existing libraries, making use of countless hours of optimization work, but Julia also makes it rather easy to write code it cannot optimize very well. Maintaining close-to-optimal performance therefore requires regular measurements and occasional fixes. High-level optimizations on the other hand require understanding the performance characteristics of different approaches and choosing the right one, often by measuring the performance of different implementations. This type of optimization work benefits from the Julia language features and the ease of integrating packages from a growing ecosystem. The impact of language choice on performance depends not so much on what optimizations are theoretically possible but which ones are simple enough that they are done in practice. It is therefore possible that the choice of Julia will be a benefit rather than a drawback for the performance of this code over its lifetime.

5 Conclusions

Turbulence-resolving ABL flow simulations are subject to a number of competing requirements that have to be considered when developing simulation code. Availability of physical and numerical models, performance and scalability, ease of use and ease of modification, safeguards against implementation and usage errors, as well as license terms may vary considerably, and trade-offs are often inevitable. The Julia programming language is promising more favorable trade-offs by offering the ergonomics of a modern high-level language without sacrificing performance.

In this paper, we have introduced a new code for turbulence-resolving flow simulations, designed for the requirements of atmospheric boundary-layer research and written in Julia. The performance is shown to be in line with a Fortran implementation of the same modeling approach. In fact, it even appears that easier experimentation with algorithmic approaches and implementation trade-offs might have a stronger impact on performance than the remaining computational overhead compared to highly optimized Fortran compilers.

The code also places a focus on automated testing and minimizing the chances for errors both during development and usage. This is particularly important in exploratory research, where the expected behavior of a new model or flow system is not known a priori, and it is difficult to discern between inconspicuous errors and novel results.

The code provides the core functionality for both direct numerical simulation and large-eddy simulation in channel-flow geometries. In the future, we expect to expand the scope by adding functionality such as more advanced subgrid-scale models, support for temperature, humidity and transport of passive scalars, and partially resolved complex terrain.