+Add cuberoot function

[Hallberg-NOAA]

Added the new functions cuberoot() and intrinsic_functions_unit_tests() to the MOM_intrinsic_functions module, and call intrinsic_functions_unit_tests from unit_tests to confirm that this new function works as intended. Separately, cuberoot was tested by replacing expressions like A**(1./3.) with cuberoot(A) in MOM_energetic_PBL and verifying that the answers only change at roundoff, but that it can give bitwise identical results when the argument is scaled by an integer power of 8 and then unscaled by the corresponding integer power of 2, but that change will occur in a subsequent commit as it can change answers depending on an ANSWER_DATE flag. With this commit, cuberoot is not yet being used so all answers are bitwise identical, although there are new publicly visible routines.

[codecov]

Codecov Report

Attention: 7 lines in your changes are missing coverage. Please review.

Comparison is base (5137442) 37.23% compared to head (655f0d5) 37.20%.

Files	Patch %	Lines
src/framework/MOM_intrinsic_functions.F90	88.88%	3 Missing and 2 partials ⚠️
src/core/MOM_unit_tests.F90	0.00%	1 Missing and 1 partial ⚠️

Additional details and impacted files

@@             Coverage Diff              @@
##           dev/gfdl     #538      +/-   ##
============================================
- Coverage     37.23%   37.20%   -0.04%     
============================================
  Files           271      271              
  Lines         80170    80355     +185     
  Branches      14978    14985       +7     
============================================
+ Hits          29855    29894      +39     
- Misses        44761    44903     +142     
- Partials       5554     5558       +4     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[marshallward]

I am seeing a small issue with the cuberoot function. If I evaluate cuberoot(0.125) then it gives me 1.0 rather than 0.5. Any other value that is a power of 0.125 (e.g. 2.**-6) shows the same problem.

I believe it is due to the elseif in this if-block:

    if (asx > 1.0) then
      asx = 0.125 * asx ; ex_3 = ex_3 + 1
    elseif (asx <= 0.125) then
      asx = 8.0 * asx
    endif

Values greater that 1.0 correct their exponent (ex_3 += 1) after scaling down, but values below 0.125 do not include this index correct after scaling up. It looks like the bug can be fixed with ex_3 = ex_3 - 1.

---

I looked at the behavior of the uncorrected asx function, scale(abs(x), -3 * (exponent(x) / 3)), and it appears to already be bounded between 0.125 and 4; only powers of 0.125 are evaluated to 0.125, so it's an extremely small number of values that would be affected by this missing term.

I considered removing the entire elseif block, but it seemed to introduce a minor roundoff that could not be removed by iteration. I don't know how important it is to reproduce the answer to machine precision, but something else to consider. (Are all values correct to machine precision?) (Edit: Nevermind, I had modified the 1e100/1e-100 convergence values. 0.125 converges to 0.5 without the elseif block.)

[Hallberg-NOAA]

Thank you for pointing out the bug with the (asx <= 0.125) case. I just fixed this in a force-push.

The Newton's method iterations do give solutions that are accurate to machine precision, but in some cases when we have arrived at the answer the last bit bounces back and forth around the exact (irrational) solution with added iterations. It is worth noting that the Newton's method convergence from above is monotonic (at least mathematically). If the first guess is an under-estimate, Newton's method will overshoot (a little) and then all subsequent iterations converge monotonically from above. So the case where the last bit is jumping around, it means that we are in a perpetual roundoff-induced two-step limit cycle.

Given this limiting behavior at effective convergence, we could choose to break the iterations when they are at the minimum of two successive solutions that are within roundoff of each other, but it would involve adding lines to store the previous value, and then to compare the solutions to make sure that they are not increasing, and this would add to the cost of the function without giving any meaningful increase in accuracy. Is it worth the (hopefully small) extra cost to avoid incorporating any arbitrary choices about the number of iterations to apply, even if it means that our version of the cube root would then always become an under-estimate of some irrational solutions?

[ashao]

What about checking the size of the Newton step as a stopping criterion if you're ok with sacrificing a little bit of extra accuracy in the solution? Something like:

if (ABS(update) < 2.0*10.0*EPSILON(root)) exit

(factor of 2.0 for the roundoff calculation and a factor of 10.0 to go one order of magnitude above roundoff).

[marshallward]

I have a suggestion to keep the value of asx bounded between 0.125 and 1. The current form has some odd behavior due to Fortran-isms related to integer division.

The current formula, asx = scale(abs(x), -3*(exponent(x) / 3)), looks something like this:

After which the if-blocks are required to bound the values above 1 to something within 0.125 and 1. (I did not look at the actual curve, but surely it will be some piecewise monstrosity.)

This is happening because Fortran integer division simply drops the decimal part of the number, rather than consistently rounding up or down. For example, 1.5 / 2 is zero, but -1.5 / 2 is also zero rather than -1. If we replace exponent(x) / 3 with ceiling(exponent(x) / 3.) and remove the if-blocks (which now do nothing), then we get a more consistent curve:

This new expression does not appear to change answers in your test.

[Hallberg-NOAA]

Thank you for all of this helpful feedback, @marshallward. I have incorporated all of this into the updated version of this PR, which is now algorithmically lean and thoroughly commented!

[Hallberg-NOAA]

@ashao, thank you for your suggestion with respect to checking the increment size, but in this case it turns out that we have another way to check for convergence and the extra iterations seem pretty cheap, so we can iterate to convergence and then take the better of two near-solutions if the iterations encounter a roundoff-level limit cycle.

[marshallward]

I've added what should be my final comments. It seems to work very well.

I have not checked whether it gives the same answers across compilers and hardware, but we should probably aspire to this at some point.

[marshallward]

Everything looks good!

[marshallward]

Gaea regression: https://gitlab.gfdl.noaa.gov/ogrp/MOM6/-/pipelines/21886 ✔️