CI

A high-performance C++ solver for incompressible turbulent flow with pluggable turbulence closures ranging from classical algebraic models to advanced transport equations and data-driven neural networks. Features a fractional-step projection method with multiple Poisson solvers, pure C++ NN inference, and comprehensive GPU acceleration via OpenMP target offload.

• Fractional-step projection method for incompressible Navier-Stokes
  • Explicit Euler time integration with adaptive CFL-based time stepping
  • Multiple Poisson solvers with automatic selection (FFT, Multigrid, HYPRE)
  • Pressure projection for divergence-free velocity
• Staggered MAC grid with second-order central finite differences
• 10 turbulence closures: algebraic, transport, EARSM, and neural network models
• Pure C++ NN inference - no Python/TensorFlow at runtime
• GPU acceleration via OpenMP target directives for NVIDIA and AMD GPUs

• Quick Start
• Governing Equations
• Numerical Methods
• Boundary Conditions
• Poisson Solvers
• Turbulence Closures
• Supported Flow Configurations
• Command-Line Options
• GPU Acceleration
• Validation
• References

mkdir build && cd build
cmake .. -DCMAKE_BUILD_TYPE=Release
make -j4

# Laminar channel flow (Poiseuille, analytical validation)
./channel --Nx 32 --Ny 64 --nu 0.01 --adaptive_dt --max_steps 10000
# Turbulent channel with SST k-omega
./channel --Nx 64 --Ny 128 --Re 5000 --model sst --adaptive_dt
# Neural network turbulence model
./channel --model nn_tbnn --nn_preset tbnn_channel_caseholdout --adaptive_dt
# 3D Taylor-Green vortex
./taylor_green_3d --Nx 64 --Ny 64 --Nz 64 --Re 100 --max_steps 1000

The solver implements the incompressible Reynolds-Averaged Navier-Stokes (RANS) equations:

$$\frac{\partial \bar{u}_i}{\partial t} + \bar{u}_j \frac{\partial \bar{u}_i}{\partial x_j} = -\frac{1}{\rho} \frac{\partial \bar{p}}{\partial x_i} + \frac{\partial}{\partial x_j}\left[(\nu + \nu_t) \frac{\partial \bar{u}_i}{\partial x_j}\right] + f_i$$

$$\nabla \cdot \mathbf{u} = 0$$

Variables:

The solver uses a three-step fractional-step method (Chorin 1968) to decouple pressure and velocity:

Step 1: Provisional Velocity (explicit momentum without pressure)

$$\frac{\mathbf{u}^* - \mathbf{u}^n}{\Delta t} = -(\mathbf{u}^n \cdot \nabla)\mathbf{u}^n + \nabla \cdot [(\nu + \nu_t) \nabla \mathbf{u}^n] + \mathbf{f}$$

Step 2: Pressure Poisson Equation (enforce incompressibility)

$$\nabla^2 p' = \frac{1}{\Delta t} \nabla \cdot \mathbf{u}^*$$

Step 3: Velocity Correction (project onto divergence-free space)

$$\mathbf{u}^{n+1} = \mathbf{u}^* - \Delta t \nabla p'$$

This ensures $\nabla \cdot \mathbf{u}^{n+1} = 0$ to machine precision.

All spatial derivatives use second-order central finite differences on a staggered Marker-and-Cell (MAC) grid:

Gradient (central difference):

$$\left.\frac{\partial u}{\partial x}\right|_{i,j} = \frac{u_{i+1,j} - u_{i-1,j}}{2\Delta x}$$

Laplacian (5-point stencil in 2D, 7-point in 3D):

$$\left.\nabla^2 u\right|_{i,j} = \frac{u_{i+1,j} - 2u_{i,j} + u_{i-1,j}}{\Delta x^2} + \frac{u_{i,j+1} - 2u_{i,j} + u_{i,j-1}}{\Delta y^2}$$

Convective Schemes (selectable via convective_scheme config):

• Central differences (default): Second-order accurate, lower numerical dissipation
• First-order upwind: More stable at high Reynolds numbers, increased numerical dissipation

Variable Viscosity Diffusion:

$$\nabla \cdot [(\nu + \nu_t) \nabla u] = \frac{1}{\Delta x^2}\left[\nu_{e}(u_{i+1,j} - u_{i,j}) - \nu_{w}(u_{i,j} - u_{i-1,j})\right] + \ldots$$

where $\nu_e, \nu_w$ are face-averaged effective viscosities.

Explicit Euler with adaptive time stepping:

$$\Delta t = \text{CFL} \cdot \min\left(\frac{\Delta x}{|u|_{\max}}, \frac{\Delta y}{|v|_{\max}}, \frac{\Delta z}{|w|_{\max}}\right)$$

• CFL number: Default 0.5 (configurable via --CFL)
• Steady-state convergence: Iterate until $|\mathbf{u}^{n+1} - \mathbf{u}^n|_\infty < \text{tol}$

Note: Inflow/Outflow boundary conditions are defined in the code interface but not yet fully implemented for all solvers.

The pressure Poisson equation supports three BC types:

Standard BC Configurations:

For problems with all Neumann or periodic pressure boundaries (no Dirichlet BC), the pressure is underdetermined up to a constant. The solver automatically:

1. Detects this condition via has_nullspace() check
2. Subtracts the mean pressure after each solve to fix the gauge

The solver provides 6 Poisson solver options with automatic selection based on grid configuration and boundary conditions:

FFT (3D) → FFT2D (2D) → FFT1D (3D partial-periodic) → HYPRE → Multigrid

The default solver implements a geometric multigrid V-cycle:

1. Pre-smooth: Apply smoothing iterations on fine grid (Chebyshev or Jacobi)
2. Restrict: Compute residual and transfer to coarse grid (full weighting)
3. Recurse: Solve on coarse grid (recursively)
4. Prolongate: Interpolate correction back to fine grid (bilinear)
5. Post-smooth: Apply smoothing iterations

Features:

• O(N) complexity (optimal)
• 5-15 V-cycles to convergence (vs 1000-10000 SOR iterations)
• CUDA Graph optimization: Entire V-cycle captured as single GPU graph (NVHPC compilers)
• Chebyshev polynomial smoother (faster than Jacobi)

Convergence Criteria (any triggers exit):

• tol_rhs: RHS-relative $|r|/|b| < \epsilon$ (recommended for projection)
• tol_rel: Initial-residual relative $|r|/|r_0| < \epsilon$
• tol_abs: Absolute $|r|_\infty < \epsilon$

For problems with periodic boundaries, FFT solvers provide spectral accuracy:

• FFT (3D): 2D FFT in x-z + batched tridiagonal solve in y (cuSPARSE)
• FFT2D: 1D FFT in x + batched tridiagonal in y
• FFT1D: 1D FFT in periodic direction + 2D Helmholtz solve per mode

GPU-accelerated parallel semicoarsening multigrid from HYPRE:

• Supports uniform AND stretched grids
• Entire solve runs on GPU via CUDA backend
• Automatic download and build via CMake FetchContent

See docs/POISSON_SOLVER_GUIDE.md for detailed documentation.

The solver supports 10 turbulence closure options:

Classical model with van Driest wall damping:

$$\nu_t = (\kappa y)^2 |\mathbf{S}| \left(1 - e^{-y^+/A^+}\right)^2$$

• $\kappa = 0.41$ (von Kármán constant)
• $A^+ \approx 26$ (van Driest damping constant)
• $|\mathbf{S}| = \sqrt{2S_{ij}S_{ij}}$ (strain rate magnitude)

Symbolic regression formula discovered by genetic algorithms (Weatheritt & Sandberg 2016):

$$\nu_t = f_{\text{GEP}}(S_{ij}, \Omega_{ij}, y, \text{Re}_\tau, \ldots)$$

Menter's Shear Stress Transport model (1994):

k-equation: $$\frac{\partial k}{\partial t} + \bar{u}_j \frac{\partial k}{\partial x_j} = P_k - \beta^* k \omega + \nabla \cdot [(\nu + \sigma_k \nu_t) \nabla k]$$

ω-equation (with cross-diffusion): $$\frac{\partial \omega}{\partial t} + \bar{u}j \frac{\partial \omega}{\partial x_j} = \alpha \frac{\omega}{k} P_k - \beta \omega^2 + \nabla \cdot [(\nu + \sigma\omega \nu_t) \nabla \omega] + CD_\omega$$

Eddy viscosity: $$\nu_t = \frac{a_1 k}{\max(a_1 \omega, S F_2)}$$

• Blending functions F1, F2 for k-ε/k-ω transition
• Production limiter for numerical stability
• Wall boundary conditions: k = 0, ω = ω_wall(y)

Wilcox (1988) formulation without blending:

$$\nu_t = \frac{k}{\omega}$$

EARSM models predict the full Reynolds stress anisotropy tensor using a tensor basis expansion:

$$b_{ij} = \sum_{n=1}^{N} G_n(\eta, \xi) , T_{ij}^{(n)}(\mathbf{S}, \mathbf{\Omega})$$

where:

• $b_{ij}$ = anisotropy tensor (traceless)
• $T_{ij}^{(n)}$ = integrity basis tensors
• $G_n$ = scalar coefficient functions
• $\eta = Sk/\epsilon$, $\xi = \Omega k/\epsilon$ = normalized invariants

Combined with SST k-ω transport for k and ω evolution.

Most sophisticated variant with cubic implicit equation for realizability.

Quadratic model without implicit solve.

Classical weak-equilibrium model using first 3 basis tensors.

Re_t-Based Blending:

EARSM models use smooth blending between linear Boussinesq (laminar) and full nonlinear (turbulent):

$$\alpha(\text{Re}_t) = \frac{1}{2}\left(1 + \tanh\left(\frac{\text{Re}_t - \text{Re}_{t,\text{center}}}{\text{Re}_{t,\text{width}}}\right)\right)$$

where $\text{Re}_t = k/(\nu\omega)$. Default transition: center at Re_t = 10, width = 5.

Multi-layer perceptron for scalar eddy viscosity:

$$\nu_t = \text{NN}_{\text{MLP}}(\lambda_1, \ldots, \lambda_5, y/\delta)$$

Inputs (invariants of strain and rotation tensors):

• $\lambda_1 = S_{ij}S_{ij}$, $\lambda_2 = \Omega_{ij}\Omega_{ij}$, $\lambda_3 = S_{ij}S_{jk}S_{ki}$, ...
• $y/\delta$ = normalized wall distance

Architecture: 6 → 32 → 32 → 1 (ReLU activations)

Tensor Basis Neural Network (Ling et al. 2016) for anisotropic Reynolds stresses:

$$b_{ij} = \sum_{n=1}^{10} g_n(\lambda_1, \ldots, \lambda_5) , T_{ij}^{(n)}$$

Architecture: 5 → 64 → 64 → 64 → 10 (outputs one coefficient per basis tensor)

Key Properties:

• Frame invariance: Guaranteed by using invariant inputs + tensor basis
• Realizability: Enforced during training
• Anisotropy: Captures different normal stresses and off-diagonal components

Pressure-driven flow between parallel plates.

Analytical solution (Poiseuille): $$u(y) = -\frac{1}{2\nu}\frac{dp}{dx}(H^2/4 - y^2)$$

Pressure-driven flow in square cross-section.

Classic benchmark for unsteady flow and energy decay.

Initial condition: $$u = \sin(x)\cos(y)\cos(z), \quad v = -\cos(x)\sin(y)\cos(z), \quad w = 0$$

Energy decay (low Re): $$KE(t) = KE(0) \cdot e^{-2\nu t}$$

All parameters can be set via command-line arguments (--param value) or config file (key-value pairs). Command-line arguments override config file values.

The solver uses the relationship: $\text{Re} = \frac{-dp/dx \cdot \delta^3}{3\nu^2}$ where $\delta$ is the channel half-height.

You should specify only TWO of (Re, nu, dp_dx). The third is computed automatically:

If all three are specified, the solver checks consistency and errors if they don't match (within 1% tolerance).

When adaptive_dt is enabled (default), the time step is computed each iteration as: $$\Delta t = \text{CFL} \cdot \min\left(\frac{\Delta x}{|u|{\max}}, \frac{\Delta y}{|v|{\max}}, \frac{\Delta z}{|w|_{\max}}\right)$$

• Steady mode: Iterates until $|\mathbf{u}^{n+1} - \mathbf{u}^n|_\infty < \text{tol}$ or max_iter reached
• Unsteady mode: Runs exactly max_iter time steps (or until T_final)

Available turbulence models:

For NN models, you must specify either:

• --nn_preset NAME (loads from data/models/<NAME>/), or
• --weights DIR --scaling DIR (explicit paths)

Available presets: tbnn_channel_caseholdout, tbnn_phll_caseholdout, example_tbnn, example_scalar_nut

Poisson solver options:

Auto-selection priority: FFT → FFT2D → FFT1D → HYPRE → MG

Convergence criteria (any triggers exit):

• tol_rhs: ‖r‖/‖b‖ < ε (recommended for projection)
• tol_rel: ‖r‖/‖r0‖ < ε
• tol_abs: ‖r‖ < ε

The --benchmark flag configures the solver for performance timing with minimal overhead:

# Run benchmark with defaults (192^3 grid, 20 iterations)
./duct --benchmark
# Override grid size
./duct --benchmark --Nx 256 --Ny 256 --Nz 256
# Override iteration count
./duct --benchmark --max_steps 100

Benchmark mode sets these defaults (all can be overridden by subsequent flags):

Config files use simple key-value syntax:

# Comment lines start with #
Nx = 128
Ny = 256
Re = 5000
turb_model = sst
adaptive_dt = true

Load a config file with --config FILE. Command-line arguments override config file values.

All solver components support GPU offload via OpenMP target directives.

# NVIDIA GPUs (NVHPC compiler)
CC=nvc CXX=nvc++ cmake .. -DCMAKE_BUILD_TYPE=Release -DUSE_GPU_OFFLOAD=ON
# With HYPRE PFMG (fastest Poisson solver)
CC=nvc CXX=nvc++ cmake .. -DUSE_GPU_OFFLOAD=ON -DUSE_HYPRE=ON

• Momentum equation (convection, diffusion)
• Pressure Poisson solver (multigrid V-cycles or HYPRE PFMG)
• Turbulence transport equations (k, ω)
• EARSM tensor basis computations
• Neural network inference

On NVIDIA GPUs with NVHPC compiler, the multigrid V-cycle is captured as a CUDA Graph:

• Eliminates per-kernel launch overhead
• Single cudaGraphLaunch() replaces O(levels ×ばつ kernels) launches
• Automatically recaptured if boundary conditions change

Train custom turbulence models on DNS/LES data:

# Setup environment
python3 -m venv venv && source venv/bin/activate
pip install -r requirements.txt
# Download dataset (~500 MB)
bash scripts/download_mcconkey_data.sh
# Train TBNN model
python scripts/train_tbnn_mcconkey.py \
 --data_dir mcconkey_data \
 --case channel \
 --output data/models/tbnn_channel \
 --epochs 100
# Use in solver
./channel --model nn_tbnn --nn_preset tbnn_channel

• Chorin, A. J. "Numerical solution of the Navier-Stokes equations." Math. Comput. 22.104 (1968): 745-762
• Briggs, W. L., Henson, V. E., & McCormick, S. F. A Multigrid Tutorial, 2nd ed. SIAM, 2000

• Menter, F. R. "Two-equation eddy-viscosity turbulence models for engineering applications." AIAA J. 32.8 (1994): 1598-1605
• Wilcox, D. C. "Reassessment of the scale-determining equation for advanced turbulence models." AIAA J. 26.11 (1988): 1299-1310
• Wallin, S., & Johansson, A. V. "An explicit algebraic Reynolds stress model..." J. Fluid Mech. 403 (2000): 89-132
• Gatski, T. B., & Speziale, C. G. "On explicit algebraic stress models..." J. Fluid Mech. 254 (1993): 59-78
• Pope, S. B. "A more general effective-viscosity hypothesis." J. Fluid Mech. 72.2 (1975): 331-340

• Ling, J., Kurzawski, A., & Templeton, J. "Reynolds averaged turbulence modelling using deep neural networks with embedded invariance." J. Fluid Mech. 807 (2016): 155-166
• Weatheritt, J., & Sandberg, R. D. "A novel evolutionary algorithm applied to algebraic modifications of the RANS stress-strain relationship." J. Comput. Phys. 325 (2016): 22-37

• McConkey, R., et al. "A curated dataset for data-driven turbulence modelling." Scientific Data 8 (2021): 255

MIT License - see license file