Bit-4-bit reproducability for AMReX C++ and Fortran code

[johnmauff]

Recently, when I compared the ocean.stats file for the double-gyre test case, I discovered that there is not a bit-4-bit agreement between the AMReX versions of the code and the Fortran. After some investigation, I narrowed it down to the following lines of code from the routine PPM_limit_pos from MOM_continuity_PPM.F90

      ...
    elseif (12.0*curv*(h_in(i,j,k) - h_min) < (curv**2 + 3.0*dh**2)) then
         scale = 12.0*curv*(h_in(i,j,k) - h_min) / (curv**2 + 3.0*dh**2)

The AMReX equivalent is

       else if (12.0_rt * curv * (h_in - h_min) < (curv * curv + 3.0_rt * dh * dh)) {
            Real const scale = 12.0_rt * curv * (h_in - h_min) / (curv * curv + 3.0_rt * dh * dh);

I was using the GCC compiler with a '-02' flag, which, as I learned, allows FMA. So I deactivated FMA by using the flag:
'-O2 -ffp-contract=off'. Strangely, Fortran and C++ still did not match. So I decided to be more specific about the order of operation, and rewrote the Fortran and C++ to be

     ((curv*curv) + 3.0 *(dh*dh))

Finally, with FMAs deactivated, I got the Fortran and C++ to match. However, when I turned FMA back on, neither Fortran nor C++ matched. Furthermore, even the original Fortran with the **2 operator and the modified Fortran did not match. So this creates a couple of problems/questions:

(1). Currently, we do not have a way to reproduce the original Fortran results in C++. Is this a problem, or does this fall into the category that different compilers will have unavoidable bit-rounding differences?

(2). I expect that the differences in floating point results for Fortran and C++ are well known to the E3SM group that did the Kokkos port of CAM. Can we leverage their knowledge base to match Fortran and C++

(3). Do we potentially sacrifice performance to match Fortran results, or should we just use some compatibility mode to simplify the verification of new AMReX functions?

[marshallward]

@johnmauff Have you tried (curv*curv) + (3.0 *(dh*dh))? The expression above will control the order of multiplication in each term, but won't control whether an FMA is used.

There is probably merit in replacing a**2 with a*a. I have observed Intel replacing a**3 with pow(a,3), which does not match a*a*a. I have never seen **2 substitution but a little paranoia might be helpful.

[johnmauff]

@marshallward, thanks for the suggestion. That indeed allows the C++ code to match that original Fortran with

curv**2 + 3.0 * dh**2