Induction differences between single and double precision

[ebranlard]

Bug description
Large differences in inductions are obtained between the single and double precision versions (in debug or release).

the structural and wind velocity at the airfoil section are the same in both versions. I suspect something happening in the BEM implementation and in particular in the high-thrust correction (from some preliminary debugging done about a year ago).

To Reproduce
Compile the AeroDyn driver with single and double precision and run in on the test case enclosed.

AD_Sine_Bug.zip

(NOTE: this is a simplification of the AeroDyn driver test case BAR_SineMotion, so solving this bug might affect the r-test)

Expected behavior
The Axial and tangential induction should match between the two versions.

Screenshots, if applicable

OpenFAST Version
Tried on Windows and Mac, with OpenFAST 3.2, but likely several years old.

[jjonkman]

Hi @ebranlard,

For some reason I can't unzip your files to take a look, but just a quick question.

Presumably you have already ruled the influence of how the "DEFAULT" value of BEM input IndToler effects the results; is that correct? Per the AeroDyn documentation: https://openfast.readthedocs.io/en/main/source/user/aerodyn/input.html#blade-element-momentum-theory-options, BEM input IndToler is DEFAULTed to different values depending on whether you are compiling in single or double precision.

Thanks,

[ebranlard]

Checking IndToler was a good idea but unfortunately setting it to the same value for single and debug (5e-5) didn't narrow the differences.

[bjonkman]

It is possible that small differences in the numerics in the inductionFactors routine can cause the BEM algorithm to find a solution in a different region. We noticed this, particularly in the vicinity of phi=0. We changed the local variables in inductionFactors to double precision (R8Ki) to fix this issue.

Not sure if that's the same issue you are seeing, but thought I'd mention it.

[deslaughter]

The difference between the Single and Double precision results is caused by the BEM algorithm finding a solution for phi in different regions. This is due to the selection of the test ranges for phi based on the BEMT_epsilon2 variable defined in BEMTUncoupled.f90. The test ranges are calculated in GetSolveRegionOrdering in BEMT.f90. The following table compares the values of the relevant variables for Single and Double precision compilation:

Variable	Single Precision	Double Precision	Calculation
BEMT_epsilon2	0.00345267	1.490116e-07	10.0_ReKi*sqrt(epsilon(1.0_ReKi))
phi_lower(1)	0.00345267	1.490116e-07	BEMT_epsilon2
phi_upper(1)	1.56734371	1.570796178	PiBy2 - BEMT_epsilon2
phi_lower(2)	-0.7853982	-0.7853982	-pi/4.0_ReKi
phi_upper(2)	-0.00345267	-1.490116e-07	-BEMT_epsilon2

When phi_lower(1) and phi_upper(1) are passed to FindTestRegion, in the first iteration of the previously attached test case (AD_Sine_Bug.zip), the residual signs return by BEMTU_InductionWithResidual are the same for Single precision, but different for Double precision (opposite signs indicate the test region is valid). So, the Single precision algorithm checks the next range, phi_lower(2) and phi_upper(2) and finds this range is valid for calculating the induction. The following table summarizes the output of the algorithm after converging on phi and calculating the axial and tangential induction for the first iteration (these values match those shown in the plots from the original post, at time zero):

Variable	Single Precision	Double Precision
phi	-0.0176572613	1.603164e-05
tanInduction	0.00120669336	-0.27350477
axInduction	1.2943496	0.88434187

I'm curious about the selection of the phi test ranges and how both positive and negative phi can yield valid inductions for this case. The BEM solution method, including test range selection, is described in A simple solution method for the blade element momentum equations with guaranteed convergence, but I didn't see a discussion of negative phi values (I may have missed it).

As for a solution, reducing BEMT_epsilon2 does change the test region and sufficiently small values cause a solution to be found in the same region for Single and Double precision. A comment near the variable definition in the source code notes that the value must be large enough so that EqualRealNos(BEMT_epsilon2, 0.0_ReKi) is false. For Single precision, setting BEMT_epsilon2 greater than 100.0_SiKi*Eps / 2.0_SiKi will satisfy the comment and using a value of 100.0_ReKi*epsilon(1.0_ReKi) will select the correct test region in the example case. However, if the tolerance in EqualRealNos4 is changed in the future, this issue could reappear.

Should the algorithm be changed to handle this case differently or is changing the value of BEMT_epsilon2 sufficient?

Also, here is a plot comparing AB1002TnInd in Single and Double precision for the test input after changing BEMT_epsilon2 as described above:

*Edited to fix values for phi_lower(2) and phi_upper(2)

[ebranlard]

Thanks a lot @deslaughter for looking into it! I personally would need to dive a bit more into the algorithm to be able to see what should be preferred, but looking into it again, the "double version" is likely the wrong one as it leads to tangential induced velocity of ~10m/s (signal Vindy).

I just tested it, and it seems the AeroDyn code from Envision returns values closer to the "single version". The Envision version has a different ways to solve for the regions, and returns 4 regions instead of 3. Yet, I've tried using their sorting of regions, and that alone didn't do the trick. For a different project, I'm looking at a step by step approach to try to reduce the differences between the Envision and current AeroDyn version. I might be able to track down what changes was introduced that solved this issue.

But in the meantime, maybe we could find a quick fix so that the "double version" matches the "single version". I'll chat with you internally.

[jjonkman]

Thanks for the detailed summary. It sounds like you have a good path forward.

I'll just respond to @deslaughter's comment:

The BEM solution method, including test range selection, is described in A simple solution method for the blade element momentum equations with guaranteed convergence, but I didn't see a discussion of negative phi values (I may have missed it).

Negative values of phi (between -pi/4 and -eps) represent the propeller brake region in Ning's paper.

Thanks,

[bjonkman]

@deslaughter , I believe the 2013 paper you reference was later updated to include the negative phi values (they were omitted in early versions of this work -- including the code). I don't readily have the new reference, but I'll bet @jjonkman does. It would probably be dated 2017 or later.

[jjonkman]

Actually, the negative phi values and the propeller brake region are discussed in section 3.2 and is in Algorithm 1 of Ning's original BEM solution paper that @deslaughter referenced.

I never reviewed Ning's updated paper, but from a quick Google search, I believe this is the one @bjonkman is referring to: https://link.springer.com/content/pdf/10.1007/s00158-021-02883-6.pdf