AeroDyn OLAF performance optimization

[luwang00]

Need to merge PR #3457 first.

Feature or improvement description
Speeds up the OLAF free-wake near-field particle-to-particle sum by

• restructuring the near-field interaction kernels so the compiler can vectorize them, and by
• tightening the treecode leaf handling.

Also added input validation for TreeBranchFactor requiring it to be >= 1 (default is 1.5) for multipole expansion convergence. (Note that the exact threshold for convergence depends on the particle distribution within a node with a worst case requirement of sqrt(3)=1.73, so taking >=1 as a reasonable requirement here.)

No new inputs; this is purely a data-flow/structure rework.

Based on one production case, the overall simulation time reduced by 41% from 2.8 days to 1.7 days.

Branch-free, SIMD-friendly near-field kernels

Rewrites the batch direct sum (ui_part_nograd) and adds a contained tree-leaf kernel (ui_leaves_nograd) in masked/branchless form with !$OMP SIMD reduction for the inner loop over source particles. (This is nested within the !$OMP PARALLEL outer loop over the query points.)

The singularity guard (cycle) is replaced by a mask (msk = merge(1,0,r²≥MINNORM²)) plus a floored divisor (r2s = max(r², MINNORM²)) so the divide is always finite (no Inf/NaN) and the loop has no data-dependent control flow. The RegFunction select case is hoisted out of the loops over particles so each inner loop is uniform. This lets the None and Compact branches vectorize; the Exponential branch is deliberately left scalar (its exp() blocks SIMD), so the Exponential kernel arithmetic itself is unchanged. (Overall Exponential kernel results still shift due to the leaf-multipole/MAC traversal change and tree branching early termination described below.)

Treecode leaf multipoles + Multipole Acceptance Criterion (MAC) on leaves + early tree branching termination

Childless nodes are now standard treecode leaves that carry multipole moments and are MAC-tested: a far query/control point uses the leaf's multipole instead of direct-summing the bucket, while a near query point direct-sums the leaf's particles via the SIMD leaf kernel.

Early termination kicks in when a node's radius drops below the largest core size in the cell (<=maxRegParam, i.e., a sub-core cell) or its particle count falls below the bucket cap (currently hardcoded to 32 as an SIMD friendly size). The first clause bounds tree depth in dense clumps; the second helps the inner loop vectorize.

These changes also simplify the tree traversal: childless nodes are no longer a force-direct special case but follow the same MAC path as internal nodes, and leaf buckets now carry multipole moments accepted when well-separated. This brings the implementation closer to a textbook treecode.

Note on terminology: This codebase uses "leaves" in a nonstandard way. It refers to the particles belonging to a node. In standard treecode terminology, leaves are the lowest-level (childless) nodes. This PR message uses the standard terminology: "leaf" = childless node, and "the leaf's particles/bucket" for the particles it holds.

Related issue, if one exists
Requires the compact polynomial core added through PR #3457.

Impacted areas of the software

FVW_BiotSavart.f90 masked/branchless SIMD rewrite of the batch ui_part_nograd direct sum.
FVW_VortexTools.f90 SIMD leaf kernel ui_leaves_nograd, leaf multipoles + MAC-on-leaves, early branching termination

No input-file, registry, or API changes.

Additional supporting information
The SIMD payoff requires a vectorizable regularization kernel (None or Compact). With the default Exponential kernel, the leaf loop stays scalar, so this PR yields little to no speedup. It pairs with the compact-C2 core from the earlier PR #3457, which removes the exp() and makes the near-field loop vectorizable. Considering the vast improvement in performance, we might want to make the compact-C2 core default in the future.

Numerical impact: SIMD reduction reordering changes results at floating-point truncation error level. The leaf-multipole/MAC change makes far childless-node interactions multipole-approximated rather than exact-direct, bounded by the existing BranchFactor acceptance criterion. Both are within regression tolerance, but not bit-for-bit. The early tree branching termination also alters the tree structure by reducing the depth of the tree.

With this optimization, we held short of M2M upward path for future work. M2M reduces multipole moment construction from O(NlogN) of direct/pure P2M to O(N), but it is not the bottleneck at the moment. Similarly, additional steps like M2L, L2L, and L2P for a full O(N) fast-multipole implementation are likely not helpful right now due to nearfield P2P dominance. These help only once far-field aggregation becomes a significant fraction (typically large N on the order of 10⁶). (P: Particle, M: Multipole expansion, L: Local expansion)

Test results, if applicable
Performance (multi-day production case with compact kernel): 41% wall-clock reduction (1.69×) from the SIMD near-field sum alone. Combining both SIMD vectorization from this PR and the Compact-C2 core from PR #3457, total speed up relative to the original OLAF with exponential cores is just over 3x on the same hardware.

Regression: matches within r-test tolerance.

Generative AI usage
Co-authored-by: Anthropic Claude claude@anthropic.com
Assisted by: GitHub Copilot support@github.com

[Copilot]

Copilot encountered an error and was unable to review this pull request. You can try again by re-requesting a review.

[Copilot]

🟡 Changes recommended

Unresolved critical numerical-correctness and input-compatibility issues must be addressed before approval.

Get a fresh assessment by requesting another Copilot review.

Review details

Suppressed comments (4)

modules/aerodyn/src/FVW_BiotSavart.f90:489

• The same r>2*rc replacement changes the scalar ui_part_nograd_11 kernel: for points just beyond the cutoff, it changes the factor from approximately 1-exp(-8) to 1. This routine is also the reference path used by callers and tests, so the PR's claimed exponential-kernel equivalence is not preserved. Keep the original exponential expression unless this numerical approximation is intentional and documented.

         if (r3 > PART_REG_CUT3*rc3) then ! r > 2*rc: mollifier -> 1 (skip exp), consistent with far-field multipole floor
            ScalarPart = r3_inv*fourpi_inv
         else
            E          = exp(-r3/rc3)
            ScalarPart = (1._ReKi-E)*r3_inv*fourpi_inv

modules/aerodyn/src/FVW_BiotSavart.f90:441

• The compact SIMD path computes tt = r2s/rc2 and the polynomial branch unconditionally before merge. FVW_IO rejects only negative regularization parameters, so a zero core can reach this code with rc2=0 and execute division by zero; the scalar ui_part_nograd_11 avoids that by selecting the outside-core branch first. Either reject zero cores for RegFunctionPart=2 or make this arithmetic safe for a zero core.

            tt  = r2s/rc2
            Cx = Alpha(2,ip)*dz - Alpha(3,ip)*dy
            Cy = Alpha(3,ip)*dx - Alpha(1,ip)*dz
            Cz = Alpha(1,ip)*dy - Alpha(2,ip)*dx
            sp = merge(fourpi_inv/r3, (35._ReKi + tt*(-42._ReKi + 15._ReKi*tt))*fourpi_inv/(8._ReKi*rc3), r2s >= rc2)

modules/aerodyn/src/FVW_VortexTools.f90:1696

• The leaf kernel repeats the same hard cutoff, so tree-accelerated direct leaf interactions also differ from the prior exponential kernel by about exp(-8) immediately beyond the boundary. This is a behavioral change beyond the stated SIMD/data-flow rewrite; either preserve the exponential expression or update the numerical-impact description and regression criteria.

               if (r3 > PART_REG_CUT3*rc3) then    ! r>2rc: mollifier->1, skip the expensive exp
                  sp = fourpi_inv/r3
               else
                  sp = (1.0_ReKi-exp(-r3/rc3))*fourpi_inv/r3
               end if

modules/aerodyn/src/FVW_VortexTools.f90:1710

• The leaf SIMD kernel has the same zero-core problem: for the accepted RegParam == 0 input, rc2 and rc3 are zero and tt = r2s/rc2 performs invalid arithmetic for every particle before MERGE selects the singular result. This should share the zero-core handling with ui_part_nograd so tree and direct particle methods remain valid under the same input.

               rc2 = (PART_REG_C2*Part%RegParam(ip))**2
               rc3 = rc2*PART_REG_C2*Part%RegParam(ip)
               tt  = r2s/rc2

• Files reviewed: 14/14 changed files
• Comments generated: 4
• Review effort level: Lite