Performance audit: inefficient compute patterns in the main time loop

# Performance audit: inefficient compute patterns in the main time loop

A code audit of the per-timestep hot path (advection → subgrid → wall functions/IBM → Poisson → integration → halos/BCs → thermodynamics) found a number of systemic performance issues. The recurring theme is **memory-bandwidth waste**: the model makes an estimated 15–25+ full-grid passes per RK3 substep on bookkeeping (zeroing, copying, snapshot/diffing) in addition to the physics, and the build configuration actively prevents SIMD vectorization.

All line references are permalinks to 21647cc. Findings marked ✅ were verified by hand in the source; the rest come from a systematic file-by-file audit.

---

## Tier 0 — build configuration (free performance, zero code risk)

- [ ] ✅ **Release builds trap floating-point exceptions.** `-ffpe-trap=invalid,zero,overflow` (GNU), `-fpe0` (Intel) and `-K trap=...` (Cray) sit in the *global* flags, so they apply to `-O3` builds:

  https://github.com/uDALES/u-dales/blob/21647cc6317a08276a39cf7dee124fb563b064bc/CMakeLists.txt#L42-L45

  FPE trapping forbids the compiler from speculating/masking FP operations, which blocks vectorization across the whole code. These flags should move to the Debug configuration only.
- [ ] **No architecture targeting.** Release is bare `-O3`; GNU then targets generic SSE2. Add `-march=native` (or a cluster-appropriate `-x`/`-xHost` flag for Intel) to enable AVX2/AVX-512, and consider `-fno-math-errno` so gfortran can vectorize `sqrt`/`log`.
- [ ] **No OpenMP.** The code is pure MPI; a hybrid mode would help at high core counts where the Poisson transposes become latency-bound. (Long-term item.)

## Tier 1 — hot-loop defects (every RK3 substep, full grid)

### Wrong loop nesting (k or j innermost over `(i,j,k)` arrays) ✅

The innermost loop must iterate `i` (column-major); these nests stride a full i×j plane per innermost iteration, so essentially every access is a cache miss:

- [ ] Vertical advection in `advecc_2nd`, `advecu_2nd`, `advecv_2nd` (`j,i,k` order). Note `advecw_2nd` is already correct *and* fused with the horizontal pass — it is the template for fixing the other three (which also do horizontal and vertical advection in two separate grid passes):

  https://github.com/uDALES/u-dales/blob/21647cc6317a08276a39cf7dee124fb563b064bc/src/advec_2nd.f90#L72-L87

- [ ] Pressure-gradient velocity correction in `tderive` (`i,j,k` order):

  https://github.com/uDALES/u-dales/blob/21647cc6317a08276a39cf7dee124fb563b064bc/src/modpois.f90#L1049-L1059

- [ ] The Smagorinsky closure runs `k,i,j` with `j` innermost across a ~36-element strain-rate stencil, and uses `** 2.` with a *real* exponent (may compile to `pow()`). The Vreman loop in the same file is correctly ordered:

  https://github.com/uDALES/u-dales/blob/21647cc6317a08276a39cf7dee124fb563b064bc/src/modsubgrid.f90#L217-L269

- [ ] `periodicEBcorr` (k innermost + 2 allreduces per call), `fixuinf1` (`j` innermost), `shiftedPBCs` (k innermost, plus a per-grid-point `sin()` of an i-only quantity) in `modforces.f90`.

### Poisson solver (`modpois.f90`) ✅

- [ ] Pencil work arrays allocated/deallocated **every call** (3× per step); should be persistent module arrays like `p`/`rhs` already are:

  https://github.com/uDALES/u-dales/blob/21647cc6317a08276a39cf7dee124fb563b064bc/src/modpois.f90#L442-L444

- [ ] FFTs execute **line-by-line 1D transforms** with a copy into a scratch vector and an unpack per pencil line; the y-direction gather `py(i,:,k)` is strided on top. FFTW's advanced (`plan_many`) interface removes the per-line call overhead and copies; the separate `px = px*fac` normalization passes can be folded into the unpack:

  https://github.com/uDALES/u-dales/blob/21647cc6317a08276a39cf7dee124fb563b064bc/src/modpois.f90#L479-L501

- [ ] Dead diagnostics paid for unconditionally every substep: `Fxy = pz` and `Fxyz = pz` full-field copies are consumed only by *commented-out* fielddump code; `rhs = p(...)` (L434) is similar; `fillps` recomputes the identical divergence a second time into `dpupdx/dpvpdy/dpwpdz`:

  https://github.com/uDALES/u-dales/blob/21647cc6317a08276a39cf7dee124fb563b064bc/src/modpois.f90#L550

- [ ] `solmpj` declares a full-3D automatic array `d` per call (also a stack hazard).

### IBM (`modibm.f90`) ✅

- [ ] `ibmwallfun` snapshots `up/vp/wp/thlp` into a freshly `allocate`d full-size array and diffs whole fields afterwards (~16 full-grid sweeps + a heap allocation per substep) — purely to recover stress/flux deltas that the wall functions apply only at sparse facet cells. Accumulating directly into `tau_x`/`thl_flux` inside the facet-section loop (the value is already in hand there) eliminates all of it:

  https://github.com/uDALES/u-dales/blob/21647cc6317a08276a39cf7dee124fb563b064bc/src/modibm.f90#L1191-L1231

- [ ] `bottom` does the same 8 full-grid sweeps even when `lbottom = .false.` (delta is nonzero only on the `kb` plane):

  https://github.com/uDALES/u-dales/blob/21647cc6317a08276a39cf7dee124fb563b064bc/src/modibm.f90#L2026-L2029

- [ ] Wall functions recompute `log`/`sqrt` of **time-constant geometry** (`log(dist/z0)` etc., 3–4 transcendentals per facet section per grid per substep — millions of calls per step) that could be precomputed per section in `initibmwallfun`.

### Time integration (`tstep.f90`) ✅

- [ ] `tstep_integrate` ends with 7 unconditional whole-array zero-fills (`qtp=0.` in dry runs, `e12p=0.` without the one-equation model, …) plus 7 whole-array copies (`um=u0` …) at substep 3 — ~14 full-field bookkeeping passes that could be fused into the integration nest or replaced by pointer swaps:

  https://github.com/uDALES/u-dales/blob/21647cc6317a08276a39cf7dee124fb563b064bc/src/tstep.f90#L317-L333

- [ ] The scalar update iterates `n` (4th, slowest index) innermost inside the `(k,j,i)` nest; hoist `do n` outside.
- [ ] `thl0c(...) = thl0(...)` is a separate full-field pass that could be fused into the `thl0` update loop.

### Kappa scheme (`advec_kappa.f90`)

- [ ] `advecc_kappa` zeroes two full-size automatic arrays three times each per call, updates `varp` over the whole halo region in three separate whole-array passes, and calls an `external` limiter function `rlim` once per face — which cannot be inlined, so the inner loops cannot vectorize. A single flux array with an inlined, branch-free (`merge`/`sign`) limiter is ~3–4× less memory traffic. `advecc_upw` similarly allocates its work array per call and makes 6 grid passes.

## Tier 2 — communication patterns (hurts at scale)

- [ ] **~18 + 2·nsv sequential halo exchanges per substep**, one variable at a time in `halos` (plus `ekm/ekh/pup/pvp/pwp/p/pres0` elsewhere). The `*m` fields (`um,vm,wm,thlm,qtm,svm`) are exchanged every substep although they change only once per full RK3 step — redundant 2 of 3 times:

  https://github.com/uDALES/u-dales/blob/21647cc6317a08276a39cf7dee124fb563b064bc/src/modboundary.f90#L79-L93

- [ ] **Scattered scalar collectives** that could be batched into single packed reductions: 2 allreduces for CFL/diffusion in `tstep_update`, up to 4 per substep in `masscorr` (one branch reduces a whole 3D field to get one scalar), 2 in `periodicEBcorr`, 2 slab-average collectives per substep with `BCtopm_pressure`, 7 facet-sized allreduces per step with `lwritefac` (ALLREDUCE where REDUCE to rank 0 suffices).
- [ ] Legacy `excjs`/`excis` allocate 4 buffers per call and post blocking ordered receives after the sends; persistent buffers + IRECV-first + WAITALL.
- [ ] The energy-balance facet solve runs entirely on rank 0 while all other ranks wait, followed by 5 separate broadcasts (`modEB.f90`); facets are independent and could be distributed. The dense longwave exchange recomputes `facem(m)*boltz*facT(m,1)**4` inside the O(nfcts²) loop instead of precomputing an O(nfcts) emission vector.

## Tier 3 — smaller / conditional

- [ ] `scalsource` sweeps the **entire grid** per source per scalar to deposit a 3σ Gaussian ball, with `sqrt`/`exp`/`erf` inside the mask — per-source index bounds computed once (as `vegetation.f90` already does; it is the model citizen here) make this O(source volume).
- [ ] ✅ Duplicated top-BC calls: `fluxtopscal`/`valuetopscal` are each called **twice**, but both already update `sv0` *and* `svm` internally — the second call is 100% duplicate work (pattern copied from the `fluxtop(qtm)`/`fluxtop(qt0)` pairs, where it is needed):

  https://github.com/uDALES/u-dales/blob/21647cc6317a08276a39cf7dee124fb563b064bc/src/modboundary.f90#L238-L248

- [ ] Vreman/Smagorinsky post-passes (`ekm = ekm + numol`, `ekh = ekm*prandtli`, `damp = max(damp,dampmin)`) are 3–5 extra whole-array sweeps (including halos) that fold into the closure loop; the `damp` array itself is dead weight (Piomelli damping is commented out, so it is written as 1.0 and multiplied everywhere). `sbshr/sbbuo/sbdiss` are stored to full 3D arrays then summed into `e12p` in a fourth pass.
- [ ] Dry runs still pay moist costs: `qtp` zeroing/updates, `qt0av`/`ql0av` slab averages (each a full-grid sweep + MPI collective) and the full `thv0` formula run even when `lmoist = .false.`.
- [ ] `nudge` operates on `(:,:,k)` sections including ghost cells; `fromztop` allocates its 1D work arrays per call; `checksim` makes 3 separate sweeps + 4 collectives where 1 sweep + 2 collectives suffice (though it is properly time-gated).

## Possible correctness bugs found in passing (may deserve separate issues)

- `wallfunheat`: `elseif (is_equal(norm, -yhat)) then flux = bctfxm` — almost certainly should be `bctfym` (https://github.com/uDALES/u-dales/blob/21647cc6317a08276a39cf7dee124fb563b064bc/src/modibm.f90#L1543)
- `modboundary.f90` `xm_periodic`/`ym_periodic`: the `loneeqn` blocks use loop variable `m` *after* the `do m = 1, ih` loop ends, so `e120/e12m` copy an out-of-range-by-one plane.
- `yqi_profile`: `do i = jb - 1, ie + 1` — lower bound should be `ib - 1` (https://github.com/uDALES/u-dales/blob/21647cc6317a08276a39cf7dee124fb563b064bc/src/modboundary.f90#L1071)
- `thermo` Newton–Raphson branch (`lqlnr`): iteration is uncapped and `niter` is used uninitialized.
- `EB`/`intqH`: `fachf/facef` are zeroed on every call including substeps 1–2, silently discarding substep-1/2 flux accumulations — worth checking whether that is intended.

## What is already good

Precomputed reciprocal metrics mean almost no divisions in inner loops; FFT plans and tridiagonal coefficients are built once at init; the batch Thomas solve is textbook-vectorized; the runtime IBM path uses sparse point lists rather than masked sweeps; `vegetation.f90` is exemplary (sparse lists, init-time precomputation, feature gating); and most loop nests are correctly ordered — the defects above are exceptions, not the rule.

## Suggested attack order

1. **Tier 0 CMake changes** (move FPE traps to Debug, add `-march=native`) — one-line edits, then rebenchmark; possibly a large fraction of the gap on its own.
2. Fix the inverted loop nests — mechanical, low-risk, golden-output-testable.
3. Kill the IBM snapshot/diff sweeps and the `tstep_integrate` zero/copy passes.
4. Poisson: persistent work arrays, delete dead diagnostic copies, batched FFTs via `plan_many`.
5. Halo aggregation + skip `*m` exchanges except on substep 1 — the biggest win at high rank counts.

Before the larger refactors, profile one representative urban case (gprof / VTune / `-pg` + a mid-size run) to confirm how time splits between Poisson, advection, subgrid and IBM.


	set(CMAKE_Fortran_FLAGS "-fdefault-real-8 -ffree-line-length-none -ffpe-trap=invalid,zero,overflow -std=f2008" CACHE STRING "Fortran flags" FORCE)
	set(CMAKE_Fortran_FLAGS_DEBUG "-g -fbacktrace -O0 -finit-real=nan -fcheck=all -Wall -Wextra -Wuninitialized -Warray-bounds -Wconversion" CACHE STRING "Fortran debug flags" FORCE)
	set(CMAKE_Fortran_FLAGS_RELEASE "-O3" CACHE STRING "Fortran release flags" FORCE)
	elseif(CMAKE_Fortran_COMPILER_ID MATCHES "Intel")

	do j = jb, je
	jm = j - 1
	jp = j + 1
	do i = ib, ie
	im = i - 1
	ip = i + 1
	do k = kb, ke
	km = k - 1
	kp = k + 1
	putout(i, j, k) = putout(i, j, k) - ( &
	w0(i, j, kp)(putin(i, j, kp)dzf(k) + putin(i, j, k)dzf(kp))dzhi(kp) &
	- w0(i, j, k)(putin(i, j, km)dzf(k) + putin(i, j, k)dzf(km))dzhi(k) &
	)*dzfi5(k)
	end do
	end do
	end do

	do i=ib,ie
	do j=jb,je
	up(i,j,kb) = up(i,j,kb)-(p(i,j,kb)-p(i-1,j,kb))*dxi
	vp(i,j,kb) = vp(i,j,kb)-(p(i,j,kb)-p(i,j-1,kb))*dyi
	do k=kb+1,ke
	up(i,j,k) = up(i,j,k)-(p(i,j,k)-p(i-1,j,k))*dxi
	vp(i,j,k) = vp(i,j,k)-(p(i,j,k)-p(i,j-1,k))*dyi
	wp(i,j,k) = wp(i,j,k)-(p(i,j,k)-p(i,j,k-1))*dzhi(k)
	end do
	end do
	end do

	do k = kb,ke
	kp=k+1
	km=k-1
	do i = ib,ie
	ip=i+1
	im=i-1
	mlen = csz(i,k) * delta(i,k)
	do j = jb,je
	jp=j+1
	jm=j-1

	! iw = wall(i,j,k,1) ! indices of closest wall
	! jw = wall(i,j,k,2)-myid*jmax ! indices of closest wall in local j-index
	! kw = wall(i,j,k,3)
	! c1 = wall(i,j,k,4) ! shear stress component
	! c2 = wall(i,j,k,5) ! shear stress component
	! if ((jw >= jb-1) .and. (jw <= je+1)) then ! check if jw is within the halo of this proc
	! !write(*,'(A,E9.2,A,E9.2,A,E9.2,A,E9.2)') 'component1:', c1, 'component2:', c2, 'shear c1:', shear(iw,jw,kw,c1), 'shear c2:', shear(iw,jw,kw,c2)
	! distplus = mindist(i,j,k)sqrt(abs(shear(iw,jw,kw,c1))+abs(shear(iw,jw,kw,c2)))numoli
	! damp(i,j,k) = sqrt(1. - exp((-distplus0.04)*3.)) ! Wall-damping according to Piomelli
	! ! write(*,'(A,2(1pE9.2))') 'damp, distplus', damp(i,j,k), distplus
	! else
	damp(i,j,k) = 1.
	! end if


	strain2 = ((u0(ip,j,k)-u0(i,j,k)) dxi )*2 &
	+ ((v0(i,jp,k)-v0(i,j,k)) dyi )*2 &
	+ ((w0(i,j,kp)-w0(i,j,k)) dzfi(k) )*2

	strain2 = strain2 + 0.125 * ( &
	+ ((w0(i ,j,kp)-w0(im,j ,kp))dxi + (u0(i ,j,kp)-u0(i ,j,k ))dzhi(kp))**2 &
	+ ((w0(i ,j,k )-w0(im,j ,k ))dxi + (u0(i ,j,k )-u0(i ,j,km))dzhi(k ))**2 &
	+ ((w0(ip,j,k )-w0(i ,j ,k ))dxi + (u0(ip,j,k )-u0(ip,j,km))dzhi(k ))**2 &
	+ ((w0(ip,j,kp)-w0(i ,j ,kp))dxi + (u0(ip,j,kp)-u0(ip,j,k ))dzhi(kp))**2 )

	strain2 = strain2 + 0.125 * ( &
	+ ((u0(i ,jp,k)-u0(i ,j ,k))dyi + (v0(i ,jp,k)-v0(im,jp,k))dxi)**2 &
	+ ((u0(i ,j ,k)-u0(i ,jm,k))dyi + (v0(i ,j ,k)-v0(im,j ,k))dxi)**2 &
	+ ((u0(ip,j ,k)-u0(ip,jm,k))dyi + (v0(ip,j ,k)-v0(i ,j ,k))dxi)**2 &
	+ ((u0(ip,jp,k)-u0(ip,j ,k))dyi + (v0(ip,jp,k)-v0(i ,jp,k))dxi)**2 )

	strain2 = strain2 + 0.125 * ( &
	+ ((v0(i,j ,kp)-v0(i,j ,k ))dzhi(kp) + (w0(i,j ,kp)-w0(i,jm,kp))dyi)**2 &
	+ ((v0(i,j ,k )-v0(i,j ,km))dzhi(k ) + (w0(i,j ,k )-w0(i,jm,k ))dyi)**2 &
	+ ((v0(i,jp,k )-v0(i,jp,km))dzhi(k ) + (w0(i,jp,k )-w0(i,j ,k ))dyi)**2 &
	+ ((v0(i,jp,kp)-v0(i,jp,k ))dzhi(kp) + (w0(i,jp,kp)-w0(i,j ,kp))dyi)**2 )

	ekm(i,j,k) = (mlendamp(i,j,k)) * 2. * sqrt(2. * strain2)
	ekh(i,j,k) = ekm(i,j,k) * prandtli
	end do
	end do
	end do

	do k=1,xsize(3)
	do j=1,xsize(2)
	Sxr = px(:,j,k)
	call dfftw_execute(plan_r2fc_x)
	px(1,j,k) = REAL(Sxfc(0))
	do i=1,itot/2-1
	px(2*i,j,k) = REAL(Sxfc(i))
	px(2*i+1,j,k) = AIMAG(Sxfc(i))
	end do
	px(itot,j,k) = REAL(Sxfc(itot/2))
	end do
	end do
	px = px*fac

	else
	fac = 1./sqrt(2.*itot)
	do k=1,xsize(3)
	do j=1,xsize(2)
	Sxr = px(:,j,k)
	call dfftw_execute(plan_r2fr_x)
	px(:,j,k) = Sxfr*fac
	end do
	end do

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Uh oh!

Performance audit: inefficient compute patterns in the main time loop #330

Performance audit: inefficient compute patterns in the main time loop

Tier 0 — build configuration (free performance, zero code risk)

Tier 1 — hot-loop defects (every RK3 substep, full grid)

Wrong loop nesting (k or j innermost over `(i,j,k)` arrays) ✅

Poisson solver (`modpois.f90`) ✅

IBM (`modibm.f90`) ✅

Time integration (`tstep.f90`) ✅

Kappa scheme (`advec_kappa.f90`)

Tier 2 — communication patterns (hurts at scale)

Tier 3 — smaller / conditional

Possible correctness bugs found in passing (may deserve separate issues)

What is already good

Suggested attack order

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

	call alloc_x(px, opt_xlevel=(/0,0,0/))
	call alloc_y(py, opt_ylevel=(/0,0,0/))
	call alloc_z(pz, opt_zlevel=(/0,0,0/))

	allocate(rhs(ib-ih:ie+ih,jb-jh:je+jh,kb:ke+kh))

	if (iwallmom > 1) then
	rhs = up
	call wallfunmom(xhat, up, bound_info_u)
	tau_x(:,:,kb:ke+kh) = tau_x(:,:,kb:ke+kh) + (up - rhs)

	rhs = vp
	call wallfunmom(yhat, vp, bound_info_v)
	tau_y(:,:,kb:ke+kh) = tau_y(:,:,kb:ke+kh) + (vp - rhs)

	rhs = wp
	call wallfunmom(zhat, wp, bound_info_w)
	tau_z(:,:,kb:ke+kh) = tau_z(:,:,kb:ke+kh) + (wp - rhs)

	! mom_flux_sum = sum(tau_x(ib:ie,jb:je,kb+1:ke) + tau_y(ib:ie,jb:je,kb+1:ke) + tau_z(ib:ie,jb:je,kb+1:ke))
	! call MPI_ALLREDUCE(mom_flux_sum, mom_flux_tot, 1, MY_REAL, MPI_SUM, comm3d, mpierr)
	! if (myid == 0) then
	! if (rk3step == 3) then
	! inquire(file="mom_flux.txt", exist=mom_flux_file_exists)
	! if (mom_flux_file_exists) then
	! open(12, file="mom_flux.txt", status="old", position="append", action="write")
	! else
	! open(12, file="mom_flux.txt", status="new", action="write")
	! end if
	! write(12, *) timee, -mom_flux_tot
	! close(12)
	! end if
	! end if
	end if

	call diffu_corr
	call diffv_corr
	call diffw_corr

	if (ltempeq .or. lmoist .or. lwritefac) then
	rhs = thlp
	totheatflux = 0 ! Reset total heat flux to zero so we only account for that in this step.
	totqflux = 0
	call wallfunheat
	thl_flux(:,:,kb:ke+kh) = thl_flux(:,:,kb:ke+kh) + (thlp - rhs)

	tau_x(:,:,kb:ke+kh) = up
	tau_y(:,:,kb:ke+kh) = vp
	tau_z(:,:,kb:ke+kh) = wp
	thl_flux(:,:,kb:ke+kh) = thlp

	up=0.
	vp=0.
	wp=0.
	thlp=0.
	svp=0.
	e12p=0.
	qtp=0.

	if(rk3step == 3) then
	um = u0
	vm = v0
	wm = w0
	thlm = thl0
	e12m = e120
	svm = sv0
	qtm = qt0
	end if

	call exchange_halo_z(u0)
	call exchange_halo_z(v0)
	call exchange_halo_z(w0)
	call exchange_halo_z(um)
	call exchange_halo_z(vm)
	call exchange_halo_z(wm)
	call exchange_halo_z(thl0)
	call exchange_halo_z(thlm)
	call exchange_halo_z(thl0c, opt_zlevel=(/ihc,jhc,khc/))
	call exchange_halo_z(qt0)
	call exchange_halo_z(qtm)
	do n = 1, nsv
	call exchange_halo_z(sv0(:, :, :, n), opt_zlevel=(/ihc,jhc,khc/))
	call exchange_halo_z(svm(:, :, :, n), opt_zlevel=(/ihc,jhc,khc/))
	enddo

	select case(BCtops)
	case(BCtops_flux)
	call fluxtopscal(wsvtop)
	call fluxtopscal(wsvtop)
	case(BCtops_value)
	call valuetopscal(sv_top)
	call valuetopscal(sv_top)
	case default
	write(0, *) "ERROR: top boundary type for scalars undefined"
	stop 1
	end select

Uh oh!

Performance audit: inefficient compute patterns in the main time loop #330

Description

Performance audit: inefficient compute patterns in the main time loop

Tier 0 — build configuration (free performance, zero code risk)

Tier 1 — hot-loop defects (every RK3 substep, full grid)

Wrong loop nesting (k or j innermost over (i,j,k) arrays) ✅

Poisson solver (modpois.f90) ✅

IBM (modibm.f90) ✅

Time integration (tstep.f90) ✅

Kappa scheme (advec_kappa.f90)

Tier 2 — communication patterns (hurts at scale)

Tier 3 — smaller / conditional

Possible correctness bugs found in passing (may deserve separate issues)

What is already good

Suggested attack order

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions

Wrong loop nesting (k or j innermost over `(i,j,k)` arrays) ✅

Poisson solver (`modpois.f90`) ✅

IBM (`modibm.f90`) ✅

Time integration (`tstep.f90`) ✅

Kappa scheme (`advec_kappa.f90`)