v4.2.0: Numba JIT compilation for TKE closure

Numba JIT compilation for the Deardorff TKE closure, eliminating the dominant performance bottleneck.

Performance

• Profiling showed the TKE closure consumed 86% of total runtime in both 2D and 3D simulations
• The inner column functions (_compute_mixing_length and compute_tke_closure_column) contained Python-level for k loops over vertical levels
• After JIT compilation: ~200x speedup per column call (from ~1 ms to ~5 μs)

Changes

• Numba @njit on TKE closure column functions — _compute_mixing_length_jit() and _compute_tke_closure_column_jit() accept only NumPy arrays and scalar primitives; Python wrappers preserve the original API for backward compatibility
• config/cbl_3d.yaml — t_end extended from 3600s (1 hour) to 7200s (2 hours)
• requirements.txt — Added numba

Backward Compatibility

• All public APIs unchanged — callers do not need any modifications
• First call incurs a one-time JIT compilation overhead (~1-2s); subsequent calls run at compiled speed
• Numerical results identical to v4.1.0

See RELEASE_NOTES.md for full details.

v4.1.0

New Features

• 5th-order upwind advection (Wicker & Skamarock 2002) for horizontal fluxes in 2D and 3D, configurable via advection_scheme: "upwind5" — provides implicit numerical dissipation, eliminating the need for explicit horizontal diffusion
• Beljaars (1994) convective velocity scale in MOST surface layer — prevents flux collapse in free-convective (low wind) conditions

Bug Fixes

• Fix MOST surface flux collapse when mean wind is weak by incorporating M_eff = sqrt(M² + (1.2·w*)²)
• Fix TKE over-initialization: use tke_min cold-start instead of convective scaling profile that caused excessive initial dissipation
• Increase theta perturbation amplitude to 0.1 K (standard LES practice)
• Set theta_surface to 307 K for consistency with target heat flux (~0.24 K m/s)

Other

• Add minimal 2D config (cbl_2d_minimal.yaml) for isolated physics testing

miniPBL v4.0.0

Addition of advanced physics parameterizations to miniPBL: vertical grid stretching, Deardorff prognostic TKE subgrid closure, Monin-Obukhov similarity theory surface layer, Rayleigh sponge damping, and large-scale subsidence. All operators updated to support variable vertical spacing. Existing 1D, 2D, and 3D modes are fully backward-compatible.

New Features

Vertical Grid Stretching

• Optional geometric stretching of the vertical grid with configurable stretch_factor (ratio between successive cell thicknesses) and nz_uniform (number of uniform cells near the surface before stretching begins)
• New grid arrays dz_center[k] (cell thickness) and dz_face[k] (distance between cell centers) replace the scalar dz throughout all operators
• When stretch_factor = 1.0 (default), the grid is identical to the previous uniform spacing

Variable-dz Support in All Operators

• Diffusion: flux F[k] = -K[k] * (f[k] - f[k-1]) / dz_face[k]; tendency (F[k+1] - F[k]) / dz_center[k]
• Advection: flux divergence uses local dz_center[k] and dz_face[k]
• Pressure Poisson solver: tridiagonal coefficients a[k] = 1/(dz_c[k]*dz_f[k]), c[k] = 1/(dz_c[k]*dz_f[k+1]) for variable spacing
• Divergence and pressure gradient correction use position-dependent spacing
• Boundary conditions (top gradient, surface stress) use local cell thickness

Deardorff TKE Subgrid-Scale Closure

• Prognostic TKE equation: de/dt = Shear + Buoyancy - Dissipation + Diffusion
• Shear production: S = K_m * (|du/dz|^2 + |dv/dz|^2)
• Buoyancy production/destruction: B = -(g/theta_ref) * K_h * dtheta/dz
• Dissipation: eps = c_eps * e^(3/2) / l with c_eps = 0.19 + 0.51*l/Delta
• Stability-dependent mixing length: l = min(c_l * sqrt(e) / N, Delta) (Delta for unstable conditions)
• Eddy viscosity: K_m = c_m * l * sqrt(e); eddy diffusivity: K_h = (1 + 2*l/Delta) * K_m
• TKE advanced by RK3 alongside theta, u, v, w; clamped to configurable minimum
• Column-by-column computation for 2D and 3D fields
• Activated by setting scheme: "deardorff-tke" in the turbulence config

Monin-Obukhov Similarity Theory (MOST) Surface Layer

• Iterative solver for friction velocity (u_star) and temperature scale (theta_star) given wind speed and temperature at the first grid level
• Businger-Dyer stability functions for momentum (psi_m) and heat (psi_h) under both stable and unstable conditions
• Surface heat flux: w'theta'_sfc = -u_star * theta_star
• Surface momentum flux: tau = -u_star^2 * (u,v) / |V|
• Replaces both the prescribed surface heat flux and the no-slip surface stress approximation when activated
• Activated by setting surface_flux_scheme: "most" in the physics config

Rayleigh Sponge Damping Layer

• Rayleigh damping in the upper domain: d(phi)/dt += -alpha(z) * (phi - phi_ref)
• Damping profile: alpha(z) = alpha_max * sin^2(pi/2 * (z - z_sponge) / (Lz - z_sponge))
• Applied to theta (relaxed toward horizontal mean), u and v (horizontal mean), w and TKE (zero)
• Configurable sponge_fraction (default 0.25) and sponge_alpha_max (default 0, i.e. off)

Large-Scale Subsidence

• Prescribed subsidence profile: w_s(z) = -D * z where D is the divergence rate
• Theta tendency: d(theta)/dt += -w_s * d(theta)/dz
• Prevents unbounded boundary layer growth in long simulations
• Activated by setting subsidence_divergence > 0 in the physics config

Configuration Changes

GridConfig

• stretch_factor: float = 1.0 — geometric stretch ratio (1.0 = uniform)
• nz_uniform: int = 0 — uniform cells near surface before stretching

PhysicsConfig

• surface_flux_scheme: str = "prescribed" — "prescribed" or "most"
• theta_surface: float = 302.0 — surface temperature for MOST (K)
• subsidence_divergence: float = 0.0 — large-scale divergence rate (1/s)
• sponge_fraction: float = 0.25 — fraction of domain for sponge layer
• sponge_alpha_max: float = 0.01 — maximum Rayleigh damping rate (1/s)

TurbulenceConfig

• scheme now supports "k-profile" (existing) and "deardorff-tke" (new)
• tke_c_m: float = 0.1 — eddy viscosity coefficient
• tke_c_eps_base: float = 0.19 — base dissipation coefficient
• tke_c_l: float = 0.76 — mixing length coefficient
• tke_min: float = 1e-4 — minimum TKE floor (m^2/s^2)

Backward Compatibility

All new features default to off. Existing configurations produce identical results:

• stretch_factor = 1.0: uniform grid (identical to v3.0.0)
• scheme = "k-profile": diagnostic K-profile closure (unchanged)
• surface_flux_scheme = "prescribed": fixed surface heat flux (unchanged)
• sponge_alpha_max = 0.0: no sponge damping
• subsidence_divergence = 0.0: no subsidence

Files Added

• minipbl/tke_closure.py — Deardorff prognostic TKE closure
• minipbl/surface_layer.py — Monin-Obukhov similarity theory surface fluxes
• minipbl/forcing.py — Rayleigh sponge damping and large-scale subsidence

Dependencies

• numpy, scipy, matplotlib, pyyaml

miniPBL v3.0.0

Extension of miniPBL from 2D x-z to full 3D x-y-z Boussinesq solver with a prognostic v velocity, Coriolis coupling on both u and v, and 3D pressure projection. The 1D and 2D modes are fully backward-compatible.

New Features

3D x-y-z Boussinesq Solver

• Prognostic v velocity on Arakawa C-grid y-faces, alongside existing theta, u, and w
• All arrays stored internally as (nx, ny, nz) with NetCDF output transposed to (time, z, y, x) for VisIT
• Automatic 3D mode activation when both nx > 1 and ny > 1 in the configuration

3D Pressure Solver

• PoissonSolver3D: 2D FFT in x and y (both periodic) + tridiagonal solve in z per wavenumber pair
• Combined eigenvalues lambda_xy[m,n] = lambda_x[m] + lambda_y[n] for each (kx, ky) mode
• kx=0, ky=0 mode: pressure pinned to zero (null-space treatment)
• project_velocity_3d: 3D divergence, pressure correction for u, v, and w

3D Advection

• Flux-form centered advection for theta, u, v, and w with y-flux terms
• Periodic y boundary conditions via np.roll on axis=1

3D Diffusion

• Horizontal Laplacian in x and y for cell-center and z-face fields
• Vertical + horizontal diffusion for theta, u, v, and w
• No-slip surface stress for both u and v

Coriolis Forcing

• Full Coriolis coupling: du/dt += f*(v - v_geo), dv/dt += -f*(u - u_geo)
• v interpolated to x-faces for u tendency; u interpolated to y-faces for v tendency

3D Turbulence

• Column-by-column K-profile closure over (nx, ny) columns
• Boundary layer height diagnosed per column: bl_height(nx, ny)

3D Boundary Conditions

• Rigid-lid: w[:,:,0] = 0, w[:,:,-1] = 0
• Fixed lapse rate enforcement at domain top for 3D theta fields

3D Output and Diagnostics

• NetCDF dimensions: (time, z_center, y_center, x_center) and (time, z_face, y_center, x_center)
• Variables: theta, u, v, w, p, heat_flux, K_h, bl_height
• Diagnostic plots (x-z cross-sections, x-y horizontal slices, xy-averaged profiles, BL height time series)

Configuration

• New grid parameters: ny, Ly
• New config file: config/cbl_3d.yaml (64x64x64, dt=0.5s, 1 hour simulation)

Backward Compatibility

• ny=1 (or omitting ny): dimensionality stays 1D or 2D; all existing code paths unchanged
• nx > 1, ny > 1: activates 3D code paths

Dependencies

• numpy, scipy, matplotlib, pyyaml

miniPBL v2.0.0

Extension of miniPBL from a 1D vertical column solver to a 2D x-z Boussinesq solver with resolved convection, momentum equations, and pressure projection. The 1D mode is fully backward-compatible with v1.0.0.

New Features

2D x-z Boussinesq Solver

• Prognostic equations for horizontal velocity (u), vertical velocity (w), and potential temperature (theta) on an Arakawa C-grid
• Periodic boundary conditions in x, rigid-lid (w=0) at top and bottom
• Automatic 2D mode activation when nx > 1 in the configuration

Pressure Solver

• FFT-based Poisson solver for enforcing incompressibility (divergence-free velocity)
• FFT in x (periodic) with tridiagonal solve in z per wavenumber
• Neumann boundary conditions in z; mean pressure pinned to zero for the kx=0 mode
• Fractional-step pressure projection applied at each RK3 sub-stage

Advection

• Second-order flux-form centered advection for theta, u, and w
• Proper C-grid staggering with periodic x boundary conditions

Momentum Physics

• Buoyancy forcing on vertical velocity from potential temperature perturbations
• Coriolis forcing with prescribed geostrophic wind
• No-slip surface stress parameterization for u
• Vertical turbulent momentum diffusion using K_m = K_m_ratio * K_h

Horizontal Diffusion

• Configurable horizontal diffusivity (K_horizontal) applied to theta, u, and w

Initialization

• Logarithmic wind profile initialization with configurable surface roughness length (z0)
• Small random theta perturbation seeded near the surface to trigger convective instability

Turbulence (2D)

• Column-by-column K-profile closure reusing the existing 1D scheme
• Separate K_h (heat) and K_m (momentum) eddy diffusivities with configurable ratio

Output and Diagnostics

• 2D NetCDF output with x_center dimension and u, w, p variables
• x-z cross-section plots (pcolormesh) of theta, u, w at selected times
• x-averaged theta profile plots for comparison with 1D results
• Progress log file (output/progress.log) for monitoring long simulations

Configuration

• New parameters: nx, Lx, g, coriolis_f, geostrophic_u, geostrophic_v, z0, K_m_ratio, K_horizontal
• New 2D config file: config/cbl_2d.yaml

Backward Compatibility

The 1D solver path is completely unchanged. Running with nx=1 (or omitting nx from the config) produces identical results to v1.0.0.

miniPBL v1.0.0

miniPBL v1.0.0

Initial release of miniPBL, a lightweight Python solver for planetary boundary layer simulation.

Features

Physics

• 1D vertical column model for the convective boundary layer (CBL)
• Prognostic potential temperature (θ) equation driven by vertical turbulent heat flux divergence
• K-profile turbulence closure with convective velocity scaling
• Automatic boundary layer height diagnosis with linear interpolation
• Prescribed surface kinematic heat flux (lower boundary)
• Zero-flux insulating lid (upper boundary)
• Optional fixed lapse rate enforcement at the domain top

Numerics

• Staggered vertical grid (scalars at cell centers, fluxes at cell faces)
• Second-order centered finite differences for diffusion
• Third-order Runge-Kutta (RK3) time integration (Wicker & Skamarock 2002)

Configuration

• YAML-based configuration for grid, physics, turbulence, time stepping, and output settings
• Input validation with descriptive error messages

Output

• NetCDF output of θ, heat flux, eddy diffusivity, and boundary layer height at configurable intervals
• Automatic diagnostic plots: θ profiles, BL height time series, and heat flux profiles

Dependencies

• numpy
• scipy
• matplotlib
• pyyaml