From c4ee310d87623abe1a3e1d6fd2c43e1f3c9c9a30 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Wed, 20 Dec 2023 17:02:40 -0500 Subject: [PATCH] +Add cuberoot function 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. --- src/core/MOM_unit_tests.F90 | 13 +-- src/framework/MOM_intrinsic_functions.F90 | 104 +++++++++++++++++++++- 2 files changed, 111 insertions(+), 6 deletions(-) diff --git a/src/core/MOM_unit_tests.F90 b/src/core/MOM_unit_tests.F90 index 89383c4936..bd449d0b39 100644 --- a/src/core/MOM_unit_tests.F90 +++ b/src/core/MOM_unit_tests.F90 @@ -4,15 +4,16 @@ module MOM_unit_tests ! This file is part of MOM6. See LICENSE.md for the license. use MOM_error_handler, only : MOM_error, FATAL, is_root_pe - -use MOM_string_functions, only : string_functions_unit_tests -use MOM_remapping, only : remapping_unit_tests +use MOM_hor_bnd_diffusion, only : near_boundary_unit_tests +use MOM_intrinsic_functions, only : intrinsic_functions_unit_tests +use MOM_mixed_layer_restrat, only : mixedlayer_restrat_unit_tests use MOM_neutral_diffusion, only : neutral_diffusion_unit_tests use MOM_random, only : random_unit_tests -use MOM_hor_bnd_diffusion, only : near_boundary_unit_tests +use MOM_remapping, only : remapping_unit_tests +use MOM_string_functions, only : string_functions_unit_tests use MOM_CFC_cap, only : CFC_cap_unit_tests use MOM_EOS, only : EOS_unit_tests -use MOM_mixed_layer_restrat, only : mixedlayer_restrat_unit_tests + implicit none ; private public unit_tests @@ -36,6 +37,8 @@ subroutine unit_tests(verbosity) "MOM_unit_tests: EOS_unit_tests FAILED") if (remapping_unit_tests(verbose)) call MOM_error(FATAL, & "MOM_unit_tests: remapping_unit_tests FAILED") + if (intrinsic_functions_unit_tests(verbose)) call MOM_error(FATAL, & + "MOM_unit_tests: intrinsic_functions_unit_tests FAILED") if (neutral_diffusion_unit_tests(verbose)) call MOM_error(FATAL, & "MOM_unit_tests: neutralDiffusionUnitTests FAILED") if (random_unit_tests(verbose)) call MOM_error(FATAL, & diff --git a/src/framework/MOM_intrinsic_functions.F90 b/src/framework/MOM_intrinsic_functions.F90 index 2439c628fc..d189d3346a 100644 --- a/src/framework/MOM_intrinsic_functions.F90 +++ b/src/framework/MOM_intrinsic_functions.F90 @@ -4,9 +4,12 @@ module MOM_intrinsic_functions ! This file is part of MOM6. See LICENSE.md for the license. +use iso_fortran_env, only : stdout=>output_unit, stderr=>error_unit + implicit none ; private -public :: invcosh +public :: invcosh, cuberoot +public :: intrinsic_functions_unit_tests contains @@ -25,4 +28,103 @@ function invcosh(x) end function invcosh +!> Returns the cube root of a real argument at roundoff accuracy, in a form that works properly with +!! rescaling of the argument by integer powers of 8. If the argument is a NaN, a NaN is returned. +elemental function cuberoot(x) result(root) + real, intent(in) :: x !< The argument of cuberoot in arbitrary units cubed [A3] + real :: root !< The real cube root of x in arbitrary units [A] + + real :: asx ! The absolute value of x rescaled by an integer power of 8 to put it into + ! the range from 0.125 < a <= 1.0, in ambiguous units cubed [B3] + real :: num ! The numerator of an expression for the evolving estimate of the cube root of asx + ! in arbitrary units that can grow or shrink with each iteration [B C] + real :: den ! The denominator of an expression for the evolving estimate of the cube root of asx + ! in arbitrary units that can grow or shrink with each iteration [C] + real :: np2 ! The previous value of num squared [B2 C2] + integer :: ex_3 ! One third of the exponent part of x, used to rescale x to get a. + integer :: itt + + if ((x >= 0.0) .eqv. (x <= 0.0)) then + ! Return 0 for an input of 0, or NaN for a NaN input. + root = x + else + ex_3 = exponent(x) / 3 + asx = scale(abs(x), -3*ex_3) + if (asx > 1.0) then + asx = 0.125 * asx ; ex_3 = ex_3 + 1 + elseif (asx <= 0.125) then + asx = 8.0 * asx + endif + + num = (2.0 + asx) + den = 3.0 + if (asx < 1.0) then + ! Iteratively determine R = asx**1/3 using Newton's method, noting that in this case Newton's + ! method converges monotonically from above and needs no bounding. For the range of asx from + ! 0.125 to 1.0 with a first guess of 1.0, 6 iterations suffice to converge to roundoff. + do itt=1,6 + ! Newton's method iterates estimates as Root = Root - (Root**3 - asx) / (3.0 * Root**2). + ! Keeping the estimates in a fractional form Root = num / den allows this calculation with + ! fewer (or no) real divisions during the iterations before doing a single real division + ! at the end, and it is therefore more computationally efficienty. + + ! Because successive estimates of the numerator and denominator tend to be the cube of their + ! predecessors, the numerator and denominator need to be rescaled by division when they get + ! too large or small to avoid overflow or underflow. Were this being done for 32 bit reals, + ! the values to compare with would be about 1.0e12 and 1.0e-12. + if ((den > 1.0e100) .or. (den < 1.0e-100)) then + num = num / den ; den = 1.0 + endif + + ! Test whether this iteration would change anything. + if (abs(num**3 - asx * den**3) < 5.e-16 * num**3) exit + + np2 = num**2 + num = (2.0 * num**3 + asx * den**3) + den = 3.0 * (den * np2) + enddo + endif + root = sign(scale(num / den, ex_3), x) + endif + +end function cuberoot + +!> Returns true if any unit test of intrinsic_functions fails, or false if they all pass. +logical function intrinsic_functions_unit_tests(verbose) + logical, intent(in) :: verbose !< If true, write results to stdout + + ! Local variables + real :: testval ! A test value for self-consistency testing [nondim] + logical :: fail, v + + fail = .false. + v = verbose + write(stdout,*) '==== MOM_intrinsic_functions: intrinsic_functions_unit_tests ===' + + fail = fail .or. Test_cuberoot(v, 1.2345678901234e9) + fail = fail .or. Test_cuberoot(v, -9.8765432109876e-21) + fail = fail .or. Test_cuberoot(v, 64.0) + fail = fail .or. Test_cuberoot(v, -0.5000000000001) + fail = fail .or. Test_cuberoot(v, 0.0) + +end function intrinsic_functions_unit_tests + +!> True if the cube of cuberoot(val) does not closely match val. False otherwise. +logical function Test_cuberoot(verbose, val) + logical, intent(in) :: verbose !< If true, write results to stdout + real, intent(in) :: val !< The real value to test, in arbitrary units [A] + ! Local variables + real :: diff ! The difference between val and the cube root of its cube. + + diff = val - cuberoot(val**3) + Test_cuberoot = (abs(diff) > 2.0e-15*abs(val)) + + if (Test_cuberoot) then + write(stdout, '("For val = ",ES22.15,", (val - cuberoot(val**3))) = ",ES9.2," <-- FAIL")') val, diff + elseif (verbose) then + write(stdout, '("For val = ",ES22.15,", (val - cuberoot(val**3))) = ",ES9.2)') val, diff + + endif +end function Test_cuberoot + end module MOM_intrinsic_functions