AD: fix OLAF treecode near-core regularization floor and TwrInfl OpenMP data race

[luwang00]

Feature or improvement description
FVW treecode: the wake path used Tree%DistanceDirect = 2*mean(RegParam), averaged over an over-allocated array whose sentinel tail drove it negative during wake buildup, disabling the near-core direct-evaluation fallback. Replace it with a per-node maxRegParam (max eps over the cell's particles); each branch now uses distDirect = BranchFactor*radius + 2*maxRegParam so control points inside a regularization core fall back to direct evaluation. Applied to both the particle and segment trees.

TwrInflArray: make FirstWarn_TowerStrike firstprivate and ErrStat2/ErrMsg2 private in the OpenMP loop to avoid a data race.

In ui_part_nograd_11, the idRegExp mollifier (1-exp(-(r/rc)^3))/r^3 differs from 1/r^3 by <exp(-8)~3.4e-4 once r>2*rc, which is exactly the accuracy the far-field multipole already accepts at its BranchFactor*radius+2*maxRegParam floor. Treat the mollifier as 1 beyond that boundary (new PART_REG_NRAD/PART_REG_CUT3 params) so the near-field direct kernel and the far-field multipole share the same 2*rc cutoff. Since most near-field tree pairs have r>>rc, this skips exp() for the majority of evaluations, giving ~25-32 percent serial speedup on the treecode path with the output unchanged to ~1.8e-5 relative. Also caches r^2/r^3/rc^3 to drop redundant calculations; the compact-support branch is refactored identically (bit-for-bit).

Potential impact
This bug fix can potentially generate better behaved rotor wake. Below are two snapshots of OLAF wake generated from the same AeroDyn driver/OLAF case provided by @RBergua. The simulation uses steady wind and all free near-wake panels. There are no frozen wake or far-wake panels. The two images show the same time instant. The only difference is the regularization fix.

Steady wind; without regularization fix

Steady wind; with regularization fix

Notice that the wake is significantly more ordered with the regularization fix. This is expected with the large CoreSpreadEddyVisc=1000 and age-based core spreading. At the same time, the simulation with regularization fix is considerably slower; however, this is expected because the large core radii, especially with the older wake elements further downstream, force the code to resort to the costly direct and exact evaluation of induction velocity more often to properly account for core regularization. Previously, the incorrect regularization floor effectively disables near-core direct-evaluation fallback in many cases, resulting in faster simulation but highly disordered wake. To speed up the simulation with the regularization fix, we can reduce CoreSpreadEddyVisc as appropriate.

Also, note that the two simulations provide very similar rotor loads. However, the core-regularization fix can be important if the far wake is of interest of if simulating multiple tandem rotors.

Two snapshots of the wake with turbulent wind are shown below. The turbulent inflow makes the wake far less organized. Nevertheless, the wake with regularization fix shows more structure than the one without.

Turbulent wind; without regularization fix

Turbulent wind; with regularization fix

The instantaneous wake streamwise velocity computed with the direct particle method and the particle tree method is compared below. The agreement is fairly good.

Related issue, if one exists
There are several reported issues with OLAF. Unclear if this fix addresses any of them.

Impacted areas of the software
OLAF

Generative AI usage
Co-authored-by: Microsoft Copilot copilot@microsoft.com
Co-authored-by: Anthropic Claude claude@anthropic.com

Test results, if applicable
One r-test needs rebasing: ad_BAR_OLAF. Only two output channels did not pass the threshold, but the changes are small, barely exceeding the threshold.

• r-test branch merging required

[Copilot]

Pull request overview

This PR adjusts OLAF/FVW induced-velocity computation to make near-core/direct-evaluation behavior robust during wake buildup, removes an OpenMP data race in tower-influence calculations, and optimizes the particle exponential-regularization kernel to reduce expensive exp() calls while keeping results within the stated tolerance.

Changes:

• Replace the global tree “direct-eval distance” floor (based on a mean regularization parameter) with a per-node maxRegParam floor used during traversal for both particle and segment trees.
• Fix an OpenMP data race in TwrInflArray by making loop-local state (FirstWarn_TowerStrike, ErrStat2/ErrMsg2) thread-safe.
• Optimize ui_part_nograd_11 by reusing r^2/r^3/rc^3 and skipping the exponential mollifier computation for r > 2*rc.

Reviewed changes

Copilot reviewed 4 out of 4 changed files in this pull request and generated no comments.

File	Description
modules/aerodyn/src/FVW_VortexTools.f90	Adds per-node maxRegParam and applies it in tree traversal to enforce a near-core direct-evaluation floor.
modules/aerodyn/src/FVW_Subs.f90	Restricts DEV NaN/sentinel checks to active particle ranges to avoid false negatives from preallocated tails.
modules/aerodyn/src/FVW_BiotSavart.f90	Speeds up particle exponential regularization by sharing the same 2*rc cutoff as the far-field floor and caching powers of r/rc.
modules/aerodyn/src/AeroDyn.f90	Removes an OpenMP data race in TwrInflArray via firstprivate/private scoping.

---

💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.

[RBergua]

To better understand the implications of this modification, I took the OLAF model from OC7 Phase III WP 3.1 and ran it with the original code and the proposed code.

The proposed modifications do not change the rotor loading:

However, there is a significant impact on the wake behavior. The OLAF model tested uses 3,400 near-wake panels (all of them free). No far-wake is included in the model. Before the fix, the CoreSpreadEddyVisc played a limited role. However, after the fix, the wake behavior is sensitive to it. Below I show the wake behavior in a horizontal plane at an elevation equal to the hub height.

Before the fix (CoreSpreadEddyVisc = 1000):

Affter the fix (CoreSpreadEddyVisc = 500):

As can be observed, the behavior is very different. Before the proposed changes, we can observe highly unstable flow in the wake (despite using CoreSpreadEddyVisc = 1000).

This can also be observed in the vertical plane at different downstream distances. The wake deficit shown below corresponds to an averaged timeframe of 10 s (from 20 to 30 s).

Before the fix:

After the fix:

It's also important to note that after the proposed fix, the computational time increases for higher CoreSpreadEddyVisc values.

Finally, this modification does not improve the issue reported here: #3429

[luwang00]

@RBergua Thanks for sharing the results! Regarding the increased computing time, I'm looking at implementing an alternative particle regularization function with true compact support. Right now, for particle/particle tree, OLAF either uses no regularization or exponential regularization. The latter is quite expensive. There is already a particle regularization function called idRegCompact in the code, but it is neither truly compact nor currently selectable. The plan is to replace it with a genuinely compact-support function that avoids transcendental calls like exp and allows the treecode near-core cutoff to be tightened. Hopefully, this will increase the speed by a nontrivial amount. A new RegFunctionPart input will be added to OLAF to let the user select particle regularization independently of the segment regularization controlled by RegFunction. RegFunctionPart will be inserted right after RegFunction. The options are 1 for exponential and 2 for compact. Default is 1. In case you want to test it, here is the work-in-progress: https://github.com/luwang00/openfast/tree/f/OLAF (PR #3457)

[ebranlard]

A good check is to compare the wake shapes between VelocityMethod =2 and 3 (particles only). The "best reference should be VelocityMethod=1 (segment N^2), but there is no one-to-one equivalence between the regularization of segments and particles. A rough overview of the wake shape for those three cases would be informative.

[ebranlard]

I agree with the changes, the mean(eps) previously introduced was indeed a temporary hack introduced as a compromise.
If we want to match the N^2 computation, we do need to include the effects within the 2 Reg Param region of each particle. But with this new implementation, for large epsilon of a particle (will occur in the far wake, especially with large CoreSpreadEddyVisc) computations will get quite expensive.

As a quick workaround, we should probably introduce a maximum reg param (either as part of the core model, or part of the tree model).

[luwang00]

Hi @ebranlard, thanks for reviewing the changes.

Yes, we can remove Tree%DistanceDirect and Tree_Seg%DistanceDirect since these will always be zero now, but maybe we should retain the local argument/variable DistanceDirect of ui_tree_part and ui_tree_segment? I think we can still get nonzero values assigned in LiftingLineInducedVelocities.

Regarding the increased computing time, it looks like the compact particle regularization (PR #3457) is helping somewhat. The same run that took >5 days is on track to complete in just over 3 days (won't get exact time until it's done). Max regularization param is a good idea. I think it should be part of the core model. Enforcing that in the tree code could cause discontinuous velocity field, and we are back to the same problem. Also, it would be better to have a dynamic bound based on minimum core overlap ratio to ensure wake smoothness.