Use a slightly better approximation of erfc() in the long-range truncation function

.gitignore

@@ -115,6 +115,7 @@ tests/testSingle
 tests/testTimeIntegration
 tests/testSPHStep
 tests/testExp
+tests/testErfc
 tests/testKernel
 tests/testKernelGrav
 tests/testKernelLongGrav

examples/EAGLE_low_z/EAGLE_6/eagle_6.yml

@@ -55,7 +55,7 @@ Statistics:
 # Parameters for the self-gravity scheme
 Gravity:
   eta:                      0.025    # Constant dimensionless multiplier for time integration.
-  MAC:                      geometric
+  MAC:                      adaptive
   epsilon_fmm:              0.001
   theta_cr:                 0.7      # Opening angle (Multipole acceptance criterion)
   use_tree_below_softening: 1

src/kernel_long_gravity.h

@@ -85,18 +85,28 @@ kernel_long_grav_derivatives(const float r, const float r_s_inv,
 
   const float exp_u2 = expf(-u2);
 
-  /* Compute erfcf(u) using eq. 7.1.25 of
+  /* Compute erfcf(u) using eq. 7.1.26 of
    * Abramowitz & Stegun, 1972.
    *
    * This has a *relative* error of less than 4e-3 over
-   * the range of interest (0 < u <  5) */
-
-  const float t = 1.f / (1.f + 0.47047f * u);
-
-  /* 0.3480242 * t - 0.0958798 * t^2 + 0.7478556 * t^3 */
-  float a = 0.7478556f;
-  a = a * t - 0.0958798f;
-  a = a * t + 0.3480242f;
+   * the range of interest (0 < u < 5)
+   *
+   * This is a good approximation to use since we already
+   * need exp(-u2) */
+
+  const float t = 1.f / (1.f + 0.3275911f * u);
+
+  const float a1 = 0.254829592f;
+  const float a2 = -0.284496736f;
+  const float a3 = 1.421413741f;
+  const float a4 = -1.453152027;
+  const float a5 = 1.061405429f;
+
+  /* a1 * t + a2 * t^2 + a3 * t^3 + a4 * t^4 + a5 * t^5 */
+  float a = a5 * t + a4;
+  a = a * t + a3;
+  a = a * t + a2;
+  a = a * t + a1;
   a = a * t;
 
   const float erfc_u = a * exp_u2;
@@ -199,18 +209,28 @@ kernel_long_grav_eval(const float r_over_r_s, float *restrict corr_f,
   const float u2 = u * u;
   const float exp_u2 = expf(-u2);
 
-  /* Compute erfcf(u) using eq. 7.1.25 of
+  /* Compute erfcf(u) using eq. 7.1.26 of
    * Abramowitz & Stegun, 1972.
    *
    * This has a *relative* error of less than 4e-3 over
-   * the range of interest (0 < u <  5) */
-
-  const float t = 1.f / (1.f + 0.47047f * u);
-
-  /* 0.3480242 * t - 0.0958798 * t^2 + 0.7478556 * t^3 */
-  float a = 0.7478556f;
-  a = a * t - 0.0958798f;
-  a = a * t + 0.3480242f;
+   * the range of interest (0 < u < 5)\
+   *
+   * This is a good approximation to use since we already
+   * need exp(-u2) */
+
+  const float t = 1.f / (1.f + 0.3275911f * u);
+
+  const float a1 = 0.254829592f;
+  const float a2 = -0.284496736f;
+  const float a3 = 1.421413741f;
+  const float a4 = -1.453152027;
+  const float a5 = 1.061405429f;
+
+  /* a1 * t + a2 * t^2 + a3 * t^3 + a4 * t^4 + a5 * t^5 */
+  float a = a5 * t + a4;
+  a = a * t + a3;
+  a = a * t + a2;
+  a = a * t + a1;
   a = a * t;
 
   const float erfc_u = a * exp_u2;

tests/Makefile.am

@@ -23,7 +23,7 @@ AM_LDFLAGS = ../src/.libs/libswiftsim.a $(HDF5_LDFLAGS) $(HDF5_LIBS) $(FFTW_LIBS
 TESTS = testGreetings testMaths testReading.sh testKernel testKernelLongGrav \
         testActivePair.sh test27cells.sh test27cellsPerturbed.sh testExp \
         testParser.sh test125cells.sh test125cellsPerturbed.sh testFFT \
-        testAdiabaticIndex testRandom testRandomSpacing \
+        testAdiabaticIndex testRandom testRandomSpacing testErfc \
         testMatrixInversion testThreadpool testDump testLogger testInteractions.sh \
         testVoronoi1D testVoronoi2D testVoronoi3D testGravityDerivatives \
 	testPeriodicBC.sh testPeriodicBCPerturbed.sh testPotentialSelf \
@@ -36,7 +36,7 @@ TESTS = testGreetings testMaths testReading.sh testKernel testKernelLongGrav \
 check_PROGRAMS = testGreetings testReading testTimeIntegration testKernelLongGrav \
 		 testActivePair test27cells test27cells_subset test125cells testParser \
                  testKernel testFFT testInteractions testMaths testRandom testExp \
-                 testSymmetry testThreadpool testRandomSpacing \
+                 testSymmetry testThreadpool testRandomSpacing testErfc \
                  testAdiabaticIndex testRiemannExact testRiemannTRRS \
                  testRiemannHLLC testMatrixInversion testDump testLogger \
 		 testVoronoi1D testVoronoi2D testVoronoi3D testPeriodicBC \
@@ -126,6 +126,8 @@ testLogger_SOURCES = testLogger.c
 
 testExp_SOURCES = testExp.c
 
+testErfc_SOURCES = testErfc.c
+
 testGravityDerivatives_SOURCES = testGravityDerivatives.c
 
 testGravitySpeed_SOURCES = testGravitySpeed.c

tests/testErfc.c

+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (C) 2020 Matthieu Schaller (schaller@strw.leidenuniv.nl)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+
+#include "../config.h"
+
+#include "swift.h"
+
+/* Standard includes */
+#include <fenv.h>
+#include <math.h>
+
+/**
+ * Compute erfcf(u) using eq. 7.1.26 of
+ * Abramowitz & Stegun, 1972.
+ *
+ * This has a *relative* error of less than 4e-3 over
+ * the range of interest (0 < x < 5)
+ */
+float optimized_erfcf(const float x) {
+
+  const float x2 = x * x;
+  const float exp_x2 = expf(-x2);
+
+  const float t = 1.f / (1.f + 0.3275911f * x);
+
+  const float a1 = 0.254829592f;
+  const float a2 = -0.284496736f;
+  const float a3 = 1.421413741f;
+  const float a4 = -1.453152027;
+  const float a5 = 1.061405429f;
+
+  /* a1 * t + a2 * t^2 + a3 * t^3 + a4 * t^4 + a5 * t^5 */
+  float a = a5 * t + a4;
+  a = a * t + a3;
+  a = a * t + a2;
+  a = a * t + a1;
+  a = a * t;
+
+  return a * exp_x2;
+}
+
+/**
+ * @brief Check that a and b are consistent (up to some relative error)
+ *
+ * @param a First value
+ * @param b Second value
+ * @param s String used to identify this check in messages
+ */
+void check_value(const double a, const double b, const double rel_tol,
+                 const double abs_tol, const double x) {
+
+  if (fabs(a - b) / fabs(a + b) > rel_tol)
+    error("Values are inconsistent: %12.15e %12.15e rel=%e (for x=%e).", a, b,
+          fabs(a - b) / fabs(a + b), x);
+  if (fabs(a - b) > abs_tol)
+    error("Values are inconsistent: %12.15e %12.15e abs=%e (for x=%e).", a, b,
+          fabs(a - b), x);
+}
+
+int main(int argc, char* argv[]) {
+
+  /* Initialize CPU frequency, this also starts time. */
+  unsigned long long cpufreq = 0;
+  clocks_set_cpufreq(cpufreq);
+
+/* Choke on FPEs */
+#ifdef HAVE_FE_ENABLE_EXCEPT
+  feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
+#endif
+
+  /* Loop over some values */
+  for (float x = 0.f; x < 5.f; x += 0.000001f) {
+
+    const double exact = erfc(x);
+    const double swift_erfcf = optimized_erfcf(x);
+
+    check_value(exact, swift_erfcf, 3.358e-3, 6e-7, x);
+  }
+
+  return 0;
+}