Added an integer grid function NaNChecker::NaNmask which is used as a bitmask

to mark locations of NaNs found in grid functions. The mask is output as an
HDF5 file and can then be visualized in OpenDX.
Setting the mask and outputting it is activated by default. It can be switched
off with NaNChecker::out_NaNmask = "no".


git-svn-id: http://svn.cactuscode.org/arrangements/CactusUtils/NaNChecker/trunk@39 ff385933-4943-42dc-877b-ffc776028de6

4 files changed, 168 insertions, 62 deletions

diff --git a/interface.ccl b/interface.ccl
index 495a949..b65bb8e 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -4,3 +4,10 @@
 implements: NaNChecker
 
 INCLUDES HEADER: NaNCheck.h in NaNChecker.h
+
+private:
+
+INT NaNmask_group TYPE = GF
+{
+  NaNmask
+} "Grid function mask for NaN locations"
diff --git a/param.ccl b/param.ccl
index ae06c46..56e1892 100644
--- a/param.ccl
+++ b/param.ccl
@@ -30,6 +30,10 @@ KEYWORD check_for "Check for NaN's and/or infinite numbers (only evaluated if fi
   "both" :: "Check for both NaN's and infinite numbers"
 } "both"
 
+BOOLEAN out_NaNmask "Dump the NaN grid function mask into an HDF5 file" STEERABLE = ALWAYS
+{
+} "yes"
+
 KEYWORD action_if_found "What to do if a NaN was found" STEERABLE = ALWAYS
 {
   "just warn" :: "Just print a level 1 warning"
@@ -41,4 +45,4 @@ KEYWORD verbose "How much information to give" STEERABLE = ALWAYS
 {
   "all"      :: "All information"
   "standard" :: "Standard information"
-} "standard"
\ No newline at end of file
+} "standard"
diff --git a/schedule.ccl b/schedule.ccl
index 9fdff83..862ca49 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -6,5 +6,6 @@ if (*check_vars && check_every > 0)
   schedule NaNChecker_NaNCheck at POSTSTEP
   {
     LANG:C
+    STORAGE: NaNmask_group
   } "Check for NaN's"
 }
diff --git a/src/NaNCheck.c b/src/NaNCheck.c
index c1176a0..4368143 100644
--- a/src/NaNCheck.c
+++ b/src/NaNCheck.c
@@ -53,7 +53,7 @@ void CCTK_FCALL CCTK_FNAME (NaNChecker_SetVarsToNaN)
 static void CheckForNaN (int vindex, const char *optstring, void *arg);
 static void SetToNaN (int vindex, const char *optstring, void *arg);
 static void PrintWarning (const char *error_type,
-			  int verbose,
+                          int verbose,
                           int linear_index,
                           int fp_type,
                           const CCTK_REAL *const coords[],
@@ -70,12 +70,20 @@ typedef struct
   int report_max;
   const char *action_if_found;
   int count;
+  CCTK_INT *NaNmask;
+  int bitmask;
 #ifdef HAVE_FINITE
   enum {CHECK_FOR_NAN = 0, CHECK_FOR_INF = 1, CHECK_FOR_BOTH = 2} check_for;
 #endif
 } t_nanchecker_info;
 
 
+/********************************************************************
+ ********************    Static Variables    ************************
+ ********************************************************************/
+static int last_iteration_output = -1;
+
+
  /*@@
    @routine    NaNChecker_NaNCheck
    @author     Thomas Radke
@@ -101,44 +109,87 @@ typedef struct
 @@*/
 int NaNChecker_NaNCheck (const cGH *GH)
 {
-  int retval;
+  int i, nelems, retval;
+  const char *method;
   t_nanchecker_info info;
   DECLARE_CCTK_PARAMETERS
 
 
-  if (GH->cctk_iteration % check_every == 0)
+  if (GH->cctk_iteration % check_every)
   {
-    info.GH = GH;
-    info.report_max = report_max;
-    info.action_if_found = action_if_found;
-    if (CCTK_EQUALS(verbose,"all"))
-    {
-      info.verbose = 1;
-    }
-    else
+    return (0);
+  }
+  if (GH->cctk_iteration == last_iteration_output)
+  {
+    CCTK_WARN (2, "Already called NaNChecker I/O method at this iteration");
+    return (0);
+  }
+
+  info.GH = GH;
+  info.count = 0;
+  info.bitmask = 0;
+  info.report_max = report_max;
+  info.action_if_found = action_if_found;
+  info.verbose = CCTK_Equals (verbose, "all");
+
+  info.NaNmask = out_NaNmask ?
+                 CCTK_VarDataPtr (GH, 0, "NaNChecker::NaNmask") : NULL;
+  if (info.NaNmask)
+  {
+    /* zero out the NaN mask */
+    for (i = 0, nelems = 1; i < GH->cctk_dim; i++)
     {
-      info.verbose = 0;
+      nelems *= GH->cctk_lsh[i];
     }
+    memset (info.NaNmask, 0, nelems * sizeof (CCTK_INT));
+  }
+
 #ifdef HAVE_FINITE
-    if (CCTK_Equals (check_for, "NaN"))
+  if (CCTK_Equals (check_for, "NaN"))
+  {
+    info.check_for = CHECK_FOR_NAN;
+  }
+  else if (CCTK_Equals (check_for, "Inf"))
+  {
+    info.check_for = CHECK_FOR_INF;
+  }
+  else
+  {
+    info.check_for = CHECK_FOR_BOTH;
+  }
+#endif
+
+  retval = CCTK_TraverseString (check_vars, CheckForNaN, &info,
+                                CCTK_GROUP_OR_VAR);
+
+  /* save the iteration of the last call */
+  last_iteration_output = GH->cctk_iteration;
+
+  if (info.count > 0)
+  {
+    /* if NaNs were found then output NaN mask with the 'IOHDF5' I/O method */
+    if (info.NaNmask && info.bitmask)
     {
-      info.check_for = CHECK_FOR_NAN;
+      if (info.verbose)
+      {
+        CCTK_INFO ("Write out NaN mask using the 'IOHDF5' I/O method");
+      }
+      CCTK_OutputVarAsByMethod (GH, "NaNChecker::NaNmask[downsample={1 1 1}]",
+                                "IOHDF5", "NaNmask");
     }
-    else if (CCTK_Equals (check_for, "Inf"))
+
+    if (CCTK_Equals (info.action_if_found, "terminate"))
     {
-      info.check_for = CHECK_FOR_INF;
+      CCTK_WARN (1, "'action_if_found' parameter is set to 'terminate' - "
+                    "scheduling graceful termination of Cactus");
+      CCTK_TerminateNext (NULL);
     }
-    else
+    else if (CCTK_Equals (info.action_if_found, "abort"))
     {
-      info.check_for = CHECK_FOR_BOTH;
+      CCTK_WARN (1, "'action_if_found' parameter is set to 'abort' - "
+                    "aborting Cactus now");
+      CCTK_Abort (NULL, 0);
     }
-#endif
-    retval = CCTK_TraverseString (check_vars, CheckForNaN, &info,
-                                  CCTK_GROUP_OR_VAR);
-  }
-  else
-  {
-    retval = 0;
   }
 
   return (retval);
@@ -191,8 +242,9 @@ int NaNChecker_NaNCheck (const cGH *GH)
    @returndesc
                the total number of NaN values found, or<BR>
                -1 if a NULL pointer was passed for 'GH' and/or 'vars',<BR>
-               -2 if an unknow keyword was passed in 'action_if_found',<BR>
-               -3 if the 'vars' string couldn't be parsed
+               -2 if NaNChecker was already called at this iteration,<BR>
+               -3 if an unknow keyword was passed in 'action_if_found',<BR>
+               -4 if the 'vars' string couldn't be parsed
    @endreturndesc
 @@*/
 int NaNChecker_CheckVarsForNaN (const cGH *GH,
@@ -208,6 +260,11 @@ int NaNChecker_CheckVarsForNaN (const cGH *GH,
   {
     CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
                 "NULL pointer passed for 'GH' and/or 'vars' argument");
+    return (-2);
+  }
+  if (GH->cctk_iteration == last_iteration_output)
+  {
+    CCTK_WARN (2, "Already called NaNChecker I/O method at this iteration");
     return (-1);
   }
 
@@ -218,13 +275,16 @@ int NaNChecker_CheckVarsForNaN (const cGH *GH,
     CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
                 "Unknown keyword '%s' for 'action_if_found' argument",
                 action_if_found);
-    return (-2);
+    return (-3);
   }
 
   info.GH = GH;
+  info.count = 0;
+  info.bitmask = 0;
+  info.verbose = 0;
   info.report_max = report_max;
   info.action_if_found = action_if_found;
-  info.count = 0;
+
 #ifdef HAVE_FINITE
   if (CCTK_Equals (check_for, "NaN"))
   {
@@ -245,11 +305,42 @@ int NaNChecker_CheckVarsForNaN (const cGH *GH,
     info.check_for = CHECK_FOR_BOTH;
   }
 #endif
+
   if (CCTK_TraverseString (vars, CheckForNaN, &info, CCTK_GROUP_OR_VAR) < 0)
   {
     CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
                 "Couldn't traverse 'vars' string '%s'", vars);
-    return (-3);
+    return (-4);
+  }
+
+  /* save the iteration of the last call */
+  last_iteration_output = GH->cctk_iteration;
+
+  if (info.count > 0)
+  {
+    /* if NaNs were found then output NaN mask with the 'IOHDF5' I/O method */
+    if (info.NaNmask && info.bitmask)
+    {
+      if (info.verbose)
+      {
+        CCTK_INFO ("Write out NaN mask using the 'IOHDF5' I/O method");
+      }
+      CCTK_OutputVarAsByMethod (GH, "NaNChecker::NaNmask[downsample={1 1 1}]",
+                                "IOHDF5", "NaNmask");
+    }
+
+    if (CCTK_Equals (info.action_if_found, "terminate"))
+    {
+      CCTK_WARN (1, "'action_if_found' parameter is set to 'terminate' - "
+                    "scheduling graceful termination of Cactus");
+      CCTK_TerminateNext (NULL);
+    }
+    else if (CCTK_Equals (info.action_if_found, "abort"))
+    {
+      CCTK_WARN (1, "'action_if_found' parameter is set to 'abort' - "
+                    "aborting Cactus now");
+      CCTK_Abort (NULL, 0);
+    }
   }
 
   return (info.count);
@@ -389,7 +480,7 @@ void CCTK_FCALL CCTK_FNAME (NaNChecker_SetVarsToNaN)
    @endvar
 @@*/
 static void PrintWarning (const char *error_type,
-			  int verbose,
+                          int verbose,
                           int linear_index,
                           int fp_type,
                           const CCTK_REAL *const coords[],
@@ -422,7 +513,7 @@ static void PrintWarning (const char *error_type,
   {
     /* assume max. 10 characters per index number and 40 characters per
        coordinate value (including separators) */
-    index_string = (char *) malloc (5 * 10 * gdata->dim);
+    index_string = malloc (5 * 10 * gdata->dim);
     coord_string = index_string + 10 * gdata->dim;
 
     amended_index = linear_index;
@@ -440,7 +531,7 @@ static void PrintWarning (const char *error_type,
                (amended_index % gdata->lsh[i]) + gdata->lbnd[i] + 1);
       if (coords)
       {
-	sprintf (coord_string, "%s, %5.3e", coord_string,
+        sprintf (coord_string, "%s, %5.3e", coord_string,
                  (double) coords[i][linear_index]);
       }
     }
@@ -449,16 +540,16 @@ static void PrintWarning (const char *error_type,
     {
       if (coords)
       {
-	CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING,
-		    "%s caught in %svariable '%s' at index (%s) with coordinates "
-		    "(%s)", error_type, complex_part, fullname, index_string,
-		    coord_string);
+        CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING,
+                    "%s caught in %svariable '%s' at index (%s) with coordinates "
+                    "(%s)", error_type, complex_part, fullname, index_string,
+                    coord_string);
       }
       else
       {
-	CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING,
-		    "%s caught in %svariable '%s' at (%s)",
-		    error_type, complex_part, fullname, index_string);
+        CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING,
+                    "%s caught in %svariable '%s' at (%s)",
+                    error_type, complex_part, fullname, index_string);
       }
 
       free (index_string);
@@ -503,7 +594,11 @@ static void PrintWarning (const char *error_type,
           (info->report_max < 0 || nans_found <= info->report_max))           \
       {                                                                       \
         PrintWarning (isnan ((double) _data[_i]) ? "NaN" : "Inf",             \
-		      info->verbose, _i, fp_type, coords, fullname, &gdata);  \
+                      info->verbose, _i, fp_type, coords, fullname, &gdata);  \
+      }                                                                       \
+      if (info->NaNmask && gtype == CCTK_GF)                                  \
+      {                                                                       \
+        info->NaNmask[_i] |= 1 << info->bitmask;                              \
       }                                                                       \
     }                                                                         \
   }                                                                           \
@@ -529,6 +624,10 @@ static void PrintWarning (const char *error_type,
         PrintWarning ("NaN", info->verbose, _i, fp_type, coords, fullname,    \
                        &gdata);                                               \
       }                                                                       \
+      if (info->NaNmask)                                                      \
+      {                                                                       \
+        info->NaNmask[_i] |= 1 << info->bitmask;                              \
+      }                                                                       \
     }                                                                         \
   }                                                                           \
 }
@@ -580,7 +679,7 @@ static void CheckForNaN (int vindex, const char *optstring, void *_info)
   const void *data;
 
 
-  info = (t_nanchecker_info *) _info;
+  info = _info;
 
   vtype = CCTK_VarTypeI (vindex);
   fullname = CCTK_FullName (vindex);
@@ -657,8 +756,7 @@ static void CheckForNaN (int vindex, const char *optstring, void *_info)
       sprintf (coord_system_name, "cart%dd", gdata.dim);
       if (CCTK_CoordSystemHandle (coord_system_name) >= 0)
       {
-        coords = (const CCTK_REAL **)
-                 malloc (gdata.dim * sizeof (const CCTK_REAL *));
+        coords = malloc (gdata.dim * sizeof (const CCTK_REAL *));
       }
     }
     for (i = 0; i < gdata.dim; i++)
@@ -668,7 +766,7 @@ static void CheckForNaN (int vindex, const char *optstring, void *_info)
       {
         coords[i] = (const CCTK_REAL *) CCTK_VarDataPtrI (info->GH, 0,
                     CCTK_CoordIndex (i + 1, NULL, coord_system_name));
-        if (! coords[i]) 
+        if (! coords[i])
         {
           free (coords);
           coords = NULL;
@@ -712,31 +810,27 @@ static void CheckForNaN (int vindex, const char *optstring, void *_info)
   }
 
   /* Do more than just print a warning ? */
-  if (nans_found > 0 && info->action_if_found)
+  if (nans_found > 0 && info->action_if_found && gdata.dim > 0)
   {
-    if (gdata.dim > 0)
-    {
-      CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
-                  "There were %d NaN/Inf value(s) found in variable '%s'",
-                  nans_found, fullname);
-    }
-
-    if (CCTK_Equals (info->action_if_found, "terminate"))
+    if (info->NaNmask && gtype == CCTK_GF)
     {
       CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
-                  "\"action_if_found\" parameter is set to \"terminate\" - "
-                  "scheduling graceful termination of Cactus");
-      CCTK_TerminateNext (NULL);
+                  "There were %d NaN/Inf value(s) found in variable '%s' "
+                  "(NaNmask bitfield %d)",
+                  nans_found, fullname, info->bitmask);
     }
-    else if (CCTK_Equals (info->action_if_found, "abort"))
+    else
     {
       CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
-                  "\"action_if_found\" parameter is set to \"abort\" - "
-                  "aborting Cactus now");
-      CCTK_Abort (NULL, 0);
+                  "There were %d NaN/Inf value(s) found in variable '%s'",
+                  nans_found, fullname);
     }
   }
   info->count += nans_found;
+  if (gtype == CCTK_GF)
+  {
+    info->bitmask++;
+  }
 
   if (coords)
   {