path: root/src/NaNCheck.c

/*@@
  @file       NaNCheck.c
   @date      Sat 21 Apr 2001
   @author    Thomas Radke
   @desc
              Routines to check CCTK real and complex variables
              against Not-a-Number values
   @enddesc
   @version   $Id$
 @@*/

#include <math.h>
#include <stdlib.h>
#include <string.h>

#include "cctk.h"
#include "cctk_WarnLevel.h"
#include "cctk_Parameters.h"
#include "cctk_Termination.h"

/* the rcsid and the macro to use it */
static char *rcsid = "$Header$";
CCTK_FILEVERSION(BetaThorns_NaNChecker_NaNCheck_c)


/********************************************************************
 ********************    External Routines   ************************
 ********************************************************************/
int NaNChecker_NaNCheck (cGH *GH);


/********************************************************************
 ********************    Internal Routines   ************************
 ********************************************************************/
static void NaNCheck (int vindex, const char *optstring, void *arg);
static void PrintWarning (const char *error_type,
                          int linear_index,
                          int fp_type,
                          const char *fullname,
                          const cGroupDynamicData *gdata);


 /*@@
   @routine    NaNChecker_NaNCheck
   @author     Thomas Radke
   @date       Sat 21 Apr 2001
   @desc
               If it is time to check, loop through all variables to check
               for NaN values
   @enddesc
   @calls      CCTK_TraverseString
               NaNCheck

   @var        GH
   @vdesc      Pointer to CCTK GH
   @vtype      cGH *
   @vio        in
   @endvar

   @returntype int
   @returndesc
               0 if checking was not done during current iteration, or
               return code of @seeroutine CCTK_TraverseString otherwise
   @endreturndesc
@@*/
int NaNChecker_NaNCheck (cGH *GH)
{
  DECLARE_CCTK_PARAMETERS
  int retval;


  if (GH->cctk_iteration % check_every == 0)
  {
    retval = CCTK_TraverseString (check_vars, NaNCheck, GH, CCTK_GROUP_OR_VAR);
  }
  else
  {
    retval = 0;
  }

  return (retval);
}


/********************************************************************
 ********************    Internal Routines   ************************
 ********************************************************************/

 /*@@
   @routine    PrintWarning
   @author     Thomas Radke
   @date       Sat 21 Apr 2001
   @desc
               Prints a level 1 warning for a Inf/NaN found in a variable
               at the given processor-local linear index and terminates
               Cactus if requested.
               The warning includes the variable's fullname along with
               the global index of the NaN element in fortran order.
   @enddesc
   @calls      CCTK_VWarn

   @var        error_type
   @vdesc      string containing the error value found (NaN or Inf)
   @vtype      const char *
   @vio        in
   @endvar
   @var        fp_type
   @vdesc      indicates if variable of of real or complex type
   @vtype      int
   @vio        in
   @endvar
   @var        linear_index
   @vdesc      processor-local linear index of the NaN/Inf
   @vtype      int
   @vio        in
   @endvar
   @var        fullname
   @vdesc      full name of the variable
   @vtype      const char *
   @vio        in
   @endvar
   @var        gdata
   @vdesc      size information on the variable to compute the global index
   @vtype      const cGroupDynamicData *
   @vio        in
   @endvar
@@*/
static void PrintWarning (const char *error_type,
                          int linear_index,
                          int fp_type,
                          const char *fullname,
                          const cGroupDynamicData *gdata)
{
  DECLARE_CCTK_PARAMETERS
  int i;
  char *index_string;
  const char *complex_part;


  if (fp_type == 2)
  {
    complex_part = linear_index & 1 ? "imag part of " : "real part of ";
    linear_index /= 2;
  }
  else
  {
    complex_part = "";
  }

  if (gdata->dim == 0)
  {
    CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
                "%s caught in %svariable '%s'",
                error_type, complex_part, fullname);
  }
  else
  {
    /* assume max. 10 characters per index number (including separators) */
    index_string = (char *) malloc (10 * gdata->dim);

    sprintf (index_string, "%d",
             (linear_index % gdata->lsh[0]) + gdata->lbnd[0] + 1);
    for (i = 1; i < gdata->dim; i++)
    {
      linear_index /= gdata->lsh[i - 1];
      sprintf (index_string, "%s, %d", index_string,
               (linear_index % gdata->lsh[i]) + gdata->lbnd[i] + 1);
    }

    CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
                "%s caught in %svariable '%s' at (%s)",
                error_type, complex_part, fullname, index_string);

    free (index_string);
  }

  /* Do more than just print a warning ? */
  if (CCTK_Equals (action_if_found, "terminate"))
  {
    CCTK_TerminateNext (NULL);
  }
  else if (CCTK_Equals (action_if_found, "abort"))
  {
    CCTK_Abort (NULL, 0);
  }
}


 /*@@
   @routine    CHECK_DATA
   @date       Sat 21 Apr 2001
   @author     Thomas Radke
   @desc
               Macro to check a given typed array against NaN's.
               If finite(3) is available on the system it will also check
               for Inf values.
   @enddesc
   @calls      PrintWarning

   @var        cctk_type
   @vdesc      CCTK variable type of variable to check
   @vtype      <CCTK vartype>
   @vio        in
   @endvar
@@*/
#ifdef HAVE_FINITE

#define CHECK_DATA(cctk_type)                                                 \
{                                                                             \
  int _i;                                                                     \
  const cctk_type *_typed_data = (const cctk_type *) data;                    \
                                                                              \
                                                                              \
  /* now loop over all elements and check against NaN's */                    \
  for (_i = 0; _i < nelems; _i++)                                             \
  {                                                                           \
    if (! finite ((double) _typed_data[_i]))                                  \
    {                                                                         \
      PrintWarning (isnan ((double) _typed_data[_i]) ? "NaN" : "Inf",         \
                    _i, fp_type, fullname, &gdata);                           \
    }                                                                         \
  }                                                                           \
}

#else

#define CHECK_DATA(cctk_type)                                                 \
{                                                                             \
  int _i;                                                                     \
  const cctk_type *_typed_data = (const cctk_type *) data;                    \
                                                                              \
                                                                              \
  /* now loop over all elements and check against NaN's */                    \
  for (_i = 0; _i < nelems; _i++)                                             \
  {                                                                           \
    if (isnan ((double) _typed_data[_i]))                                     \
    {                                                                         \
      PrintWarning ("NaN", _i, fp_type, fullname, &gdata);                    \
    }                                                                         \
  }                                                                           \
}

#endif /* HAVE_FINITE */


 /*@@
   @routine    NaNCheck
   @date       Sat 21 Apr 2001
   @author     Thomas Radke
   @desc
               Checks a CCTK variable given by its index against NaN's.
               Only floating point typed variables are checked.
               <BR>
               This routine is called as a callback via CCTK_TraverseString().
   @enddesc
   @calls      CHECK_DATA

   @var        vindex
   @vdesc      index of variable to check
   @vtype      int
   @vio        in
   @endvar
   @var        optstring
   @vdesc      optional string appended to the group/variable name
   @vtype      unused
   @vio        in
   @endvar
   @var        _GH
   @vdesc      Pointer to CCTK GH
   @vtype      void *
   @vio        in
   @endvar
@@*/
static void NaNCheck (int vindex, const char *optstring, void *_GH)
{
  cGH *GH;
  int i, fp_type;
  int vtype, gindex, nelems;
  char *fullname;
  const char *vtypename;
  cGroupDynamicData gdata;
  const void *data;


  /* avoid compiler warning about unused parameters */
  optstring = optstring;

  GH = (cGH *) _GH;

  vtype = CCTK_VarTypeI (vindex);
  fullname = CCTK_FullName (vindex);

  /* check if the variable type is some floating point */
  vtypename = CCTK_VarTypeName (vtype);
  if (strncmp (vtypename, "CCTK_VARIABLE_REAL", 18) == 0)
  {
    fp_type = 1;
  }
  else if (strncmp (vtypename, "CCTK_VARIABLE_COMPLEX", 22) == 0)
  {
    fp_type = 2;
  }
  else
  {
    fp_type = 0;
  }

  if (fp_type)
  {
    gindex = CCTK_GroupIndexFromVarI (vindex);

    /* check if variable has storage assigned */
    if (CCTK_QueryGroupStorageI (GH, gindex))
    {
      /* get the number of elements to check for this variable */
      nelems = 1;
      gdata.dim = 0;
      if (CCTK_GroupTypeI (vindex) != CCTK_SCALAR)
      {
        CCTK_GroupDynamicData (GH, gindex, &gdata);
        for (i = 0; i < gdata.dim; i++)
        {
          nelems *= gdata.lsh[i];
        }
      }

      /* simply double the number of elements for a CCTK_COMPLEX variable */
      if (fp_type == 2)
      {
        nelems *= 2;
      }

      /* get the pointer to the data (current time level) */
      data = CCTK_VarDataPtrI (GH, CCTK_NumTimeLevelsFromVarI(vindex)-1,vindex);

      /* do the checking according to the variable's type */
      if (vtype == CCTK_VARIABLE_REAL || CCTK_VARIABLE_COMPLEX)
      {
        CHECK_DATA (CCTK_REAL);
      }
#ifdef CCTK_REAL4
      else if (vtype == CCTK_VARIABLE_REAL4 || CCTK_VARIABLE_COMPLEX8)
      {
        CHECK_DATA (CCTK_REAL4);
      }
#endif
#ifdef CCTK_REAL8
      else if (vtype == CCTK_VARIABLE_REAL8 || CCTK_VARIABLE_COMPLEX16)
      {
        CHECK_DATA (CCTK_REAL8);
      }
#endif
#ifdef CCTK_REAL16
      else if (vtype == CCTK_VARIABLE_REAL16 || CCTK_VARIABLE_COMPLEX32)
      {
        CHECK_DATA (CCTK_REAL16);
      }
#endif
      else
      {
        CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
                    "NanCheck: Unknown variable type '%s' for variable '%s'",
                    vtypename, fullname);
      }
    }
    else
    {
      CCTK_VWarn (3, __LINE__, __FILE__, CCTK_THORNSTRING,
                  "NaNCheck: Ignoring variable '%s' (no storage)", fullname);
    }
  }
  else
  {
    CCTK_VWarn (9, __LINE__, __FILE__, CCTK_THORNSTRING,
                "NaNCheck: Ignoring variable '%s' "
                "(not a floating-point variable)", fullname);
  }

  free (fullname);
}