Simplify logic in adaptive step size control

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

1 files changed, 56 insertions, 36 deletions

diff --git a/src/StepSize.c b/src/StepSize.c
index b20c752..4ea17ff 100644
--- a/src/StepSize.c
+++ b/src/StepSize.c
@@ -10,6 +10,7 @@
 
 #include <assert.h>
 #include <math.h>
+#include <stdio.h>
 #include <stdlib.h>
 
 #include "cctk.h"
@@ -137,53 +138,53 @@ void MoL_FindAdaptiveError(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
 
-  CCTK_REAL const * restrict UpdateVar;
-  CCTK_REAL const * restrict RHSVar;
-  CCTK_REAL const * restrict ErrorVar;
-
-  int imin[3], imax[3];
-
-  int var;
-  int i, j, k;
-  int d;
-
+  /* TODO: exclude symmetry boundaries */
   assert (cctk_dim <= 3);
-  for (d = 0; d < cctk_dim; d++)
+  int imin[3], imax[3];
+  for (int d = 0; d < cctk_dim; d++)
   {
     imin[d] = cctk_bbox[2*d] ? 0 : cctk_nghostzones[d];
     imax[d] = cctk_lsh[d] - (cctk_bbox[2*d+1] ? 0 : cctk_nghostzones[d]);
   }
 
   /* Calculate absolute error */
-  for (var = 0; var < MoLNumEvolvedVariables; var++)
+  CCTK_REAL local_error = 0.0;
+  for (int var = 0; var < MoLNumEvolvedVariables; var++)
   {
 
-    UpdateVar = CCTK_VarDataPtrI(cctkGH, 0, EvolvedVariableIndex[var]);
-    RHSVar = CCTK_VarDataPtrI(cctkGH, 0, RHSVariableIndex[var]);
-    ErrorVar
+    CCTK_REAL const * restrict const
+      UpdateVar = CCTK_VarDataPtrI(cctkGH, 0, EvolvedVariableIndex[var]);
+    CCTK_REAL const * restrict const
+      RHSVar = CCTK_VarDataPtrI(cctkGH, 0, RHSVariableIndex[var]);
+    CCTK_REAL const * restrict const
+      ErrorVar
       = CCTK_VarDataPtrI(cctkGH, 0, 
                          CCTK_FirstVarIndex("MOL::ERRORESTIMATE") + var);
 
+    CCTK_REAL const rhs_relative_error
+      = maximum_relative_error * RHS_error_weight * (*Original_Delta_Time);
+
     assert (cctk_dim == 3);
-    for (k = imin[2]; k < imax[2]; k++)
+#pragma omp parallel for reduction(+: local_error)
+    for (int k = imin[2]; k < imax[2]; k++)
     {
-      for (j = imin[1]; j < imax[1]; j++)
+      for (int j = imin[1]; j < imax[1]; j++)
       {
-        for (i = imin[0]; i < imax[0]; i++)
+        for (int i = imin[0]; i < imax[0]; i++)
         {
           int const index = CCTK_GFINDEX3D(cctkGH, i, j, k);
           CCTK_REAL const scale
             = (square(maximum_absolute_error)
                + square(maximum_relative_error * UpdateVar[index])
-               + square(maximum_relative_error * RHS_error_weight
-                        * (*Original_Delta_Time) * RHSVar[index]));
-          *Error += square(ErrorVar[index]) / scale;
+               + square(rhs_relative_error * RHSVar[index]));
+          local_error += square(ErrorVar[index]) / scale;
         }
       }
     }
 
   } /* for var */
 
+  *Error += local_error;
   *Count
     += (MoLNumEvolvedVariables
         * (imax[0] - imin[0]) * (imax[1] - imin[1]) * (imax[2] - imin[2]));
@@ -234,27 +235,34 @@ void MoL_ReduceAdaptiveError(CCTK_ARGUMENTS)
   {
     CCTK_VInfo (CCTK_THORNSTRING, "Integration accuracy quotient is %g", (double)*Error);
   }
+  if (! isfinite(*Error))
+  {
+    CCTK_VWarn (CCTK_WARN_ALERT,__LINE__,__FILE__,CCTK_THORNSTRING,
+                "Integration accuracy quotient is %g, which is not a finite number -- reducing the step size", (double)*Error);
+  }
 
   if ( CCTK_EQUALS(ODE_Method,"RK45") )
   {
-     p1 = 0.2;
-     p2 = 0.25;
+    p1 = 1.0/5.0;
+    p2 = 1.0/4.0;
   }
-
-  if ( CCTK_EQUALS(ODE_Method,"RK65") )
+  else if ( CCTK_EQUALS(ODE_Method,"RK65") )
   {
-     p1 = 1.0/6.0;
-     p2 = 0.2;
+    p1 = 1.0/6.0;
+    p2 = 1.0/5.0;
   }
-
-  if ( CCTK_EQUALS(ODE_Method,"RK87") )
+  else if ( CCTK_EQUALS(ODE_Method,"RK87") )
   {
-     p1 = 0.125;
-     p2 = 1.0/7.0;
+    p1 = 1.0/8.0;
+    p2 = 1.0/7.0;
+  }
+  else
+  {
+    CCTK_WARN (CCTK_WARN_ABORT, "unsupported ODE_Method in stepsize control");
   }
 
   /* Decide whether to accept this step */
-  *MoL_Stepsize_Bad = *Error > 1;
+  *MoL_Stepsize_Bad = ! isfinite(*Error) || *Error > 1;
 
   if (*MoL_Stepsize_Bad)
   {
@@ -266,12 +274,19 @@ void MoL_ReduceAdaptiveError(CCTK_ARGUMENTS)
       CCTK_VInfo (CCTK_THORNSTRING, "*** REJECTING TIME STEP ***");
     }
 
-    cctkGH->cctk_delta_time
-      = safety_factor * (*Original_Delta_Time) / pow(*Error, p1);
-    if (! CCTK_EQUALS(verbose, "none"))
+    if (isfinite(*Error))
     {
-      CCTK_VInfo (CCTK_THORNSTRING, "Setting time step to %g", (double)cctkGH->cctk_delta_time);
+      cctkGH->cctk_delta_time
+        = safety_factor * (*Original_Delta_Time) / pow(*Error, p1);
     }
+    else
+    {
+      cctkGH->cctk_delta_time = (*Original_Delta_Time) / maximum_decrease;
+    }
+    /* if (! CCTK_EQUALS(verbose, "none")) */
+    /* { */
+    /*   CCTK_VInfo (CCTK_THORNSTRING, "Setting time step to %g", (double)cctkGH->cctk_delta_time); */
+    /* } */
 
     if (cctkGH->cctk_delta_time < (*Original_Delta_Time) / maximum_decrease)
     {
@@ -282,6 +297,11 @@ void MoL_ReduceAdaptiveError(CCTK_ARGUMENTS)
         CCTK_VInfo (CCTK_THORNSTRING, "   Time step reduction too large; clamping time step to %g", (double)cctkGH->cctk_delta_time);
       }
     }
+    if (cctkGH->cctk_delta_time == (CCTK_REAL)0.0)
+      /* yes, we want to compare to zero exactly, to catch underflows */
+      {
+        CCTK_WARN (CCTK_WARN_ABORT, "New step size would be zero -- aborting");
+      }
 
     cctkGH->cctk_time += cctkGH->cctk_delta_time;
   }