Rounding errors on mac gfortran compilers

[conradtchan]

Running tests (e.g. polytrope) on the mac versions of gfortran produces different results compared to on Linux. Running on an Ubuntu container on mac gives results consistent with Linux.

These errors are typically of order 1.d-15, but when FMR is enabled on hydrostatic solutions, this can result in significant velocity errors due to the coordinate transformation required for vectors along the theta and phi directions.

[dliptai]

Relative errors after 1 timestep with smearing, polytrope problem: Git sha: b8e35d6cbb8d7b2775877d096afb72641ef08dfc (MPI=no, OMP_NUM_THREADS=1)

[dliptai]

Relative errors after setup+gravity relax (i.e. before any hydro fluxes)

[dliptai]

The following changes eliminate the roundoff error in the smeared region. (Converting from spherical to cartesian coordinates using quad precision for the trigonometric functions: sin and cos).

diff --git a/src/utils.f90 b/src/utils.f90
index e182d4e..2b8dca7 100644
--- a/src/utils.f90
+++ b/src/utils.f90
@@ -8,9 +8,9 @@ contains
   implicit none
   real(8),intent(in):: xp(1:3)
   real(8):: x(1:3)
-  x(1) = xp(1)*sin(xp(2))*cos(xp(3))
-  x(2) = xp(1)*sin(xp(2))*sin(xp(3))
-  x(3) = xp(1)*cos(xp(2))
+  x(1) = xp(1)*real(sin(real(xp(2),kind=16))*cos(real(xp(3),kind=16)),kind=8)
+  x(2) = xp(1)*real(sin(real(xp(2),kind=16))*sin(real(xp(3),kind=16)),kind=8)
+  x(3) = xp(1)*real(cos(real(xp(2),kind=16)),kind=8)
  end function polcar

 ! convert cartesian to polar coordinates

[dliptai]

Intel compilers (ifort/ifx) also differ. This is the error after setup between ifort+mac and gfortran+linux

[ryosuke-hirai]

After some investigation, I think this is not just about coordinate transformation but it's also simply a roundoff error. As a test, I changed the following line in subroutine angular_smear

momtot = momtot + vcar*dvol(i,j,k)! add up momenta

to

momtot = momtot + abs(vcar)*dvol(i,j,k)! add up momenta

When comparing the momtot values between the previous and new results, they are off by >16 orders of magnitude, meaning that basically the momenta are perfectly cancelling each other. So the results can look wildly different with even the smallest changes.