Improve roundoff accuracy of efficient RK4

I just encountered a weird case where evolving Minkowski with the
efficient RK4 implementation was not static. It turns out that the
reason was floating-point round-off error in RK4. During each of the
substeps, RK4 adds terms with factors of 1/3 and 2/3, and for the
final step, adds a term with a factor of 4/3. Due to round-off, these
terms do not cancel exactly, so that the time-evolved lapse is
slightly different from 1 (although all the lapse RHS terms are 0).
The same goes for the diagonal terms of the metric.

The attached patch keeps track of the round-off error when adding
these terms, leading to a more accurate integration result.

The problem described above is not always visible, but also depends on
optimisation settings.


git-svn-id: http://svn.cactuscode.org/arrangements/CactusNumerical/MoL/trunk@145 578cdeb0-5ea1-4b81-8215-5a3b8777ee0b

1 files changed, 15 insertions, 4 deletions

diff --git a/src/RK4.c b/src/RK4.c
index 250c847..55a4a4e 100644
--- a/src/RK4.c
+++ b/src/RK4.c
@@ -86,8 +86,8 @@ CCTK_WARN(0, "not implemented");
 #endif
 
   static CCTK_INT scratchspace_firstindex = -99;
-  CCTK_INT index, var, scratchstep, alphaindex, scratchindex;
-  CCTK_INT totalsize, singlearraysize;
+  CCTK_INT index, var;
+  CCTK_INT totalsize;
   CCTK_REAL alpha, beta;
   CCTK_REAL * restrict UpdateVar;
   CCTK_REAL * restrict OldVar;
@@ -96,6 +96,11 @@ CCTK_WARN(0, "not implemented");
 
   CCTK_INT arrayscratchlocation;
 
+  /* Keep a running total of alpha as we perform the substeps, so that
+     we know the "real" alpha (including round-off errors) when we
+     calculate the final result. */
+  static CCTK_REAL sum_alpha;
+
   totalsize = 1;
   for (arraydim = 0; arraydim < cctk_dim; arraydim++)
   {
@@ -126,6 +131,12 @@ CCTK_WARN(0, "not implemented");
       beta  = 1.0 / 6.0;
   }
 
+  if (MoL_Intermediate_Steps == (*MoL_Intermediate_Step))
+  {
+    sum_alpha = 0.0;
+  }
+  sum_alpha += alpha;
+
   /* FIXME */
 
 #ifdef MOLDOESCOMPLEX
@@ -197,7 +208,7 @@ CCTK_WARN(0, "not implemented");
 #pragma omp parallel for
       for (index = 0; index < totalsize; index++)
       {
-        UpdateVar[index] += ScratchVar[index] - 4.0 / 3.0 * OldVar[index];
+        UpdateVar[index] += ScratchVar[index] - (sum_alpha - 1.0) * OldVar[index];
       }
     }
 
@@ -276,7 +287,7 @@ CCTK_WARN(0, "not done");
         UpdateVar[index] += ScratchVar[index] - 4.0 / 3.0 * OldVar[index];
       }
     }
-    arrayscratchlocation +=  arraytotalsize;
+    arrayscratchlocation += arraytotalsize;
   }
 
   return;