Towards exascale for wind energy simulations

Abstract

We examine large-eddy-simulation modeling approaches and computational performance of two open-source computational fluid dynamics codes for the simulation of atmospheric boundary layer flows that are of direct relevance to wind energy production. The first code, NekRS, is a high-order, unstructured-grid, spectral element code. The second code, AMR-Wind, is a second-order, block-structured, finite-volume code with adaptive mesh refinement capabilities. The objective of this study is to co-develop these codes in order to improve model fidelity and performance for each. These features will be critical for running ABL-based applications such as wind farm analysis on advanced computing architectures. To this end, we investigate the performance of NekRS and AMR-Wind on the Oak Ridge Leadership Facility supercomputers Summit, using 4 to 800 nodes (24 to 4,800 NVIDIA V100 GPUs), and Crusher, the testbed for the Frontier exascale system, using 18 to 384 Graphics Compute Dies on AMD MI250X GPUs. We compare strong- and weak-scaling capabilities, linear solver performance, and time to solution. We also identify leading inhibitors to parallel scaling.

Keywords

Exascale scalability large-eddy simulation

1. Introduction

Atmospheric boundary layer (ABL) flows are an important part of everyday life. Aside from being a primary driver of vertical exchanges in moisture, aerosols, and atmospheric gases, the ABL affects practical aspects of life—including the transportation system, renewable energy generation, pollution dispersion, noise propagation, and transmission of electromagnetic signals. ABL flows are turbulent, and the state of the turbulence is affected by density stratification that arises in large part from surface heating and cooling. Additionally, Coriolis effects caused by planetary rotation and curvature complicate the flow. Furthermore, regional-scale weather patterns and terrain add complexity to the ABL. Significant research effort is applied to ABL flows because of their importance and complexity (Moeng, 1984; Berg and Kelly, 2020; Beare et al., 2006; Sullivan et al., 2008; Kosović and Curry, 2000; Pedersen et al., 2014; Mirocha, 2020; Churchfield and Moriarty, 2020). This work focuses on numerical computation of ABL flows using large-eddy simulation (LES), where the governing physics equations are solved in filtered form such that the larger, energy-containing eddies are directly resolved, and the remaining “subgrid-scale” (SGS) turbulence is modeled. LES was born out of ABL research roughly five decades ago (Lilly, 1962; Smagorinsky, 1963), and continues to evolve and improve.

Wind energy is a prime example of an application driven by the ABL. Generation of electrical energy from farms of wind turbines at night in the stable ABL is a particularly interesting situation. The winds tend to be stronger, so generation is higher. With decreased turbulence, wind turbine wakes persist for longer distances, significantly affecting wind farm efficiency and fatigue loads on waked wind turbines. With this example in mind, researchers wish to increase grid resolution to reduce reliance on the SGS turbulence model, but they also wish to increase the overall domain size to encompass the wind farm, which commonly extends many kilometers horizontally. Increased domain size is desirable in many other applications besides wind energy. For example, LES can be used to study deep convection, which happens over a large geographical extent many kilometers into the atmosphere, and there is a push toward LES of regional-scale weather.

High-fidelity LES of the turbulent ABL is dependent on massively parallel high-performance computing (HPC). HPC architectures are evolving from traditional homogeneous x86-CPU-based computing. For example, the world’s second fastest computer (as of June 2022), the supercomputer Fugaku at the RIKEN Center for Computational Science, is built around Fujitsu’s custom ARM A64FX processor and does not use a GPU. Alternatively, the U.S. Department of Energy (DOE) has embraced a hybrid CPU-GPU approach for its leadership-class computing. Summit, the world’s seventh fastest computer (November 2023 TOP500 list (TOP500, 2023)), has nodes that house two IBM POWER9 CPUs (each with 22 cores) and six NVIDIA V100 GPUs and is capable of 200 × 10¹⁵ floating-point operations per second (FLOPS). Similarly, DOE’s first exascale-class supercomputers, those capable of at least 10¹⁸ FLOPS, will be hybrid CPU-GPU–based systems. Frontier, the world’s first exascale-class supercomputer at the Oak Ridge Leadership Computing facility, has nodes that house one AMD CPU and 4 AMD MI250X GPUs.

This collaborative effort had its genesis in a 2015 DOE workshop report (Sprague et al., 2017), where it was noted that the transition to exascale computing brings many opportunities in computational fluid dynamics (CFD), as well as significant challenges. Two of the grand challenge opportunities described in that workshop report are relevant to this paper: the simulation of boundary layer turbulence over large areas and the simulation of an entire wind farm under realistic atmospheric flow conditions. The transition to GPU-accelerated computing is significant for those performing CFD for numerical weather forecasting or for LES of ABL flows. While hybrid CPU-GPU processing potentially enables dramatically faster computing (at low power), legacy CPU-based codes need significant overhauls or rewrites to run effectively in a CPU-GPU environment (Robinson and Sprague, 2020; Sprague et al., 2017).

Several groups have introduced CFD codes for LES of atmospheric flows with an emphasis on GPU calculations. While early efforts in weather forecasting on GPUs focused at O(100) m resolution (see, e.g., Schalkwijk et al. (2015)), more recent efforts have performed high-fidelity GPU-based LES on O(1) m grid sizes. Van Heerwaarden et al. (2017) introduced the MicroHH 1.0 incompressible-flow solver directed at atmospheric flow; the solver is based on finite-difference spatial discretization and a split-operator time integration. The authors showed that for problem sizes that fit, a single GPU had performance similar to that of 32 CPU cores. Sauer and Muñoz-Esparza (2020) introduced the FastEddy LES model that was created for CPU and GPU systems. FastEddy solves the fully compressible Navier–Stokes (NS) equations using finite-difference spatial discretization and explicit Runge–Kutta time integration. The authors showed excellent scaling on up to 32 GPUs and argued that one GPU provides performance similar to that of 256 CPU cores. Recent high-order incompressible flow codes using fast tensor-product operator evaluation include ExaDG (Arndt et al., 2020), SPECHT_FS (Huismann et al., 2019), and deal. ii (Krank et al., 2017). ExaDG and deal. ii use a discontinuous Galerkin formulation, whereas SPECHT_FS uses a continuous Galerkin formulation similar to that in Nek5000/RS, which is described below.

In this paper we examine LES modeling approaches and computational performance of two open-source, incompressible-flow, GPU-oriented CFD codes that employ fundamentally different spatial discretization and data structures. The first is NekRS (Fischer et al., 2022), which is an updated version of the Nek5000 code. Nek5000/RS is a high-order, unstructured-grid, spectral-finite-element CFD code. The second is AMR-Wind, which is part of the ExaWind simulation suite (Sharma et al., 2024; Sprague et al., 2020). AMR-Wind is a second-order, block-structured, finite-volume-method CFD code with adaptive mesh refinement (AMR) capabilities and is built on the AMReX library, a software framework for massively parallel, block-structured applications. Both of these codes are part of the DOE Exascale Computing Project, which is supporting the development of GPU-ready applications for exascale-class supercomputers (ECP, 2017; Alexander et al., 2020).

We compare NekRS and AMR-Wind predictions and performance on the well-known GEWEX (Global Energy and Water Cycle Experiment) Atmospheric Boundary Layer Study (GABLS) stably stratified benchmark LES case (Beare et al., 2006), which is illustrated in Figure 1. The flow comprises a quiescent mean flow above $\approx 200$ m (going to the right and into the page in Figure 1) with a sheared turbulent boundary layer (going to the right and out of the page) over the 0–200 m vertical range. The flow is coupled with thermal buoyancy effects. The computational domain is 400 m × 400 m × 400 m and doubly periodic in the streamwise (x) and cross-flow (y) directions. (Boundary conditions on the upper and lower surfaces are detailed in the Simulations section.) Potential temperature distributions at Reynolds number Re = 50 × 10⁶ and t = 6 h are illustrated in Figure 2 along with profiles in Figure 3.

In addition to investigating these codes’ ability to represent ABL physics, an objective of this study is to co-develop these codes in order to improve model fidelity and performance, features that will be critical for running ABL-based applications such as wind farm analysis on advanced computing architectures. By careful cross-comparison, both codes have made significant advances. This article focuses on performance. A future article will address subgrid-scale modeling for LES of the ABL. Here we investigate the scalability of NekRS and AMR-Wind on the Oak Ridge Leadership Computing Facility supercomputer Summit, using 4 to 800 nodes (p = 24 to 4,800 NVIDIA V100 GPUs). We provide iteration counts, average time per step, and the real-time ratio (ratio of wall-clock time to physical time), as well as detailed performance metrics. We additionally include a limited set of timing data for the two codes on Crusher, using up to 384 Graphics Compute Dies (GCDs) on AMD MI250X GPUs (one MPI rank per GCD).

The paper is organized as follows. The next section describes the codes and gives an overview of the numerical approaches used. We then discuss the numerical setup of the simulations. We provide detailed studies comparing the codes’ performance and scaling. We conclude with a brief summary.

2. Description of codes

The numerical results are based on LES, which requires enhanced dissipation to provide an energy drain at the grid scale. Consequently, the incompressible NS and potential temperature equations are solved in a spatially filtered resolved-scale formulation, expressed in nondimensional form as

\begin{array}{l} \frac{\partial {\bar{u}}_{i}}{\partial t} + {\bar{u}}_{j} \frac{\partial {\bar{u}}_{i}}{\partial x_{j}} & = & - \frac{1}{\bar{ρ}} \frac{\partial \bar{p}}{\partial x_{i}} - \frac{\partial τ_{i j}}{\partial x_{j}} + f_{i} - \frac{θ^{'}}{θ_{0}} g_{i}, \end{array}

(1)

\begin{array}{l} \frac{\partial {\bar{u}}_{j}}{\partial x_{j}} & = & 0, \end{array}

(2)

\begin{array}{l} \frac{\partial \bar{θ}}{\partial t} + {\bar{u}}_{j} \frac{\partial \bar{θ}}{\partial x_{j}} & = & - \frac{\partial τ_{θ j}}{\partial x_{j}} \end{array},

(3)

where an overbar denotes the LES filtering operation such that

{\bar{u}}_{i}

is the ith component of the resolved-scale velocity vector,

\bar{ρ}

is the density,

\bar{p}

is the pressure, g_i is the gravity acceleration vector, and

\bar{θ}

is the potential temperature in the resolved scale. The scalar θ′/θ₀ that dictates the sign and strength of the buoyancy force is obtained from

\frac{θ^{'}}{θ_{0}} = \frac{\bar{θ} - θ_{0}}{θ_{0}},

(4)

where θ₀ is the reference potential temperature and f_i includes the Coriolis acceleration. In addition, τ_ij and τ_θj are the stress tensors in the momentum and energy equations, respectively, which include (and are dominated by) SGS modeling terms

τ_{i j} = - \frac{2}{R e} S_{i j} + τ_{i j}^{s g s} = - \frac{1}{R e} (\frac{\partial {\bar{u}}_{i}}{\partial x_{j}} + \frac{\partial {\bar{u}}_{j}}{\partial x_{i}}) + τ_{i j}^{s g s},

(5)

and

τ_{θ j} = - \frac{1}{P e} \frac{\partial \bar{θ}}{\partial x_{j}} + τ_{θ j}^{s g s},

(6)

where Re is the Reynolds number, Pe is the Peclet number, S_ij is the resolved-scale strain-rate tensor, and

τ_{i j}^{s g s}

and

τ_{θ j}^{s g s}

are the subgrid-scale stress tensors.

The SGS modeling in Nek5000/RS is based on the work of Sullivan et al. (1994), where the SGS stress tensors $τ_{i j}^{s g s}$ and $τ_{θ j}^{s g s}$ are expressed in terms of a nonisotropic, mean-field eddy viscosity (MFEV) obtained by the horizontally averaged mean strain rate, and an isotropic, fluctuating part. (Horizontal averaging refers to averaging along both horizontal directions, x and y, and is performed on each timestep.) The SGS model of Sullivan et al. (1994) is based on the following expression:

τ_{i j}^{s g s} = - 2 ν_{t} γ S_{i j} - 2 ν_{T} 〈 S_{i j} 〉,

(7)

where the angle brackets

〈 〉

denote averaging over the homogeneous directions and ν_T is an average eddy viscosity that is expressed in terms of mean flow quantities. In equation (7) γ is an “isotropy factor,” which accounts for variability in the SGS constants due to anisotropy of the mean flow. Sullivan et al. (1994) obtain the fluctuating eddy viscosity, ν_t, by using an eddy viscosity model based on the SGS turbulent kinetic energy equation, in which the shear production term is computed from the fluctuating velocities as suggested by Schumann (1975a).

Here the fluctuating (isotropic) part is taken into account through the use of either a high-pass filter (HPF) (Stolz et al., 2005) or a Smagorinsky (SMG) model based on the fluctuating strain rate. For the former model, which is not eddy-viscosity based, ν_t in equation (7) is by definition equal to zero. On the other hand, Sullivan et al. (1994) derive the expression for ν_T so that the law-of-the-wall behavior can be recovered in the absence of any resolved turbulence.

The SGS modeling in AMR-Wind is based on the work of Smagorinsky (1963), where a single partial differential equation for subgrid-scale kinetic energy is solved, and from that a subgrid-scale eddy viscosity is computed. The Boussinesq eddy viscosity hypothesis is then invoked to obtain the subgrid-scale stress tensor and heat flux vector.

In the following subsections we discuss the details of the numerical approaches of Nek5000/RS and AMR-Wind. For simplicity, we use u_i, p, and θ, dropping the overbar notation from ${\bar{u}}_{i}$ , $\bar{p}$ , and $\bar{θ}$ in the remaining sections.

2.1. Nek5000/RS

Nek5000 is a spectral element code that is used for a wide range of thermal-fluids applications. It employs high-order spectral elements (Patera, 1984) in which the solution, data, and test functions are represented as locally structured Nth-order tensor-product polynomials on a set of E globally unstructured curvilinear hexahedral brick elements. The approach yields two principal benefits. First, for smooth functions such as solutions to the incompressible NS equations, high-order polynomial expansions exhibit rapid convergence with approximation order, often yielding a significant reduction in the number of unknowns (n ≈ EN³) required to reach engineering tolerances. Second, the locally structured forms permit local lexicographical ordering with minimal indirect addressing and, crucially, the use of tensor-product sum factorization to yield low O(n) storage costs and O(nN) work complexities (Orszag, 1980).

NekRS (Fischer et al., 2022) is a GPU-accelerated version of Nek5000 that targets high performance on forthcoming exascale platforms. For performance portability, NekRS is written in C++/OCCA (Medina et al., 2014). Several key kernels are based on highly tuned OCCA kernels coming from the development work of Warburton and co-workers in the libParanumal library (Chalmers et al., 2020). Specific attention in NekRS has been given to ensure scalability to p = 10⁴–10⁵ ranks and beyond (Fischer et al., 2021). NekRS retains access to the standard Nek5000 interface, which allows users to leverage existing user-specific source code such as statistical analysis tools for turbulence.

Time integration in Nek5000/RS is based on a semi-implicit splitting scheme using kth-order backward differences (BDFk) to approximate the time derivative coupled with implicit treatment of the viscous and pressure terms and kth-order extrapolation (EXTk) for the remaining advection and forcing terms. This approach leads to independent elliptic subproblems comprising a Poisson equation for the pressure, a coupled system of Helmholtz equations for the three velocity components, and an additional Helmholtz equation for the potential temperature. The pressure Poisson equation is obtained by taking the divergence of the momentum equation and forcing $\partial u_{i}^{n} / \partial x_{i} = 0$ at time tⁿ = nΔt. The velocity and temperature Helmholtz equations are obtained once pⁿ is known:

\begin{array}{l} - \frac{\partial^{2}}{\partial x_{j} \partial x_{j}} p^{n} & = & q^{n}, \end{array}

(8)

\begin{array}{l} \frac{β_{0}}{Δ t} u_{i}^{n} - \frac{\partial}{\partial x_{j}} (\frac{1}{R e} + γ ν_{t}) 2 S_{i j}^{n} & = & - \frac{\partial}{\partial x_{i}} p^{n} + r_{i}^{n}, \end{array}

(9)

\begin{array}{l} [\frac{β_{0}}{Δ t} - \frac{\partial}{\partial x_{j}} (\frac{1}{P e} + γ ν_{t}) \frac{\partial}{\partial x_{j}}] θ^{n} & = & s^{n} \end{array},

(10)

where β₀ is an order-unity constant associated with BDFk (Fischer et al., 2017); S_ij is the resolved-scale strain-rate tensor as described in (5); and qⁿ,

r_{i}^{n}

, and sⁿ represent the sum of the values from the previous timesteps for the contributions from BDFk and EXTk. Also included in

r_{i}^{n}

and sⁿ are eddy diffusivity terms coming from the mean-field eddy diffusivity mentioned above. The fully coupled system of Helmholtz equations for the three velocity components (9) is used only when the fluctuating (isotropic) part of the SGS stress tensor is modeled using an SMG model based on the fluctuating strain rate. When this part is modeled through the use of a high-pass filter (HPF) (Stolz et al., 2005), the resulting Helmholtz equations for the three velocity components are not coupled.

With the given time-splitting, we recast (8)–(10) into weak form and derive the spatial discretization by restricting the trial and test spaces to be in the finite-dimensional space spanned by the spectral-element basis. The discretization leads to a sequence of symmetric positive definite linear systems for pressure, velocity, and temperature. Velocity and temperature are diagonally dominant and readily addressed with Jacobi-precondition conjugate gradient iteration.¹

The pressure Poisson solve is treated with GMRES using p-multigrid as a preconditioner. Details of the formulation can be found in Fischer and Lottes (2004); Fischer et al. (2017); Phillips et al. (2022).

NekRS supports several features to accelerate performance, including overlapped communication and computation during operator evaluation, which yields a 10%–15% performance gain for NS simulations; FP32 local-operator inversion and residual evaluation for the Chebyshev-accelerated Schwarz-based p-multigrid; and projection of the velocity and pressure solutions onto the space of prior solutions to generate an initial guess, which can yield a 1.5–2-fold NS performance gain (Fischer, 1998). On the NVIDIA A100, the OCCA-based kernels are close to the bandwidth-limited roofline and sustain 2.1–2.2 TFLOPS (FP64) for the Poisson operator, and 3.1–3.8 TFLOPS (FP64) for the advection operator. In the pressure preconditioner, the forward Poisson operator on the coarser multigrid levels realizes 2.5–3.9 TFLOPS (FP32), and the Schwarz smoother sustains 2.5–5.1 TFLOPS (FP32). (The lower values are for smaller values of N that are used in the p-multigrid V-cycle.) Comparable values are realized on the NVIDIA V100s on Summit, save that they are ≈1.5 times lower than those on the A100. Sustained flop rates for the full NS solver are $\approx 470$ GFLOPS per V100 on Summit.

2.2. AMR-wind

AMR-Wind is a spatially and temporally second-order accurate finite-volume code. Important aspects of the discretization are discussed below; for more details readers can see the work of Almgren et al. (1998) since the discretization is similar to the incompressible-flow solver IAMR. Velocity, scalar quantities, and gradients of pressure are located at cell centers, whereas pressure is located at nodes. Partial staggering combined with an approximate projection method yields linear systems that are well studied, have small-bandwidth stencils, and can be efficiently solved with standard techniques such as geometric multigrid. These discretization choices give a well-balanced mix of both efficiency and accuracy. In addition to the spatial staggering there is staggering in time similar to a Crank–Nicolson formulation. The time discretization is

\begin{array}{l} \frac{c_{k}^{n + 1} - c_{k}^{n}}{Δ t} & + {[\frac{\partial c_{k} u_{j}}{\partial x_{j}}]}^{n + 1 / 2} \\ = \frac{1}{ρ^{n + 1 / 2}} \frac{\partial q_{k j}^{n + 1}}{\partial x_{j}} + G_{k}^{n + 1 / 2}, \end{array}

(11)

\begin{array}{l} \frac{u_{i}^{*} - u_{i}^{n}}{Δ t} & + {[\frac{\partial u_{i} u_{j}}{\partial x_{j}}]}^{n + 1 / 2} \\ = \frac{1}{ρ^{n + 1 / 2}} (\frac{\partial τ_{i j}^{n + 1}}{\partial x_{j}} - {\frac{\partial p}{\partial x_{i}}}^{n - 1 / 2}) + F_{i}^{n + 1 / 2} \end{array},

(12)

where c indicates a scalar quantity and u denotes velocity. The index n represents a timestep, the index i runs over the three momentum equations, the index k runs over the scalar equations, and the repeated index j indicates summation.

F_{i}^{n + 1 / 2}

and

G_{k}^{n + 1 / 2}

are source terms and are evaluated at timestep n + 1/2. The term ρ is density, which in these simulations is constant in time and space, but we leave the time level to show at what point in time density is evaluated if the code is run in variable-density mode.

The advection term is formed by extrapolating in time by using a Godunov method (Almgren et al., 1998). Specifically, the velocity is first extrapolated in space (to the faces) and in time to n + 1/2 in a predictor step. MAC projection (Bell et al., 1991) is applied to ensure that the face velocities are divergence free, which takes the form

P^{MAC} (u_{i}^{f}) = u_{i}^{f} - \frac{1}{ρ^{n}} (\frac{\partial ψ}{\partial x_{i}}),

(13)

where

u_{i}^{f}

represents a face velocity and ψ is a Lagrange multiplier located on the cells. Setting the divergence equal to zero forms a variable coefficient Poisson equation:

\frac{\partial}{\partial x_{j}} (\frac{1}{ρ^{n}} \frac{\partial ψ}{\partial x_{j}}) = \frac{\partial u_{j}^{f}}{\partial x_{j}} .

(14)

The Poisson equation is discretized by using a cell-centered 7-point stencil in 3D, which is efficiently solved by using multilevel multigrid as a linear solver (Zhang et al., 2019). Once the velocity on the faces is divergence free, the advection terms are formed. Options for discretizing these advection terms include Godunov PLM (Van Leer, 1977), PPM (Colella and Woodward, 1984), and WENO-Z (Motheau and Wakefield, 2020). The Godunov schemes are high-order accurate and use an extended stencil. The scalar equations (e.g., potential temperature) are advanced one at a time by solving a Helmholtz problem. This Helmholtz problem is discretized by using a cell-centered finite-difference method forming a 7-point stencil in 3D (Almgren et al., 1998). The momentum equations are saved for last to allow the source terms to be evaluated at n + 1/2 by using the previously updated scalar equations. A good example is the Boussinesq buoyancy term, which adds a source term at timestep n + 1/2 to the momentum equation. The Boussinesq buoyancy term is a function of the already advanced potential temperature, which is averaged to time t^n+1/2 = (n + 1/2)Δt. The scalar equations and momentum equations can be solved by using geometric multigrid, or, if diagonally dominant enough, a Krylov method such as bicgstab (biconjugate gradient stabilized) is sufficient. In this work, because of the small timestep, we use only bicgstab in all of the Helmholtz solves.

The intermediate velocity $u_{i}^{*}$ is advanced by solving a Helmholtz problem in tensor form. However, this velocity vector $u_{i}^{*}$ is not guaranteed to be divergence free. An approximate projection method is used to solve for the velocity at time tⁿ⁺¹ = (n + 1)Δt:

u^{n + 1} = P (u_{i}^{*}),

(15)

where the nodal projection P is defined to be

P (u_{i}^{*}) = u_{i}^{*} + \frac{Δ t}{ρ^{n + 1 / 2}} ({\frac{\partial p}{\partial x_{i}}}^{n - 1 / 2} - \frac{\partial ϕ}{\partial x_{i}}) .

(16)

This approximate projection is different from the algorithm of Almgren et al. (1998) and more similar to the projection of Almgren et al. (2000) labeled version 2. Taking the divergence of equation (16) and setting it equal to zero, we have

\frac{\partial}{\partial x_{j}} (\frac{Δ t}{ρ^{n + 1 / 2}} \frac{\partial ϕ}{\partial x_{j}}) = \frac{\partial}{\partial x_{j}} (u_{j}^{*} + \frac{Δ t}{ρ^{n + 1 / 2}} {\frac{\partial p}{\partial x_{j}}}^{n - 1 / 2}),

(17)

where ϕ is a Lagrange multiplier related to the pressure field and solved on the nodes. To solve the nodal projection in Equation (17), a variational form is used, which leads to a 27-point stencil (Almgren et al., 1998). The right-hand side of (17) is formed by taking finite differences across cells and averaging them to each node. Each node has four contributions in each coordinate direction. Once the solution ϕ is obtained, the velocity is updated by using (16), and the pressure and its gradient are updated by using

p^{n + 1 / 2} \leftarrow ϕ, {\frac{\partial p}{\partial x_{i}}}^{n + 1 / 2} \leftarrow \frac{\partial ϕ}{\partial x_{i}} .

The gradient ∂ϕ/∂x_i is approximated by taking finite differences of ϕ along edges and averaging each edge to the cell center. For each coordinate direction, four edges contribute to the value at the cell center.

3. Simulations

We consider a stable ABL in which the ground temperature (at z = 0) is cooler than the air temperature and where the ground temperature continues to cool over the duration of the simulation. Here, we give full details of the numerical setup. The simulation domain is Ω = L_x × L_y × L_z = 400 m × 400 m × 400 m, with x the streamwise direction, y the spanwise direction, and z the vertical direction. Simulations are initialized (at t = 0) with constant velocity in the streamwise direction equal to geostrophic wind speed of U = 8 m/s. The initial potential temperature is 265 K in 0 ≤ z ≤ 100 m and linearly increased at a rate of 0.01 K/m in 100 m $\leq z \leq 400$ m. The reference potential temperature is 263.5 K. The Reynolds number is Re = UL_b/ν, where L_b = 100 m is the thickness of the initial thermal boundary layer and ν is the molecular viscosity. An initial random perturbation is added to the potential temperature field, based on the work of Beare et al. (2006) for the GABLS case, with an amplitude of 0.1 K and zero mean for 0 < z < 50 m to stimulate turbulence.

Periodic boundary conditions (BCs) are used in the streamwise and spanwise directions. At the top boundary (z = 400 m), a stress-free, rigid lid is applied for momentum, and the heat flux for the energy equation is set consistent with the 0.01 K/m temperature gradient initially prescribed in the upper region of the flow. At the bottom boundary, we perform simulations with impenetrable traction BCs for the velocity where the specified shear stress comes from Monin–Obukhov similarity theory (Monin and Obukhov, 1954). For the energy equation, a heat flux is applied that is derived from the same theory and a specified potential temperature difference between the flow at a height, z₁, and the surface. The surface temperature is from the GABLS specification following the rule θ_b(t) = 265 − 0.25 t, where t is in hours. Because the boundary conditions are periodic (lateral) or the mass flow rate through the boundaries is zero (top and bottom), pressure boundary conditions are not needed.

In Nek5000/RS, the implementation of the traction BCs for the horizontal velocity components is performed in the context of the log-law for which we follow the approach of Grotjans and Menter (1998) and Kuzmin et al. (2007), which is appropriate for finite element methods based on a weighted residual formulation. The traction BCs imposed on the tangential velocity are based on the horizontally averaged slip velocity that develops at the boundary and the law of the wall and is effected through the use of the mean-field eddy viscosity model of Sullivan et al. (1994). The approach originally used by Schumann (1975b) is used to convert the horizontally averaged traction to local values based on the local slip velocity in each of the horizontal directions.

In AMR-Wind, the periodic BCs in the x and y directions and the slip boundary on the top wall are applied. On the bottom wall, Dirichlet BCs in the normal direction and inhomogeneous Neumann BCs in the x and y directions are applied. The inhomogenous Neumann BC is set using the expression for τ, the total wall shear stress, and q, the total waIn ll heat flux, in (18). The stresses are specified at the terrain boundary following Moeng (1984). The wall stress vector is defined as

τ_{i 3} = \frac{\bar{u_{i}} s + \bar{s} (u_{i} - \bar{u_{i}})}{{\bar{s}}^{2}} u_{τ}^{2},

(18)

where u_i is the velocity at the first cell height, s is the wind speed

s = \sqrt{u_{1}^{2} + u_{2}^{2}}

, and u_τ is the friction velocity computed by using the Monin–Obukhov similarity law (Etling, 1996; Monin and Obukhov, 1954). The overbar indicates a horizontal plane average at the first cell height. The heat flux is defined as

q = [(θ - \bar{θ}) \bar{s} + (\bar{θ} - θ_{w}) s] \frac{u_{τ} κ}{\bar{s} ϕ_{h}},

(19)

where θ is the temperature, θ_w is the wall temperature, κ is the von Karman constant, and ϕ_h comes from the Monin–Obukhov similarity law.

The traction boundary conditions used in the two codes, NekRS and AMR-Wind, for the momentum equation along the lower wall were implemented in the same way as described in equations (18) and (19); the only difference is that the MFEV decomposition is used in NekRS but not in AMR-Wind.

The range of scales in these simulations is evident in Figure 2, which shows variations in potential temperature on a horizontal x-y slice at the height z = 100 m for resolutions Δx = 3.12 m, 1.56 m, and 0.78 m, respectively from left to right, for the differing codes/models. For AMR-Wind, the grid size is uniform in each direction. For NekRS, the element size is uniform in each direction; and Δx represents the average grid spacing given by 400 m/(E_∗N), where E_∗ is the number of elements in the x-, y-, or z-direction and N is the local polynomial order. The number of elements is $E = E_{*}^{3} = 16^{3}$ , 32³, and 64³ for the stated resolutions. The top row in Figure 2 shows NekRS results using MFEV and HPF; the middle row shows NekRS results using MFEV and Smagorinsky, as described in the Description of Codes section; and the bottom row shows AMR-Wind results. At a height of z = 100 m, the temperature variations around the horizontally averaged value are small, between 264.40 K and 264.80 K. One can see that as the grid scale is decreased from Δx = 3.12 m to 0.78 m, both codes capture increasingly finer scales. We remark that direct numerical simulation at the given Reynolds number, Re = 5 × 10⁷, would require $\approx 10^{15}$ grid points, which is two orders of magnitude beyond current state-of-the-art simulations of n = 18000³ for isotropic turbulence (Ravikumar et al., 2019). The importance of the SGS model is that it potentially allows one to account for the effects of small-scale motions without needing to resolve all of them.

Figure 3 shows profiles of the horizontally averaged streamwise, $〈 u 〉$ , and spanwise, $〈 v 〉$ , wind velocities (top) and potential temperature, $〈 θ 〉$ , (bottom) at 1-h time intervals between 1 h and 10 h for a mesh resolution of n = 512³ (Δx = 0.78 m) for each code. The left figures show NekRS results using MFEV and HPF, the center figures show NekRS using MFEV and Smagorinsky, and the right figures show the results for AMR-Wind. As can be observed, the time evolution of the mean velocity and temperature profiles obtained from the two codes agree well despite using very different numerical methods and subgrid-scale models. The agreement between AMR-Wind and NekRS improves when using MFEV and Smagorinsky (i.e., between the center and right figures). Specifically, the height of the low-level jet peak velocity during quasi-steady evolution in the GABLS problem (after approximately 7 h) is between 150 and 160 m, and its maximum value is between 9.5 and 9.7 m/s.

We remark that, while the initial focus for the current work was on co-development of performant GPU-based solvers for ABL simulations, the careful validation and cross-checking of the SGS models have also led to significant developments. For example, the improvement in the NekRS SGS model, which was realized by moving away from the HPF model in favor of the Smagorinsky model for the isotropic stress term, is a direct outcome of this collaborative effort. A more detailed exploration of resolution and model influence on the results is ongoing and will be the subject of a forthcoming paper.

4. Performance

Here we compare performance and tuning for the two codes. For each case, the codes use the same spatial resolution, Δx, and timestep size, Δt. Each code uses iteration tolerances of 10⁻⁴ and 10⁻⁶ (absolute, in nondimensional units) for the respective 2-norm residuals of the pressure-Poisson and velocity-Helmholtz problems. For purposes of timings, we use the solution at 6 h as an initial condition in each case in order to ensure that performance studies are done over a time frame in which the solutions have a representative turbulent flow. Table 1 provides a summary of the test parameters, in physical units, that are used for the strong- and weak-scaling studies. The spectral element cases use 8th-order polynomial basis (N = 8) with a number of grid points given by n = EN³. For these cases we take Δx to be the average grid spacing in each direction (i.e., 400 m/n^1/3). For the weak-scale study, the domain height is fixed at 400 m while the dimensions are increased in the x and y directions as n is increased. Mesh partitioning in NekRS is based on a parallel implementation of recursive spectral bisection (Pothen et al., 1990). Being a structured solver, AMR-Wind mesh partitioning is to divide the domain into equal-sized “boxes” of grid cells. Because mesh refinement is acheived by placing refined “levels” of mesh on top of less refined levels, this box partitioning can be used even with adaptive mesh refinement. An MPI rank then owns some number of boxes. For each box, work that can efficiently be done in parallel with GPU threads is offloaded to the GPU. Each GPU thread works on only a few cells of the box (Zhang et al., 2020). In order to avoid initial transient behavior, the average (wall) time per step, t_step, in seconds is measured over steps 101–200.

Table 1.

Problem setup for strong- and weak-scaling studies.

Strong-scaling test sets
Domain size	Grid points (n)	Δx (m)	Δt (s)
[400m]³	512 × 512 × 512	0.78	.062500
[400m]³	1024 × 1024 × 1024	0.39	.031250
[400m]³	2048 × 2048 × 2048	0.19	.015625
—	—	—	—

Weak-scaling test sets
Domain size	Grid points (n)	Δx (m)	Δt (s)
400 m × 400 m × 400 m	512 × 512 × 512	0.78	.0625
800 m × 800 m × 400 m	1024 × 1024 × 512	0.78	.0625
1600 m × 1600 m × 400 m	2048 × 2048 × 512	0.78	.0625
3200 m × 3200 m × 400 m	4096 × 4096 × 512	0.78	.0625

4.1. Performance tuning and profiling

We begin with performance optimization, profiling analysis, and CPU versus GPU comparisons. NekRS GPU performance tuning on Summit is demonstrated in detail in the work of (Fischer et al., 2022; Min et al., 2022). The base libParanumal kernels have their origins in the work of Warburton and co-workers (Chalmers et al., 2020; Świrydowicz et al., 2019; Chalmers and Warburton, 2020; Abdelfattah et al., 2021). A key algorithmic component is the Chebyshev-accelerated Schwarz-based p-multigrid for the pressure solve (Phillips et al., 2022), which is performed in 32-bit precision (e.g., as done by Fehn et al. (2018)) to reduce injection-bandwidth pressure on the Summit network interface cards. Communication for the nearest-neighbor communication (direct-stiffness summation in the finite element or spectral element context (Deville et al., 2002)) is overlapped with computation whenever it proves to be effective, which can yield as much as 10%–15% savings in NS applications. At the strong-scale limit of $\approx 2$ M points per GPU, there are enough points interior to each rank’s subdomain to balance out the communication overhead for the gather-scatter exchanges, at least at the fine-mesh level evaluations. For the coarser p-multigrid levels, one cannot always cover the communication with work. When initializing the communication kernels for each level of the p-multigrid solver, NekRS selects the fastest of several available communication strategies (e.g., overlapping, pack-on device or host, GPU direct or via the host), which are determined by timing tests during runtime setup. These unit tests also report the observed messaging bandwidth and thus provide insight into possible anomalous system behavior, which is useful when porting relatively new code to relatively new and unknown HPC platforms. For example, from existing logfiles, we were able to compare the observed bandwidth for several gather-scatter exchanges on NERSC’s Perlmutter platform before and after a network update from Cray’s Slingshot 10 to Slingshot 11, as shown below.

SS10:

pw+device MPI: 7.37e-05s / bi-bw: 18.2GB/s/rank

pw+device MPI: 5.16e-05s / bi-bw: 50.1GB/s/rank

pw+device MPI: 3.84e-05s / bi-bw: 16.8GB/s/rank

pw+host MPI: 2.46e-05s / bi-bw: 1.8GB/s/rank

SS11:

pw+device MPI: 4.38e-05s / bi-bw: 30.6GB/s/rank

pw+device MPI: 3.47e-05s / bi-bw: 74.4GB/s/rank

pw+device MPI: 2.74e-05s / bi-bw: 23.6GB/s/rank

pw+host MPI: 1.66e-05s / bi-bw: 2.7GB/s/rank

Here, SS10 indicates Slingshot 10, and SS11 indicates Slingshot 11, which shows about a 1.5× improvement over SS10. The listings also show which communication mode was used. We see that pw + device, which stands for pairwise device-to-device exchange (i.e., via GPU-direct), is used in most instances. The pw + host, which indicates the use of pairwise exchanges via the host, is used only in the case of many short messages, which is typically the scenario at the coarsest levels of the p-multigrid solver.

Over the course of the collaboration, AMR-Wind realized a 1.4× speedup with some improvements derived through AMReX library updates. The performance progress is demonstrated in Table 2, where the rows present a timing breakdown of a typical flow timestep. Advection involves predicting and forming the advection term using Godunov PPM WENO. MAC projection is a Poisson equation linear solve with a 7-point stencil that ensures that the face velocities are divergence free. The pressure solve is a Poisson equation linear solve with a 27-point stencil that approximately corrects the cell velocity to be divergence free at the end of the timestep. Velocity and scalar solve are Helmholtz equations with a 7-point stencil, and Fillpatch performs all communication within and across processors outside of the linear solver communication. In the table, the old version is AMR-Wind using AMReX from April 2021. The intermediate version is the same source code but with improvements to the linear solver settings. In particular, the components of the momentum equations are solved separately instead of as a coupled tensor solve. The velocity and scalar (temperature) linear systems are solved by using bi-conjugate gradient iteration instead of a full geometric multigrid approach. In Table 2 we see that these optimizations reduce the velocity solve time by almost 3× (.114 s to .039 s) and the scalar solve time by 1.5× (.035 s to .022 s). The new version is AMR-Wind based on AMReX from 2022 with the same improved linear system settings. Here, the scalar solve improves by another factor of 1.5. And the pressure solve is reduced from .0073 s to .0063 s per step.

Table 2.

AMR-Wind performance optimization.

	AMR-wind: Performance progress with AMReX library updates
	Last 200 steps averaged from 1000-step run on 8 nodes using n = 512³
	Old version		Intermmediate version		New verision
	(s)	(%)	(s)	(%)	(s)	(%)
Wall time per timestep	3.3200e-01	100	2.4100e-01	100	2.2500e-01	100
Advection	2.9814e-02	8.98	2.6920e-02	11.17	2.8687e-02	12.75
MAC projection	6.2582e-02	18.85	6.2660e-02	26.00	6.3135e-02	28.06
Pressure solve	7.3671e-02	22.19	7.3481e-02	30.49	6.3180e-02	28.08
Velocity solve	1.1401e-01	34.34	3.9669e-02	16.46	3.9307e-02	17.47
Scalar solve	3.4827e-02	10.49	2.2148e-02	9.19	1.5930e-02	7.08
Fillpatch	1.5538e-02	4.68	1.4990e-02	6.22	1.5030e-02	6.68

The computational cost contributors for AMR-Wind are shown graphically in Figure 4 in red colors. This figure shows the time spent performing the various equation solves and other operations during a timestep both as absolute time per timestep or percent of a timestep as a function of number of GPU processors used.

The key observations are as follows:

• The elliptic solves of the MAC projection and approximate pressure projection are leading computational cost contributors. Both solves use geometric multigrid with some differences. For instance, MAC projection is a cell-centered solve with a 7-point stencil whereas the approximate pressure projection is a node-centered solve with a 27-point stencil. This leads to differences in the coarsening and refining algorithms, smoothers, and boundary conditions.

• The MAC projection step in AMR-Wind is not completely necessary. It is used mainly in computing the advection term in the momentum and scalar predictors, which requires velocities on cell faces at the n + 1/2 time level. Godunov extrapolation is used to obtain these velocities, but then MAC projection makes them divergence free. This MAC projection step is not necessary but extends the stability range to CFL numbers up to about 2. Excluding the MAC projection step, the AMR-Wind computation time becomes more in line with NekRS. In fact, for the case using 4 GPU nodes on Crusher, excluding the MAC projection step would make AMR-Wind faster than NekRS. The same is not true on more nodes or on Summit, however.

• Figure 4 also shows that the percentage of total time spent on the elliptic solves increases as the number of nodes increases because the elliptic solve is highly communication dependent.

• The other linear system solves are the predictor for the velocity transport equation and the solves for scalar transport equations, such as potential temperature or subgrid-scale kinetic energy. There is a separate solve for each of the three components of the vector velocity equation. Each of these is a Helmholtz solve using a 7-point stencil. The solver is BiCG. Typically, when running atmospheric LES, there are five Helmholtz solves per timestep.

• Other function calls comprise source term calculations, boundary conditions, planar averaging, communication (excluding linear solve communication), linear solve setup, and copying solution arrays.

• Overall, Crusher is providing better performance than Summit. This is partly because there are more GPUs per node (8 versus 6) but also because AMR-Wind’s mesh partitioning strategy is more ammenable to GPU counts that are powers of two. For example, a problem size of 512³ is more easily partitioned by 8 GPUs/node than 6 GPUs/node. With 16 Crusher nodes (128 GPUs) a time per timestep of t_step = 0.11 s is achieved with AMR-Wind. With further scaling out with Summit the lowest time per timestep was 0.128 s on 128 Summit nodes (768 GPUs), as discussed below.

For NekRS, we start the GPU analysis with NVIDIA’s profiling tools. NVIDIA’s Nsight System summarizes kernel-level metrics for critical kernels. This tool identified that a bulk of total solution time for NekRS is spent evaluating the dealiased advection operator (subCycleStrongCubatureVolumeHex3D) both for the velocity vectors and for the temperature. Other leading consumers are the gather-scatter operations. Largely missing from this kernel-level analysis, however, is the time spent in the pressure preconditioner, which is separated across many kernels for the various levels of p-multigrid.

Each NekRS job tracks basic runtime statistics using a combination of MPI Wtime and cudaDeviceSynchronize or CUDA events. These are output every 500 timesteps unless the user specifies otherwise. From these, we collect aggregate timing breakdowns, roughly following the physical substeps of advection, pressure, and viscous thermal updates, plus tracking of known communication bottlenecks such as the pMG coarse-grid solve for the pressure preconditioner. Figure 4 shows the cost breakdown for this type of analysis over node counts ranging from 4 to 16. At lower node counts, the bulk of the NekRS time is spent in the makef and makeq (advection) routines, which are respectively responsible for setting up the right-hand sides of the viscous thermal updates for the momemtum and energy equations. To allow a larger CFL, the ABL simulations of viscous thermal updates use characteristics-based timestepping, which involves solving a sequence of hyperbolic subproblems on the interval [tⁿ⁻², tⁿ] (one for each velocity component and one for temperature) (Maday et al., 1990; Patel et al., 2019). Each subproblem takes several substeps using the dealiased advection operator, which performs quadrature on a 11 × 11 × 11 grid in each element. These substeps are thus compute intensive but not communication intensive, so they scale relatively well. The velocitySolve and scalarSolve, which involve communication-free diagonal preconditioning for conjugate gradient solution of (9)–(10), show similar scaling behavior. As with AMR-Wind, we see clearly in Figure 4 that the pressure solve does not scale as well as the other components.

We remark that Figure 4 indicates a significant amount of time is spent in udfExecuteStep. The majority of that cost results from the recently adopted mean-field eddy viscosity model (7), which requires several planar averages per timestep and is currently implemented as a user-defined function. For these calculations, which have low pressure and velocity iteration counts, the frequently called planar average utility has a significant impact on runtime (about 20%). Planar averaging is typically a postprocessing operation that is not performed on every step, but clearly it will need to be optimized in this LES application.

Table 3 gives a detailed breakdown of the per rank² kernel performance for NekRS in the 4-node case of Figure 4. The kernels are the advection subroutine (advSub), with either 3 components (velocity) or 1 (temperature); the Poisson/Helmholtz matrix-vector product (Ax) for the elliptic solves; and the fast-diagonalization method (fdm) for the Schwarz smoother (Lottes and Fischer, 2005). For each kernel, N indicates the polynomial order. For each kernel, the leading order complexity is CE_p(N + 1)⁴ for the number of floating-point operations and CE(N + 1)³ for the number of memory references, where E_p is the number of elements on each rank and C ≈ 12–30 is a kernel-dependent constant. FP indicates the working floating-point precision; GB/s the sustained streaming bandwidth on the device; GFLOPS the number of billions of floating-point operations per second sustained on the device for that particular kernel; and KV the kernel version identified as the fastest entry in each runtime benchmark test. The 32-bit precision kernels are used in the lower levels of p-multigrid. Because the Schwarz smoother operates on an extended domain, fdm executes on data that is extended to N + 2 in each direction compared with its corresponding Ax operation. The advection operation is dealiased (i.e., integration is on a finer mesh than the underlying velocity representation), so that the kernel also operates on a relatively large data set. We see that the work-intensive (high-N) and 32-bit kernels achieve impressive floating-point performance, well in excess of 1 TFLOPS (incidentally, the speed of ASCI Red, the world’s fastest computer just 25 years ago). NekRS also provides a conservative estimate of the overall FP64 floating-point rate per rank—here, close to 1 TFLOPS—which includes the message-passing overhead. (For this overall rate, each 32-bit operation is viscous thermal updates counted as half a flop).

Table 3.

NekRS runtime benchmark results associated with Figure 4.

Platform	Kernel	N	FP	GB/s	GFLOPS	KV
Summit	advSub (3)	10	64	613	3773	7
NVIDIA	advSub (1)	10	64	1137	2446	8
V100	Ax	8	64	844	1622	5
(p = 24)	Ax	8	64	900	1731	4
	Ax	8	64	901	1732	4
	Ax	8	32	859	3303	4
	Ax	4	64	832	975	5
	Ax	4	32	711	1667	6
	fdm	10	32	611	6210	3
	fdm	6	32	713	4422	3
Crusher	advSub (3)	10	64	491	3018	11
AMD	advSub (1)	10	64	868	1867	8
MI250X	Ax	8	64	662	1272	2
(p = 32)	Ax	8	64	736	1416	2
	Ax	8	64	736	1415	2
	Ax	8	32	742	2854	2
	Ax	4	64	708	830	6
	Ax	4	32	658	1543	0
	fdm	10	32	546	5551	4
	fdm	6	32	521	3234	4
Summit	Sustained		64		833
Crusher	Sustained		64		937

Figure 5 shows CPU and GPU strong-scaling performance for each code on Summit. The upper figures show standard time vs. node-count plots, which clearly indicate that it is easier to strong-scale on the CPU. On Summit, however, that point is moot given that one needs 128 nodes using a CPU-only configuration in order to get to the same time per step as using 4 nodes with 6 GPUs each (i.e., a factor of 32 difference in required node-hours to do the same work).

4.2. Strong- and weak-scaling performance

We next consider GPU-only performance on Summit using a single V100 per MPI rank. Figure 6, top, shows performance in terms of t_step for strong scaling as a function of the number of of GPUs, P, in the left column and as a function of number of points per rank, n/P, in the center column. Weak-scaling performance is presented in the right column. The wall-time figure also shows the ideal speedup curves scaling as P⁻¹. The lower plots show parallel efficiency,

P_{eff} : = \frac{t_{0} P_{0}}{t_{s t e p} P},

(20)

where P₀ is the smallest value of P that will hold the given problem and t₀ is the t_step value corresponding to P₀.

We see that at the lower resolution of n = 512³, the performance of the two codes is within a factor of 2 of each other out to p = 78. From the efficiency figures we can see that both curves have dropped below 80% efficiency by that point, so a more realistic point of comparison would be at p = 66 given that users would typically not run this relatively small case on p > 66. We note that p = 66 corresponds to n/P = 2M, which is a typical strong-scaling limit for NekRS on current-generation GPU platforms.

The center column in Figure 6 replots the strong-scaling information with n/P as the independent variable. Here we see a collapse of each code’s strong-scaling data into a single curve, particularly for NekRS. The efficiency plot, lower center, clearly shows the n/P = 2M mark as the 80% parallel efficiency point for NekRS. The AMR-Wind wall-time curves, upper center, are not as tightly grouped, particularly for the large problem sizes on large processor counts. It is tempting to speculate that this increased cost is due to an increase in iteration count, but Tables 4 and 5 show that is not the case, since each solver requires only two iterations per timestep for each of the problems. An important feature of AMR-Wind is that it generally performs better if P is a power of 2. At the critical point of n/P = 2M, NekRS is only a factor of 1.6 faster than AMR-Wind for the n = 512³ case.

Table 4.

NekRS GPU vs. AMR-Wind GPU strong-scaling performance study.

	Strong-scaling on Summit GPU, n = 512³, Δx = 0.78 m, Δt = 6.25e-2s, Ω = [400 m × 400 m × 400 m]
			NekRS						AMR-wind
Node	gpu	n/gpu	v _i	p _i	T _i	t _step	P _eff	r _t	v _i	p _i	T _i	t _step	P _eff	r _t
4	24	5.5924e+06	2	1.81	1	2.44e-01	100	3.90	2	2	2	3.19e-01	100	5.10
8	48	2.7962e+06	2	1.82	1	1.39e-01	87	2.22	2	2	2	2.37e-01	67	3.80
11	66	2.0336e+06	2	1.85	1	1.11e-01	79	1.78	2	2	2	1.79e-01	64	2.87
16	96	1.3981e+06	2	1.90	1	8.66e-02	70	1.38	2	2	2	1.75e-01	45	2.80
24	144	9.3207e+05	2	2.00	1	6.87e-02	59	1.09	2	2	2	1.60e-01	33	2.56
32	192	6.9905e+05	2	2.00	1	6.77e-02	45	1.08	2	2	2	1.46e-01	27	2.34
64	384	3.4953e+05	2	2.00	1	4.40e-02	34	0.70	2	2	2	1.43e-01	13	2.30
128	768	1.7476e+05	2	2.00	1	4.02e-02	18	0.64	2	2	2	1.28e-01	7.7	2.05
256	1536	8.7381e+04	2	2.00	1	3.60e-02	10	0.57	2	2	2	1.41e-01	3.5	2.26

	Strong-scaling on Summit GPU, n = 1024³, Δx = 0.39 m, Δt = 3.125e-2s, Ω = [400 m × 400 m × 400 m]
			NekRS						AMR-wind
node	gpu	n/gpu	v _i	p _i	T _i	t _step	P _eff	r _t	v _i	p _i	T _i	t _step	P _eff	r _t
32	192	5.5924e+06	1	1.4	1	2.34e-01	100	7.50	2	2	2	0.369	100	11.82
40	240	4.4739e+06	1	1.3	1	2.04e-01	91	6.54	2	2	2	0.402	73	12.86
50	300	3.5791e+06	1	1.4	1	1.72e-01	87	5.50	2	2	2	0.288	82	9.21
70	420	2.5565e+06	1	1.4	1	1.27e-01	84	4.06	2	2	2	0.303	56	9.72
80	480	2.2370e+06	1	1.3	1	1.15e-01	81	3.68	2	2	2	0.301	49	9.65
86	512	2.0972e+06		-		-	-		2	2	2	0.206	67	6.60
90	540	1.9884e+06	1	1.3	1	1.08e-01	76	3.47	2	2	2	0.217	60	6.95
100	600	1.7896e+06	1	1.3	1	9.57e-02	78	3.06	2	2	2	0.219	54	7.01

	Strong-scaling on Summit GPU, n = 2048³, Δx = 0.39 m, Δt = 1.5625e-2s, Ω = [400 m × 400 m × 400 m]
			NekRS						AMR-wind
node	gpu	n/gpu	v _i	p _i	T _i	t _step	P _eff	r _t	v _i	p _i	T _i	t _step	P _eff	r _t
256	1536	5.5924e+06		-		-	-		2	2	2	0.437	100	34.99
320	1920	4.4739e+06	1	1.16	1	2.10e-01	100	16.8	2	2	2	0.485	72	38.80
400	2400	3.5791e+06	1	1.20	1	1.80e-01	93	14.4	2	2	2	0.370	76	29.62
480	2880	2.9826e+06	1	1.22	1	1.54e-01	91	12.3	2	2	2	0.402	58	32.16
640	3840	2.2370e+06	1	1.25	1	1.28e-01	81	10.3	2	2	2	0.440	40	35.26
683	4096	2.0972e+06		-		-	-		2	2	2	0.321	51	25.69
800	4800	1.7896e+06	1	1.18	1	1.05e-01	80	8.4	2	2	2	0.390	36	31.23

Table 5.

NekRS GPU vs. AMR-Wind GPU weak-scaling performance study with fixed mesh density and resolution per GPU.

	Weak-scaling on Summit GPU, Δx = 0.78 m, Δt = 6.25e-2s
					NekRS						AMR-wind
Node	gpu	n	Ω	n/gpu	v _i	p _i	T _i	t _step	P _eff	r _t	v _i	p _i	T _i	t _step	P _eff	r _t
5	30	512²× 512	400 m² × 400 m	4.4739e+06	1	1.5	1	0.191	100	3.8	2	2	2	0.303	100	4.8
20	120	1024²× 512	800 m² × 400 m	4.4739e+06	1	1.7	1	0.200	95	4.0	2	2	2	0.326	93	5.2
80	480	2048²× 512	1600 m² × 400 m	4.4739e+06	1	1.9	1	0.218	87	3.8	2	2	2	0.344	88	5.5
320	1920	4096²× 512	3200 m² × 400 m	4.4739e+06	1	2.4	1	0.235	81	4.7	2	2	2	0.386	78	6.2

					NekRS						AMR-wind
Node	gpu	n	Ω	n/gpu	v _i	p _i	T _i	t _step	P _eff	r _t	v _i	p _i	T _i	t _step	P _eff	r _t
10	60	512²× 512	400 m² × 400 m	2.2370e+06	1	1.5	1	0.112	100	2.2	2	2	2	0.223	100	3.6
40	240	1024²× 512	800 m² × 400 m	2.2370e+06	1	1.7	1	0.118	95	2.4	2	2	2	0.231	96	3.7
160	960	2048²× 512	1600 m² × 400 m	2.2370e+06	1	2.1	1	0.127	88	2.5	2	2	2	0.269	83	4.3
640	3840	4096²× 512	3200 m² × 400 m	2.2370e+06	1	2.5	1	0.147	76	2.9	2	2	2	0.352	63	5.6

Figure 6, right, shows weak-scaling results for n/p = 2.2 M and 4.4 M. For the heavily loaded cases, AMR-Wind is within a factor of 1.6 of NekRS, but this figure increases to roughly a factor of 2 for the 2.2 M points-per-GPU case. The weak-scaling efficiency reaches 80% at around p = 2000 GPUs for all the cases save the AMR-Wind case with n/p = 2.2 M, which crosses the 80% mark at p ≈ 1100.

Tables 4 and 5 provide a detailed breakdown of several of the key metrics for the code performance, including iteration counts (v_i, p_i, T_i, for the respective velocity, pressure, and temperature iterative solvers), t_step, parallel efficiency (P_eff), and the wall-time to physical-time ratio (r_t). This last quantity is of particular interest since it must be smaller than unity for weather modeling applications. We also note that P is denoted by gpu in the tables. We see from Table 4 that, for a fixed value of n/P, r_t effectively doubles with each doubling of (linear) resolution. The reason for this increase is that the number of timesteps must also double whenever the number of points in each direction is doubled (for fixed domain size). Throughout the table, we see that roughly two iterations are required per timestep for each of the linear solvers, indicating that the preconditioners are robust with respect to mesh size, although NekRS does show some increase in iteration count in the weak-scaling results.

We remark that AMR-Wind supports block-structured adaptive mesh refinement, which means that static grids do not leverage one of its main features. It is nonetheless highly performant on this problem. Moreover, AMR-Wind has a significant performance boost when the number of ranks is a power of 2, as seen in Table 4 for the n = 1024³ case for p = 512 and in the n = 2048³ case for p = 4096. In the former case, the parallel efficiency jumps from 49% to 67% as P changes from 480 to 512. In the latter, it jumps from 40% for p = 3840 to 51% for p = 4096. These performance gains derive from the block decompositions used in AMR-Wind, which favor block sizes (and, thus, processor counts) that are powers of 2.

We close with a scaling comparison of Summit and Crusher performance for NekRS in Figure 7. The upper figures show standard strong scaling as a function of the number of ranks on the left (one GPU or GCD per rank) and as a function of n/P on the right. The lower plots show the timing for the makef kernel (left), which evaluates the nonlinear advection term and does not require communication, and for the coarse-grid solve (right), which is communication dominated. The coarse-grid problem, which has roughly E degrees of freedom (with E = 262144 in this case), is solved by using algebraic multigrid (hypre) on the host CPUs (one CPU per MPI rank). The performance for these two platforms is remarkably similar.

5. Conclusion

We presented detailed profiling and timing results for two CFD codes, NekRS and AMR-Wind, applied to the GABLS atmospheric boundary layer test problem, which is of direct relevance to wind farm modeling and weather forecasting. Strong and weak scaling were demonstrated on up to p = 4800 NVIDIA V100 GPUs on OLCF’s Summit. For NekRS, wall-clock times of about 0.1 s per timestep were observed for n/P = 2M, which is the 80% efficiency point across a range of problem sizes. For the 512³, 1024³, and 2048³ resolutions, 0.1 s per timestep corresponds to ratios of wall-clock to simulated time of 1.8, 3.6, and 7.1, respectively. AMR-Wind was generally within a factor of 1.4–2.0 of the performance of NekRS over the range of interest, which is partly attributal to the fact that the incompressibility constraint for AMR-Wind requires two elliptic-system solves per timestep whereas NekRS requires only one. For both codes, strong scaling is ultimately inhibited by these intrinsically communication-intensive elliptic solves. For NekRS it was shown that a single GCD of the MI250X on Crusher is delivering performance that is comparable to a single V100 on Summit.

Simulation differences between NekRS and AMR-Wind were observed to follow the expected results for high- versus low-order discretizations, namely, that at a given refinement level the high-order case retains more high-frequency signal than does its low-order counterpart. We note that the second-order AMR-Wind algorithms were chosen for their robustness in addressing two challenging problems: two-phase air-water flows where the density ratio at the interface is a thousand (Kuhn et al., 2023) and coupling with an unstructured-grid solver through overset meshes (Sharma et al., 2024).

While both codes solve the acoustically incompressible Navier–Stokes equations, the codes use fundamentally different discretizations, different linear-system solvers, and different GPU-abstraction layers. Despite these differences, the strong-scaling results demonstrate the same fundamental opportunities and challenges on GPUs. First, if there are enough degrees of freedom on each GPU rank, GPU simulations can be dramatically faster than CPUs on a per-node basis. However, whereas both codes show a more typical (and good) range of ideal strong scaling on CPUs, GPU strong-scaling efficiency quickly degrades with the move to more GPU ranks. As described, this poor strong-scaling performance is tied to the elliptic-system linear solves and shows that further research is required to improve strong scaling of incompressible-flow computational fluid dynamics on GPUs.

In this work we performed careful subgrid-scale model and wall model comparisons between the two codes to realize comparable results. We employed state-of-the-art, physically accurate models that are currently available in the literature. More detailed modeling questions are part of an ongoing investigation and will be presented in a future paper.

Footnotes

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This research was supported in part by the Exascale Computing Project (17-SC-20-SC), a joint project of the U.S. Department of Energy Office of Science and the National Nuclear Security Administration, responsible for delivering a capable exascale ecosystem, including software, applications, and hardware technology, to support the nation’s exascale computing imperative. Funding was also provided by the U.S. Department of Energy, Office of Energy Efficiency and Renewable Energy, Wind Energy Technologies Office. This work was authored in part by the National Renewable Energy Laboratory, operated by Alliance for Sustainable Energy, LLC, for the U.S. Department of Energy (DOE) under Contract No. DE-AC36-08GO28308. The research used resources at the Oak Ridge Leadership Computing Facility at Oak Ridge National Laboratory, which is supported by the Office of Science of the U.S. Department of Energy under Contract DE-AC05-00OR22725. The research also used computational resources sponsored by the DOE Office of Energy Efficiency and Renewable Energy and located at the National Renewable Energy Laboratory. The views expressed in the article do not necessarily represent the views of the DOE or the U.S. Government. The U.S. Government retains and the publisher, by accepting the article for publication, acknowledges that the U.S. Government retains a nonexclusive, paid-up, irrevocable, worldwide license to publish or reproduce the published form of this work, or allow others to do so, for U.S. Government purposes.

ORCID iD

Misun Min

Notes

Author biographies

Misun Min is a Computational Mathematician in the Mathematics and Computer Science Division at Argonne National Laboratory. She holds a Ph.D. in Applied Mathematics from Brown University, where she studied spectral methods for PDEs. Dr. Min is the Argonne PI for the Center for Efficient Exascale Discretizations, which is supported as part of DOE’s Exascale Computing Project. Her research focuses on high-order methods for challenging PDE-based applications. A major part of her effort is directed towards performance and scalability, whether using millions of CPUs or tens of thousands of GPUs. Dr. Min’s software developments include spectral element codes for electromagnetics (NekCEM), lattice-Boltzmann methods (NekLBM), and fluid-thermal transport (Nek5000/RS). She developed the first GPU-enabled variant of the Nek code suite using OpenACC for NekCEM simulations on OLCF’s Titan. Dr. Min has led the scaling efforts for Nek5000/RS, with problems exceeding two billion spectral elements on DOE’s leadership computers.

Michael Brazell is a Computational Scientist at the National Renewable Energy Laboratory (NREL) working in the Computational Science Center. He holds a Ph.D. in Mechanical Engineering from Clarkson University where he studied higher-order finite-element methods for computational fluid dynamics. As a postdoctoral researcher at the CFD Laboratory at the University of Wyoming, Dr. Brazell continued research along these lines applying higher-order methods to wind turbine and rotorcraft flows. At NREL, he is a key developer of the AMR-Wind code funded by the U.S. Department of Energy (DOE) Exascale Computing Project and the DOE Office of Energy Efficiency and Renewable Energy. He has also worked with NREL’s advanced wind turbine and wind farm controls group to couple real time wind data with CFD for controls applications.

Ananias Tomboulides is a Professor of Mechanical Engineering at the Aristotle University of Thessaloniki and also an Argonne affiliate. He has a long-standing collaboration with the Mathematics and Computer Science Division at Argonne National Laboratory and is one of the main developers of the open-source CFD code Nek5000/NekRS, with contributions in the areas of low Mach number combustion, multiphase boiling flows, high-order RANS models and moving mesh developments for the simulation of flows in internal combustion engines. He has served as Chairman of ERCOFTAC (the European Research Community On Flow, Turbulence And Combustion) and is a member of the Editorial Board of the Journal Flow Turbulence and Combustion (FTAC). His group has been involved in two European Centers of Excellence funded by EuroHPC, one in Exascale Combustion (CoEC) and another one in Exascale CFD (CEEC).

Matthew Churchfield is a Distinguished Member of the Research Staff at the National Renewable Energy Laboratory (NREL) working in the National Wind Technology Center. He holds a Ph.D. in aeronautical engineering from Purdue University where he studied turbulence modeling for aircraft trailing vortex flow fields. Dr. Churchfield specializes in computational fluid dynamics of wind farm and atmospheric flows, focusing mainly on large-eddy simulations. He currently is the NREL PI for a U.S. Department of Energy (DOE) project involving hurricane boundary layers and hurricane impacts to offshore wind farms and in the recent past, he was the NREL PI for a DOE project developing method to couple mesoscale and microscale atmospheric flow models. Dr. Churchfield was a chief member of NREL’s effort to develop its first-generation wind-farm large-eddy simulation tool, the Open FOAM-based Simulator fOr Wind Farm Applications (SOWFA). More recently, he has been involved in the multi-laboratory effort to develop the ExaWind next-generation wind farm simulation tools, which was supported by the DOE Exascale Computing Project and the DOE Office of Energy Efficiency and Renewable Energy.

Paul Fischer is a Blue Water Professor in Computer Science and in Mechanical Science and Engineering at the University of Illinois, Urbana Champaign, and is also a senior scientist at Argonne National Laboratory. He has pioneered the development of spectral element methods for high-performance simulations of turbulence, including the development of Nekton 2.0, which was the first commercial software for distributed-memory parallel computers. The research variant, Nek5000, is a prior Gordon Bell winner and has scaled to millions of ranks on Mira and Sequoia. It is used by over 500 researchers in industry and academia. Fischer’s current research is focused on advanced preconditioners for GPU-based solutions of PDEs and reduced-order models for turbulent flows with applications to industrial problems.

Michael Sprague is a Chief Wind Computational Scientist and Distinguished Member of Research Staff at the National Renewable Energy Laboratory (NREL) near Boulder, Colorado. Dr. Sprague has been Principal Investigator for several U.S. Department of Energy (DOE) supported projects in high-fidelity modeling and high-performance computing for wind energy, including ExaWind, which was part of the DOE Exascale Computing Project. He is currently Director of the DOE Office of Science FLOWMAS Energy Earthshot Research Center. Before coming to NREL in 2010, Dr. Sprague was an assistant professor and founding faculty member in applied mathematics at the University of California, Merced. He graduated with a PhD in Mechanical Engineering from the University of Colorado at Boulder and Bachelor of Science degree from the University of Wisconsin-Madison.

References

Abdelfattah

Barra

Beams

, et al. (2021) GPU algorithms for efficient exascale discretizations. Parallel Computing 108: 102841.

Alexander

Almgren

Bell

, et al. (2020) Exascale applications: skin in the game. Philosophical transactions. Series A, Mathematical, physical, and engineering sciences 378(2166): 20190056.

Almgren

Bell

Colella

, et al. (1998) A conservative adaptive projection method for the variable density incompressible Navier-Stokes equations. Journal of Computational Physics 142(1): 1–46. DOI: 10.1006/jcph.1998.5890.

Almgren

Bell

Crutchfield

(2000) Approximate projection methods: Part I. Inviscid Analysis. SIAM Journal on Scientific Computing 22(4): 1139–1159. DOI: 10.1137/S1064827599357024

Arndt

Fehn

Kanschat

, et al. (2020) ExaDG: high-order discontinuous Galerkin for the exa-scale. In: Bungartz

Reiz

Uekermann

, et al. (eds) Software for Exascale Computing - SPPEXA 2016-2019. Berlin: Springer, 189–224.

Beare

Macvean

Holtslag

, et al. (2006) An intercomparison of large-eddy simulations of the stable boundary layer. Boundary-Layer Meteorology 118(2): 247–272.

Bell

Howell

Colella

(1991) An efficient second-order projection method for viscous incompressible flow. 10th computational fluid dynamics conference, Honolulu,HI, USA, 24 June 1991

Berg

Kelly

(2020) Atmospheric turbulence modelling, synthesis, and simulation Wind Energy Modeling and Simulation: Volume 1: Atmosphere and Plant 1(34): 183–216.

Chalmers

Warburton

(2020) Portable High-Order Finite Element Kernels I: Streaming Operations. CoRR abs/2009. https://arxiv.org/abs/2009.10917.

10.

Chalmers

Karakus

Austin

, et al. (2020) Libparanumal: A Performance Portable High-Order Finite Element Library. https://github.com/paranumal/libparanumal

11.

Churchfield

Moriarty

(2020) Modeling and simulation of wind-farm flows. Wind Energy Modeling and Simulation: Volume 1: Atmosphere and Plant 125: 217–271.

12.

Colella

Woodward

(1984) The piecewise parabolic method (ppm) for gas-dynamical simulations. Journal of Computational Physics 54(1): 174–201. DOI: 10.1016/0021-9991(84)90143-8

13.

Deville

Fischer

Mund

(2002) High-order Methods for Incompressible Fluid Flow. Cambridge: Cambridge University Press.

14.

ECP (2017) Exascale Computing Project. https://www.exascaleproject.org/.

15.

Etling

(1996) Modelling the vertical ABL structure. Modelling of Atmospheric Flow Fields 45.

16.

Fehn

Wall

Kronbichler

(2018) Efficiency of high-performance discontinuous Galerkin spectral element methods for under-resolved turbulent incompressible flows. International Journal for Numerical Methods in Fluids 88: 32–54. DOI: 10.1002/fld.4511.

17.

Fischer

(1998) Projection techniques for iterative solution of with successive right-hand sides Ax = b with successive right-hand sides. Computer Methods in Applied Mechanics and Engineering 163: 193–204.

18.

Fischer

Lottes

(2004) Hybrid Schwarz-multigrid methods for the spectral element method: extensions to Navier-Stokes. In: Kornhuber

Hoppe

Périaux

, et al. (eds) Domain Decomposition Methods in Science and Engineering Series. Berlin: Springer.

19.

Fischer

Schmitt

Tomboulides

(2017) Recent developments in spectral element simulations of moving-domain problems, volume 79. Fields Institute Communications. Berlin: Springer, pp. 213–244.

20.

Fischer

Merzari

Min

, et al. (2021) Highly Optimized Full-Core Reactor Simulations on Summit. arXiv preprint arXiv:2110.01716

21.

Fischer

Kerkemeier

Min

, et al. (2022) NekRS, a GPU-accelerated spectral element Navier-Stokes solver. Parallel Computing 114: 102982.

22.

Grotjans

Menter

(1998) Wall functions for general application CFD codes. In: ECCOMAS

(ed) Proc. 4th European Computational Fluid Dynamics Conference. Hoboken, NJ: John Wiley & Sons, 1112.

23.

Huismann

Stiller

Fröhlich

(2019) Scaling to the stars – a linearly scaling elliptic solver for p-multigrid. Journal of Computational Physics 398: 108868.

24.

Kosović

Curry

(2000) A large eddy simulation study of a quasi-steady, stably stratified atmospheric boundary layer. Journal of the Atmospheric Sciences 57(8): 1052–1068.

25.

Krank

Fehn

Wall

, et al. (2017) A high-order semi-explicit discontinuous Galerkin solver for 3D incompressible flow with application to DNS and LES of turbulent channel flow. Journal of Computational Physics 348(1): 634–659. DOI: 10.1016/j.jcp.2017.07.039.

26.

Kuhn

Deskos

Sprague

(2023) A mass-momentum consistent coupling for mesh-adaptive two-phase flow simulations. Computers & Fluids 252: 105770.

27.

Kuzmin

Mierka

Turek

(2007) On the implementation of the κ-ε turbulence model in incompressible flow solvers based on a finite element discretisation. International Journal of Computing Science and Mathematics 1: 193–206. DOI: 10.1504/IJCSM.2007.016531.

28.

Lilly

(1962) On the numerical simulation of buoyant convection. Tellus A: Dynamic Meteorology and Oceanography 14(2): 148–172.

29.

Lindquist

Min

Fischer

(2021) Scalable Interpolation on GPUs for Thermal Fluids Applications. Lemont, IL: Argonne National Laboratory. Technical Report ANL-21/55.

30.

Lottes

Fischer

(2005) Hybrid multigrid/Schwarz algorithms for the spectral element method. Journal of Scientific Computing 24: 45–78.

31.

Maday

Patera

Rønquist

(1990) An operator-integration-factor splitting method for time-dependent problems: application to incompressible fluid flow. Journal of Scientific Computing 5: 263–292.

32.

Medina

St-Cyr

Warburton

(2014) OCCA: A Unified Approach to Multi-Threading Languages. arXiv preprint arXiv:1403.0968

33.

Min

Lan

Fischer

, et al. (2022) Optimization of full-core reactor simulations on Summit. SC22: International Conference for High Performance Computing, Networking, Storage and Analysis. Piscataway, NJ: IEEE, 1–11.

34.

Mirocha

(2020) Mesoscale to microscale coupling for high-fidelity wind plant simulation. Wind Energy Modeling and Simulation: Volume 1: Atmosphere and Plant. London, UK: IET Digital Library, 117–181.

35.

Moeng

(1984) A large-eddy-simulation model for the study of planetary boundary-layer turbulence. Journal of the Atmospheric Sciences 41(13): 2052–2062.

36.

Monin

Obukhov

(1954) Basic laws of turbulent mixing in the surface layer of the atmosphere. Tr. Akad. Nauk. SSSR Geophiz. Inst 24(151): 163—–187.

37.

Motheau

Wakefield

(2020) Investigation of finite-volume methods to capture shocks and turbulence spectra in compressible flows. Communications in Applied Mathematics and Computational Science 15(1): 1–36.

38.

Orszag

(1980) Spectral methods for problems in complex geometries. Journal of Computational Physics 37: 70–92.

39.

Patel

Fischer

Min

, et al. (2019) A characteristic-based, spectral element method for moving-domain problems. Journal of Scientific Computing 79: 564–592.

40.

Patera

(1984) A spectral element method for fluid dynamics: laminar flow in a channel expansion. Journal of Computational Physics 54: 468–488.

41.

Pedersen

Gryning

Kelly

(2014) On the structure and adjustment of inversion-capped neutral atmospheric boundary-layer flows: large-eddy simulation study. Boundary-Layer Meteorology 153(1): 43–62. DOI: 10.1007/s10546-014-9937-z

42.

Phillips

Kerkemeier

Fischer

(2022) Tuning spectral element preconditioners for parallel scalability on GPUs. Proceedings of the 2022 SIAM conference on parallel processing for scientific computing, Seattle, WA, USA, 23–26 February 2022.

43.

Pothen

Simon

Liou

(1990) Partitioning sparse matrices with eigenvectors of graphs. SIAM Journal on Matrix Analysis and Applications 11: 430–452.

44.

Ravikumar

Appelhans

Yeung

(2019) Acceleration of extreme scale pseudo-spectral simulations of turbulence using asynchronism. In: Proceedings of the international conference for high performance computing, networking, storage and analysis, Denver, CO, USA, 17–22 November 2019

45.

Robinson

Sprague

(2020) Looking forward: the promise and challenge of exascale computing. Wind Energy Modeling and Simulation: Volume 1: Atmosphere and Plant. London: IET Digital Library, 1–22.

46.

Sauer

Muñoz-Esparza

(2020) The FastEddy® resident‐GPU accelerated large‐eddy simulation framework: model formulation, dynamical‐core validation and performance benchmarks. Journal of Advances in Modeling Earth Systems 12(11): e2020MS002100.

47.

Schalkwijk

Jonker

Siebesma

, et al. (2015) Weather forecasting using GPU-based large-eddy simulations. Bulletin of the American Meteorological Society 96(5): 715–723. DOI: 10.1175/BAMS-D-14-00114.1.

48.

Schumann

(1975a) Subgrid scale model for finite difference simulations of turbulent flows in plane channels and annuli. Journal of Computational Physics 18: 376–404.

49.

Schumann

(1975b) Subgrid scale model for finite difference simulations of turbulent flows in plane channels and annuli. Journal of Computational Physics 18(4): 376–404. DOI: 10.1016/0021-9991(75)90093-5.

50.

Sharma

Brazell

Vijayakumar

, et al. (2024) ExaWind: open-source CFD for hybrid-RANS/LES geometry-resolved wind turbine simulations in atmospheric flows. Wind Energy 27(3): 225–257. DOI: 10.1002/we.2886.

51.

Smagorinsky

(1963) General circulation experiments with the primitive equations. Monthly Weather Review 91(3): 99–164. DOI: 10.1175/1520-0493(1963)091<0099:GCEWTP>2.3.CO;2.

52.

Sprague

Boldyrev

Fischer

, et al. (2017) Turbulent Flow Simulation at the Exascale: Opportunities and Challenges Workshop. Technical report, Washington, DC: U.S. Department of Energy, Office of Science, Advanced Scientific Computing Research. https://www.nrel.gov/docs/fy17osti/67648.pdf.PublishedasTech.Rep.NREL/TP-2C00-67648bytheNationalRenewableEnergyLaboratory

53.

Sprague

Ananthan

Vijayakumar

, et al. (2020) ExaWind: a multifidelity modeling and simulation environment for wind energy. Journal of Physics: Conference Series 1452: 012071. DOI: 10.1088/1742-6596/1452/1/012071.

54.

Stolz

Schlatter

Kleiser

(2005) High-pass filtered eddy-viscosity models for large-eddy simulations of transitional and turbulent flow. Physics of Fluids 17(6): 065103. DOI: 10.1063/1.1923048.

55.

Sullivan

McWilliams

Moeng

(1994) A subgrid-scale model for large-eddy simulation of planetary boundary-layer flows. Boundary-Layer Meteorology 71: 247–276.

56.

Sullivan

Edson

Hristov

, et al. (2008) Large-eddy simulations and observations of atmospheric marine boundary layers above nonequilibrium surface waves. Journal of the Atmospheric Sciences 65(4): 1225–1245.

57.

Świrydowicz

Chalmers

Karakus

, et al. (2019) Acceleration of tensor-product operations for high-order finite element methods. The International Journal of High Performance Computing Applications 33(4): 735–757.

58.

TOP500 (2023). https://www.top500.org/.

59.

Van Heerwaarden

Van Stratum

Heus

, et al. (2017) MicroHH 1.0: a computational fluid dynamics code for direct numerical simulation and large-eddy simulation of atmospheric boundary layer flows. Geoscientific Model Development 10(8): 3145–3165.

60.

Van Leer

(1977) Towards the ultimate conservative difference scheme. iv. a new approach to numerical convection. Journal of Computational Physics 23(3): 276–299. DOI: 10.1016/0021-9991(77)90095-X.

61.

Zhang

Almgren

Beckner

, et al. (2019) AMReX: a framework for block-structured adaptive mesh refinement. Journal of Open Source Software 4(37): 1370. DOI: 10.21105/joss.01370.

62.

Zhang

Myers

Gott

, et al. (2020) AMReX: Block-Structured Adaptive Mesh Refinement For Multiphysics Applications. arXiv Preprint. https://arxiv.org/abs/2009.12009