Replace abs() by fabs()

Replace calls to abs() with floating arguments by calls to fabs(), to ensure the results are not accidentally truncated to integer.

12 files changed, 44 insertions, 49 deletions

diff --git a/Carpet/CarpetLib/src/gdata.cc b/Carpet/CarpetLib/src/gdata.cc
index c26d374a1..ecc8c4ba8 100644
--- a/Carpet/CarpetLib/src/gdata.cc
+++ b/Carpet/CarpetLib/src/gdata.cc
@@ -6,6 +6,7 @@
 #include <vectors.h>
 
 #include <cassert>
+#include <cmath>
 #include <cstdlib>
 #include <iomanip>
 #include <iostream>
@@ -401,7 +402,7 @@ find_source_timelevel (vector <CCTK_REAL> const & times,
   CCTK_REAL const min_time = * min_element (times.begin(), times.end());
   CCTK_REAL const max_time = * max_element (times.begin(), times.end());
   // TODO: Use a real delta-time from somewhere instead of 1.0
-  CCTK_REAL const some_time = abs (min_time) + abs (max_time) + 1.0;
+  CCTK_REAL const some_time = fabs (min_time) + fabs (max_time) + 1.0;
   if (op != op_copy) {
     if (time < min_time - eps * some_time or
         time > max_time + eps * some_time)
@@ -436,9 +437,7 @@ find_source_timelevel (vector <CCTK_REAL> const & times,
   }
   if (timelevel == -1) {
     for (size_t tl=0; tl<times.size(); ++tl) {
-      static_assert (abs(0.1) > 0,
-                     "Function CarpetLib::abs has wrong signature");
-      if (abs (times.AT(tl) - time) < eps * some_time) {
+      if (fabs (times.AT(tl) - time) < eps * some_time) {
         timelevel = tl;
         break;
       }
diff --git a/Carpet/CarpetLib/src/ggf.cc b/Carpet/CarpetLib/src/ggf.cc
index 17f9d9406..0f5015897 100644
--- a/Carpet/CarpetLib/src/ggf.cc
+++ b/Carpet/CarpetLib/src/ggf.cc
@@ -401,9 +401,8 @@ mg_restrict_all (comm_state & state,
   static Timer timer ("mg_restrict_all");
   timer.start ();
   // Require same times
-  static_assert (abs(0.1) > 0, "Function CarpetLib::abs has wrong signature");
-  assert (abs(t.get_time(ml,rl,0) - t.get_time(ml-1,rl,0))
-	  <= 1.0e-8 * (1.0 + abs(t.get_time(ml,rl,0))));
+  assert (fabs(t.get_time(ml,rl,0) - t.get_time(ml-1,rl,0))
+	  <= 1.0e-8 * (1.0 + fabs(t.get_time(ml,rl,0))));
   vector<int> const tl2s(1,tl);
   transfer_from_all (state,
                      tl  ,rl,ml,
@@ -425,9 +424,8 @@ mg_prolongate_all (comm_state & state,
   static Timer timer ("mg_prolongate_all");
   timer.start ();
   // Require same times
-  static_assert (abs(0.1) > 0, "Function CarpetLib::abs has wrong signature");
-  assert (abs(t.get_time(ml,rl,0) - t.get_time(ml+1,rl,0))
-	  <= 1.0e-8 * (1.0 + abs(t.get_time(ml,rl,0))));
+  assert (fabs(t.get_time(ml,rl,0) - t.get_time(ml+1,rl,0))
+	  <= 1.0e-8 * (1.0 + fabs(t.get_time(ml,rl,0))));
   vector<int> const tl2s(1,tl);
   transfer_from_all (state,
                      tl  ,rl,ml,
@@ -449,9 +447,8 @@ ref_restrict_all (comm_state & state,
   static Timer timer ("ref_restrict_all");
   timer.start ();
   // Require same times
-  static_assert (abs(0.1) > 0, "Function CarpetLib::abs has wrong signature");
-  assert (abs(t.get_time(ml,rl,tl) - t.get_time(ml,rl+1,tl)) <=
-          1.0e-8 * (1.0 + abs(t.get_time(ml,rl,tl))));
+  assert (fabs(t.get_time(ml,rl,tl) - t.get_time(ml,rl+1,tl)) <=
+          1.0e-8 * (1.0 + fabs(t.get_time(ml,rl,tl))));
   transfer_from_all (state,
                      tl,rl  ,ml,
                      & dh::fast_dboxes::fast_ref_rest_sendrecv,
@@ -500,9 +497,8 @@ ref_reflux_all (comm_state & state,
   static Timer timer ("ref_reflux_all");
   timer.start ();
   // Require same times
-  static_assert (abs(0.1) > 0, "Function CarpetLib::abs has wrong signature");
-  assert (abs(t.get_time(ml,rl,tl) - t.get_time(ml,rl+1,tl)) <=
-          1.0e-8 * (1.0 + abs(t.get_time(ml,rl,tl))));
+  assert (fabs(t.get_time(ml,rl,tl) - t.get_time(ml,rl+1,tl)) <=
+          1.0e-8 * (1.0 + fabs(t.get_time(ml,rl,tl))));
   islab slabinfo;
   slabinfo.is_centered = 1 - ivect::dir(dir);
   transfer_from_all (state,
@@ -526,9 +522,8 @@ ref_reflux_prolongate_all (comm_state & state,
   static Timer timer ("ref_reflux_prolongate_all");
   timer.start ();
   // Require same times
-  static_assert (abs(0.1) > 0, "Function CarpetLib::abs has wrong signature");
-  assert (abs(t.get_time(ml,rl,tl) - t.get_time(ml,rl-1,tl)) <=
-          1.0e-8 * (1.0 + abs(t.get_time(ml,rl,tl))));
+  assert (fabs(t.get_time(ml,rl,tl) - t.get_time(ml,rl-1,tl)) <=
+          1.0e-8 * (1.0 + fabs(t.get_time(ml,rl,tl))));
   islab slabinfo;
   slabinfo.is_centered = 1 - ivect::dir(dir);
   transfer_from_all (state,
diff --git a/Carpet/CarpetLib/src/interpolate_3d_2tl.cc b/Carpet/CarpetLib/src/interpolate_3d_2tl.cc
index ddfe15254..33c00d8fd 100644
--- a/Carpet/CarpetLib/src/interpolate_3d_2tl.cc
+++ b/Carpet/CarpetLib/src/interpolate_3d_2tl.cc
@@ -119,7 +119,7 @@ namespace CarpetLib {
     // Linear (first order) interpolation
     
     RT const eps = 1.0e-10;
-    if (abs (t1 - t2) < eps) {
+    if (fabs (t1 - t2) < eps) {
       CCTK_WARN (0, "Internal error: arrays have same time");
     }
     if (t < min (t1, t2) - eps or t > max (t1, t2) + eps) {
diff --git a/Carpet/CarpetLib/src/interpolate_3d_3tl.cc b/Carpet/CarpetLib/src/interpolate_3d_3tl.cc
index ce8249998..f6dc7790d 100644
--- a/Carpet/CarpetLib/src/interpolate_3d_3tl.cc
+++ b/Carpet/CarpetLib/src/interpolate_3d_3tl.cc
@@ -122,7 +122,7 @@ namespace CarpetLib {
     
     RT const eps = 1.0e-10;
     
-    if (abs (t1 - t2) < eps or abs (t1 - t3) < eps or abs (t2 - t3) < eps) {
+    if (fabs (t1 - t2) < eps or fabs (t1 - t3) < eps or fabs (t2 - t3) < eps) {
       CCTK_WARN (0, "Internal error: arrays have same time");
     }
     if (t < min (min (t1, t2), t3) - eps or t > max (max (t1, t2), t3) + eps) {
diff --git a/Carpet/CarpetLib/src/interpolate_3d_4tl.cc b/Carpet/CarpetLib/src/interpolate_3d_4tl.cc
index 74bca1837..4f52f5723 100644
--- a/Carpet/CarpetLib/src/interpolate_3d_4tl.cc
+++ b/Carpet/CarpetLib/src/interpolate_3d_4tl.cc
@@ -124,8 +124,8 @@ namespace CarpetLib {
     
     RT const eps = 1.0e-10;
     
-    if (abs (t1 - t2) < eps or abs (t1 - t3) < eps or abs (t1 - t4) < eps or
-        abs (t2 - t3) < eps or abs (t2 - t4) < eps or abs (t3 - t4) < eps)
+    if (fabs (t1 - t2) < eps or fabs (t1 - t3) < eps or fabs (t1 - t4) < eps or
+        fabs (t2 - t3) < eps or fabs (t2 - t4) < eps or fabs (t3 - t4) < eps)
     {
       CCTK_WARN (0, "Internal error: arrays have same time");
     }
diff --git a/Carpet/CarpetLib/src/interpolate_3d_5tl.cc b/Carpet/CarpetLib/src/interpolate_3d_5tl.cc
index 74913a131..0cfa144ca 100644
--- a/Carpet/CarpetLib/src/interpolate_3d_5tl.cc
+++ b/Carpet/CarpetLib/src/interpolate_3d_5tl.cc
@@ -126,10 +126,10 @@ namespace CarpetLib {
     
     RT const eps = 1.0e-10;
     
-    if (abs (t1 - t2) < eps or abs (t1 - t3) < eps or abs (t1 - t4) < eps or
-        abs (t1 - t5) < eps or abs (t2 - t3) < eps or abs (t2 - t4) < eps or
-        abs (t2 - t5) < eps or abs (t3 - t4) < eps or abs (t3 - t5) < eps or
-        abs (t4 - t5) < eps)
+    if (fabs (t1 - t2) < eps or fabs (t1 - t3) < eps or fabs (t1 - t4) < eps or
+        fabs (t1 - t5) < eps or fabs (t2 - t3) < eps or fabs (t2 - t4) < eps or
+        fabs (t2 - t5) < eps or fabs (t3 - t4) < eps or fabs (t3 - t5) < eps or
+        fabs (t4 - t5) < eps)
     {
       CCTK_WARN (0, "Internal error: arrays have same time");
     }
diff --git a/Carpet/CarpetLib/src/interpolate_eno_3d_3tl.cc b/Carpet/CarpetLib/src/interpolate_eno_3d_3tl.cc
index 5d56ed3e8..ac04121c0 100644
--- a/Carpet/CarpetLib/src/interpolate_eno_3d_3tl.cc
+++ b/Carpet/CarpetLib/src/interpolate_eno_3d_3tl.cc
@@ -142,7 +142,7 @@ namespace CarpetLib {
     RT const tmax = max3 (t1, t2, t3);
     RT const eps = 1.0e-10 * (tmax - tmin);
     
-    if (abs (t1 - t2) < eps or abs (t1 - t3) < eps or abs (t2 - t3) < eps) {
+    if (fabs (t1 - t2) < eps or fabs (t1 - t3) < eps or fabs (t2 - t3) < eps) {
       CCTK_WARN (0, "Internal error: arrays have same time");
     }
     if (t < min3 (t1, t2, t3) - eps or t > max3 (t1, t2, t3) + eps) {
diff --git a/Carpet/CarpetLib/src/prolongate_3d_cc_eno_rf2.cc b/Carpet/CarpetLib/src/prolongate_3d_cc_eno_rf2.cc
index 7b5c6cee2..320635fd4 100644
--- a/Carpet/CarpetLib/src/prolongate_3d_cc_eno_rf2.cc
+++ b/Carpet/CarpetLib/src/prolongate_3d_cc_eno_rf2.cc
@@ -150,7 +150,7 @@ namespace CarpetLib {
         RT const x0 = RT(0.25) + di * RT(0.5);
         //cout << "x0=" << x0 << endl;
         RT const y0 = ipow (x0, n);
-        if (not (abs (res - y0) < 1.0e-12)) {
+        if (not (fabs (res - y0) < 1.0e-12)) {
           RT rt;
           ostringstream buf;
           buf << "Error in prolongate_3d_cc_rf2::coeffs_3d_cc_rf2\n"
@@ -260,7 +260,7 @@ namespace CarpetLib {
   {
      if (a * b < 0)
         return T(0);
-     else if (abs(a) < abs(b))
+     else if (fabs(a) < fabs(b))
         return a;
      else
         return b;
@@ -327,7 +327,7 @@ namespace CarpetLib {
 	 // check that divided differences do not change sign: if so go back to first order!
 	 if (lV*rV <= 0)
 	 // if minmod linear slope is smaller than high-order left and right undivided differences, use lowest-order TVD interpolation!
-	 //if (abs(slope) < abs(lV) || abs(slope) < abs(rV))
+	 //if (fabs(slope) < fabs(lV) || fabs(slope) < fabs(rV))
 	 {
 	    // switch back to first order TVD scheme!
 	    res = 0;
@@ -348,7 +348,7 @@ namespace CarpetLib {
 	    break;
 	 }
 	 
-	 if (abs(lV) < abs(rV)) {
+	 if (fabs(lV) < fabs(rV)) {
 	    //cout << "left ";
 	    // use left-shifted stencil since it is smoother
 	    for (ptrdiff_t i=lcoeffs::imin; i<lcoeffs::imax; ++i) {
@@ -390,7 +390,7 @@ namespace CarpetLib {
 	 // check that divided differences do not change sign: if so go back to first order!
 	 if (V[0]*V[2] <= 0 || V[0]*V[1] <= 0 || V[1]*V[2] <= 0)
 	 // if minmod linear slope is smaller than high-order left and right undivided differences, use lowest-order TVD interpolation!
-	 //if (abs(slope) < abs(V[0]) || abs(slope) < abs(V[1]) || abs(slope) < abs(V[2]))
+	 //if (fabs(slope) < fabs(V[0]) || fabs(slope) < fabs(V[1]) || fabs(slope) < fabs(V[2]))
 	 {
 	    // switch back to first order!
 	    res = 0;
@@ -406,7 +406,7 @@ namespace CarpetLib {
 	 
 	 int min = 1;  // start off with centered stencil
 	 for (int i=0; i < 3; ++i)
-	    if (abs(V[i]) < abs(V[min])) min = i;
+	    if (fabs(V[i]) < fabs(V[min])) min = i;
 	 
 	 switch (min) {
 	    case 0:
@@ -532,7 +532,7 @@ namespace CarpetLib {
 	 // check that divided differences do not change sign: if so go back to first order!
 	 if (lV*rV <= 0)
 	 // if minmod linear slope is smaller than high-order left and right undivided differences, use lowest-order TVD interpolation!
-	 //if (abs(slope) < abs(lV) || abs(slope) < abs(rV))
+	 //if (fabs(slope) < fabs(lV) || fabs(slope) < fabs(rV))
 	 {
 	    //res = 0;
             //typedef coeffs1d<RT,1,dj,0> coeffs1;
@@ -555,7 +555,7 @@ namespace CarpetLib {
 	    break;
 	 }
 	 
-	 if (abs(lV) < abs(rV)) {
+	 if (fabs(lV) < fabs(rV)) {
 	    // use left-shifted stencil since it is smoother
 	    for (ptrdiff_t i=lcoeffs::imin; i<lcoeffs::imax; ++i) {
 	       res += lcoeffs::get(i) * f[i-lcoeffs::minimin]; //interp1<T,ORDER,di> (p + i*d2, d1);
@@ -593,7 +593,7 @@ namespace CarpetLib {
 	 // check that divided differences do not change sign: if so go back to first order!
 	 if (V[0]*V[2] <= 0 || V[0]*V[1] <= 0 || V[1]*V[2] <= 0)
 	 // if minmod linear slope is smaller than high-order left and right undivided differences, use lowest-order TVD interpolation!
-	 //if (abs(slope) < abs(V[0]) || abs(slope) < abs(V[1]) || abs(slope) < abs(V[2]))
+	 //if (fabs(slope) < fabs(V[0]) || fabs(slope) < fabs(V[1]) || fabs(slope) < fabs(V[2]))
 	 {
 	    // switch back to first order!
 	    res = 0;
@@ -609,7 +609,7 @@ namespace CarpetLib {
 	 
 	 int min = 1;
 	 for (int i=0; i < 3; ++i)
-	    if (abs(V[i]) < abs(V[min])) min = i;
+	    if (fabs(V[i]) < fabs(V[min])) min = i;
 	 
 	 switch (min) {
 	    case 0:
@@ -736,7 +736,7 @@ namespace CarpetLib {
 	 // check that divided differences do not change sign: if so go back to first order!
 	 if (lV*rV <= 0)
 	 // if minmod linear slope is smaller than high-order left and right undivided differences, use lowest-order TVD interpolation!
-	 //if (abs(slope) < abs(lV) || abs(slope) < abs(rV))
+	 //if (fabs(slope) < fabs(lV) || fabs(slope) < fabs(rV))
 	 {
 	    //res = 0;
             //typedef coeffs1d<RT,1,dk,0> coeffs1;
@@ -759,7 +759,7 @@ namespace CarpetLib {
 	    break;
 	 }
 	 
-	 if (abs(lV) < abs(rV)) {
+	 if (fabs(lV) < fabs(rV)) {
 	    // use left-shifted stencil since it is smoother
 	    for (ptrdiff_t i=lcoeffs::imin; i<lcoeffs::imax; ++i) {
 	       res += lcoeffs::get(i) * f[i-lcoeffs::minimin]; //interp2<T,ORDER,di,dj> (p + i*d3, d1, d2);
@@ -797,7 +797,7 @@ namespace CarpetLib {
 	 // check that divided differences do not change sign: if so go back to first order!
 	 if (V[0]*V[2] <= 0 || V[0]*V[1] <= 0 || V[1]*V[2] <= 0)
 	 // if minmod linear slope is smaller than high-order left and right undivided differences, use lowest-order TVD interpolation!
-	 //if (abs(slope) < abs(V[0]) || abs(slope) < abs(V[1]) || abs(slope) < abs(V[2]))
+	 //if (fabs(slope) < fabs(V[0]) || fabs(slope) < fabs(V[1]) || fabs(slope) < fabs(V[2]))
 	 {
 	    // switch back to first order!
 	    res = 0;
@@ -813,7 +813,7 @@ namespace CarpetLib {
 	 
 	 int min = 1;
 	 for (int i=0; i < 3; ++i)
-	    if (abs(V[i]) < abs(V[min])) min = i;
+	    if (fabs(V[i]) < fabs(V[min])) min = i;
 	 
 	 switch (min) {
 	    case 0:
diff --git a/Carpet/CarpetLib/src/prolongate_3d_cc_enovol_rf2.cc b/Carpet/CarpetLib/src/prolongate_3d_cc_enovol_rf2.cc
index 45a6ad25b..1b230d387 100644
--- a/Carpet/CarpetLib/src/prolongate_3d_cc_enovol_rf2.cc
+++ b/Carpet/CarpetLib/src/prolongate_3d_cc_enovol_rf2.cc
@@ -153,7 +153,7 @@ namespace CarpetLib {
         RT const x1 = RT(0.5);
         //cout << "x0=" << x0 << endl;
         RT const y0 = (ipow (x1, n+1) - ipow (x0, n+1)) / ((x1-x0)*(n+1));
-        if (not (abs (res - y0) < 1.0e-12)) {
+        if (not (fabs (res - y0) < 1.0e-12)) {
           RT rt;
           ostringstream buf;
           buf << "Error in prolongate_3d_cc_rf2::coeffs_3d_cc_rf2\n"
@@ -224,7 +224,7 @@ namespace CarpetLib {
   {
      if (a * b < 0)
         return T(0);
-     else if (abs(a) < abs(b))
+     else if (fabs(a) < fabs(b))
         return a;
      else
         return b;
@@ -285,7 +285,7 @@ namespace CarpetLib {
 	 } else {
            int min = 1;
            for (int i=0; i < 3; ++i)
-              if (abs(V[i]) < abs(V[min])) min = i;
+              if (fabs(V[i]) < fabs(V[min])) min = i;
            
            switch (min) {
               case 0:
@@ -418,7 +418,7 @@ namespace CarpetLib {
     RT const x1 = di * RT(0.5) + RT(0.5);
     //cout << "x0=" << x0 << endl;
     RT const y0 = (ipow (x1, ORDER) - ipow (x0, ORDER)) / ((x1-x0)*(ORDER));
-    if (not (abs (res - y0) < 1.0e-12)) {
+    if (not (fabs (res - y0) < 1.0e-12)) {
       RT rt;
       ostringstream buf;
       buf << "Error in prolongate_3d_cc_rf2::coeffs_3d_cc_rf2\n"
diff --git a/Carpet/CarpetLib/src/prolongate_3d_cc_rf2.cc b/Carpet/CarpetLib/src/prolongate_3d_cc_rf2.cc
index 269b34947..cf8246556 100644
--- a/Carpet/CarpetLib/src/prolongate_3d_cc_rf2.cc
+++ b/Carpet/CarpetLib/src/prolongate_3d_cc_rf2.cc
@@ -121,7 +121,7 @@ namespace CarpetLib {
         RT const y0 = ipow (x0, n);
         // Allow losing 3 digits:
         CCTK_REAL const eps = RT(1.0e+3) * numeric_limits<RT>::epsilon();
-        if (not (abs (res - y0) < eps)) {
+        if (not (fabs (res - y0) < eps)) {
           RT rt;
           ostringstream buf;
           buf << "Error in prolongate_3d_cc_rf2::coeffs_3d_cc_rf2\n"
diff --git a/Carpet/CarpetLib/src/prolongate_3d_rf2.cc b/Carpet/CarpetLib/src/prolongate_3d_rf2.cc
index 93cbb52eb..ee6652aaf 100644
--- a/Carpet/CarpetLib/src/prolongate_3d_rf2.cc
+++ b/Carpet/CarpetLib/src/prolongate_3d_rf2.cc
@@ -80,7 +80,7 @@ namespace CarpetLib {
       
       // Do not test integer operators (they should be disabled
       // anyway)
-      if (abs (RT(0.5) - 0.5) > 1.0e-5) return;
+      if (fabs (RT(0.5) - 0.5) > 1.0e-5) return;
       
       // Test all orders
       bool error = false;
@@ -95,7 +95,7 @@ namespace CarpetLib {
         RT const y0 = ipow (x0, n);
         // Allow losing 3 digits:
         CCTK_REAL const eps = RT(1.0e+3) * numeric_limits<RT>::epsilon();
-        if (not (abs (res - y0) < eps)) {
+        if (not (fabs (res - y0) < eps)) {
           RT rt;
           ostringstream buf;
           buf << "Error in prolongate_3d_rf2::coeffs_3d_rf2\n"
diff --git a/Carpet/CarpetLib/src/vect.hh b/Carpet/CarpetLib/src/vect.hh
index d39327130..8a1e1b6f8 100644
--- a/Carpet/CarpetLib/src/vect.hh
+++ b/Carpet/CarpetLib/src/vect.hh
@@ -373,6 +373,7 @@ inline vect<T,D> ipow (const vect<T,D>& a, const vect<int,D>& b)
 
 DECLARE_FUNCTION_1 (abs)
 DECLARE_FUNCTION_1 (ceil)
+DECLARE_FUNCTION_1 (fabs)
 DECLARE_FUNCTION_1 (floor)
 DECLARE_FUNCTION_1 (sqrt)
 namespace std {