path: root/src/Utils.c

/*@@
   @file      Utils.c
   @date      Tue 4th July 2000
   @author    Gabrielle Allen
   @desc
              Utility routines which may be called by other IO thorns
   @enddesc 
   @version   $Id$
 @@*/


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

#include "cctk.h"
#include "cctk_Parameters.h"
#include "ioutil_Utils.h"
#include "ioutil_CheckpointRecovery.h"

/* the rcs ID and its dummy function to use it */
static char *rcsid = "$Header$";
CCTK_FILEVERSION(CactusBase_IOUtil_Utils_c)


/* uncomment this if you want some debugging output */
/* #define DEBUG_IOUTIL 1 */


/*@@
   @routine    IOUtil_1DLines
   @date       July 4 2000
   @author     Gabrielle Allen, Gerd Lanfermann, Thomas Radke
   @desc
               Fills out an array determining where to position 1D lines
               for output on a multidimensional regular Cartesian unigrid.
               The first slot of the array specifies the 1D line direction,
               the second slot fixes the indices of the starting point
               of that line on the grid.
   @enddesc

   @calls      CCTK_CoordRange

   @var        GH
   @vdesc      Pointer to CCTK GH
   @vtype      cGH *
   @vio        in
   @endvar
   @var        num_dims
   @vdesc      number of dimensions of underlying grid
   @vtype      int
   @vio        in
   @endvar
   @var        origin_index
   @vdesc      origin for 1D lines given by grid point indices
   @vtype      const int [num_dims][num_dims]
   @vio        in
   @endvar
   @var        origin_phys
   @vdesc      origin for 1D lines given by physical coordinates
   @vtype      const CCTK_REAL [num_dims][num_dims]
   @vio        in
   @endvar
   @var        slice_center
   @vdesc      resulting 1D line slice center
   @vtype      int [num_dims][num_dims]
   @vio        out
   @endvar

   @returntype int
   @returndesc
               0  for success
               -1 no coordinate system of given dimensions found
   @endreturndesc
 @@*/
int IOUtil_1DLines (cGH *GH,
                    int num_dims,
                    int *const origin_index[],
                    CCTK_REAL *const origin_phys[],
                    int *slice_center[])
{
  DECLARE_CCTK_PARAMETERS
  int dim, dir;
  char system_name[20];
  CCTK_REAL *lower_range, *upper_range;


  /* get the appropriate coordinate system name */
  sprintf (system_name, "cart%dd", num_dims);

  if (CCTK_CoordSystemHandle (system_name) < 0)
  {
    CCTK_VWarn (4, __LINE__, __FILE__, CCTK_THORNSTRING,
                "IOUtil_1DLines: Cartesian coordinate system '%s' not found",
                system_name);
    memset (slice_center, 0, num_dims * num_dims * sizeof (int));
    return (-1);
  }

  /* allocate arrays for ranges */
  lower_range = (CCTK_REAL *) calloc (2 * num_dims, sizeof (CCTK_REAL));
  upper_range = lower_range + num_dims;

  /* get the ranges in every direction */
  for (dir = 0; dir < num_dims; dir++)
  {
    if (CCTK_CoordRange (GH, &lower_range[dir], &upper_range[dir],
                         dir + 1, NULL, system_name) < 0)
    {
      CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
                  "IOUtil_1DLines: Could not get ranges for %c-direction "
                  "of coordinate system '%s'",
                  'x' + dir, system_name);
    }
  }

  /* now set the slice center for each line
     according to origin_index[] or origin_phys[] */
  for (dir = 0; dir < num_dims; dir++)
  {
    for (dim = 0; dim < num_dims; dim++)
    {
      if (dim == dir)
      {
        slice_center[dir][dim] = 0;
      }
      else if (origin_index && origin_index[dir][dim] >= 0)
      {
        slice_center[dir][dim] = origin_index[dir][dim];
      }
      else if (lower_range[dim] > origin_phys[dir][dim] || 
               upper_range[dim] < origin_phys[dir][dim])
      {
        CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING,
                    "IOUtil_1DLines: Coordinate in %c-direction (%f, %f) "
                    "doesn't contain %f",
                    'x' + dim, (double) lower_range[dim],
                    (double) upper_range[dim], (double) origin_phys[dir][dim]);
        slice_center[dir][dim] = 0;
      }
      else 
      {
        /* Find index for first point above the chosen coordinate */
        slice_center[dir][dim] =
          ceil ((origin_phys[dir][dim] - lower_range[dim]) /
                GH->cctk_delta_space[dim] - 1e-6);

#ifdef DEBUG_IOUTIL
        printf("spxyz for %d-coord of lines in %d-direction is %d\n",
               dim,dir,slice_center[dir][dim]);
#endif
      }
    }
  }

  /* free allocated resources */
  free (lower_range);

  return (0);
}


/*@@
   @routine    IOUtil_2DPlanes
   @date       July 4 2000
   @author     Gabrielle Allen, Gerd Lanfermann, Thomas Radke
   @desc
               Fills out an array determining where to position 2D planes
               for output on a multidimensional regular Cartesian unigrid.
   @enddesc

   @calls      CCTK_CoordRange

   @var        GH
   @vdesc      Pointer to CCTK GH
   @vtype      cGH *
   @vio        in
   @endvar
   @var        num_dims
   @vdesc      number of dimensions of underlying grid
   @vtype      int
   @vio        in
   @endvar
   @var        origin_index
   @vdesc      origin for 2D planes given by grid point indices
   @vtype      const int [num_dims]
   @vio        in
   @endvar
   @var        origin_phys
   @vdesc      origin for 2D planes given by physical coordinates
   @vtype      const CCTK_REAL [num_dims]
   @vio        in
   @endvar
   @var        slice_center
   @vdesc      resulting 2D plane slice center
   @vtype      int [num_dims]
   @vio        out
   @endvar

   @returntype int
   @returndesc
               0  for success
               -1 no coordinate system of given dimensions found
   @endreturndesc
 @@*/
int IOUtil_2DPlanes (cGH *GH,
                     int num_dims,
                     const int origin_index[],
                     const CCTK_REAL origin_phys[],
                     int slice_center[])
{
  DECLARE_CCTK_PARAMETERS
  int dir;
  char system_name[20];
  CCTK_REAL *lower_range, *upper_range;


  /* get the appropriate coordinate system name */
  sprintf (system_name, "cart%dd", num_dims);

  if (CCTK_CoordSystemHandle (system_name) < 0)
  {
    CCTK_VWarn (4, __LINE__, __FILE__, CCTK_THORNSTRING,
                "IOUtil_1DLines: Cartesian coordinate system '%s' not found",
                system_name);
    memset (slice_center, 0, num_dims * sizeof (int));
    return (-1);
  }

  /* allocate arrays for ranges */
  lower_range = (CCTK_REAL *) calloc (2 * num_dims, sizeof (CCTK_REAL));
  upper_range = lower_range + num_dims;

  /* get the ranges in every direction */
  for (dir = 0; dir < num_dims; dir++)
  {
    if (CCTK_CoordRange (GH, &lower_range[dir], &upper_range[dir],
                         dir + 1, NULL, system_name) < 0)
    {
      CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
                  "IOUtil_2DPlanes: Could not get ranges for %c-direction "
                  "of coordinate system '%s'",
                  'x' + dir, system_name);
    }
  }

  /* now set the slice center for each line
     according to origin_index[] or origin_phys[] */
  for (dir = 0; dir < num_dims; dir++)
  {
    if (origin_index && origin_index[dir] >= 0)
    {
      slice_center[dir] = origin_index[dir];
    }
    else if (lower_range[dir] > origin_phys[dir] || 
             upper_range[dir] < origin_phys[dir])
    {
      CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING,
                  "IOUtil_2DPlanes: Coordinate in %c-direction (%f, %f) "
                  "doesn't contain %f",
                  'x' + dir, (double) lower_range[dir],
                  (double) upper_range[dir], (double) origin_phys[dir]);
      slice_center[dir] = 0;
    }
    else
    {
      /* Find index for first point above the chosen coordinate */
      slice_center[dir] = ceil ((origin_phys[dir] - lower_range[dir]) /
                                GH->cctk_delta_space[dir] - 1e-6);

#ifdef DEBUG_IOUTIL
      printf("sp2xyz for planes perpendicular to"
             " %d-direction is %d\n",
             dir,slice_center[dir]);
#endif
    }
  }

  /* free allocated resources */
  free (lower_range);

  return (0);
}


 /*@@
   @routine    IOUtil_PrintTimings
   @date       Wed Jun 28 2000
   @author     Thomas Radke
   @desc
               Gets the timing information for the given timers and prints it
               as INFO messages to screen.
   @enddesc
   @history
   @endhistory
   @var        description
   @vdesc      description of the timers
   @vtype      const char *
   @vio        in
   @endvar
   @var        ntimers
   @vdesc      number of timers passed in
   @vtype      int
   @vio        in
   @endvar
   @var        timers
   @vdesc      array of timers
   @vtype      int [ntimers]
   @vio        in
   @endvar
   @var        timer_descriptions
   @vdesc      array of timer descriptions
   @vtype      const char * [ntimers]
   @vio        in
   @endvar
@@*/

void IOUtil_PrintTimings (const char *description,
                          int ntimers,
                          int timers[/* ntimers */],
                          const char *timer_descriptions[/* ntimers */])
{
  int i, j;
  cTimerData *info;


  info = CCTK_TimerCreateData ();
  if (info)
  {
    CCTK_INFO (description);

    for (i = 0; i < info->n_vals; i++)
    {
      for (j = 0; j < ntimers; j++)
      {
        CCTK_TimerI (timers[j], info);
        if (j == 0)
        {
          CCTK_VInfo (CCTK_THORNSTRING, "  %s:", info->vals[i].heading);
        }
        switch (info->vals[i].type)
        {
          case val_int:
            CCTK_VInfo (CCTK_THORNSTRING, "    %s %5d %s",
                        timer_descriptions [j],
                        info->vals[i].val.i, info->vals[i].units);
            break;

          case val_long:
            CCTK_VInfo (CCTK_THORNSTRING, "    %s %5d %s",
                        timer_descriptions [j],
                        (int) info->vals[i].val.l, info->vals[i].units);
            break;

          case val_double:
            CCTK_VInfo (CCTK_THORNSTRING, "    %s %5.1f %s",
                        timer_descriptions [j],
                        info->vals[i].val.d, info->vals[i].units);
            break;

          default:
            CCTK_WARN(1, "Unknown data type for timer info");
            break;
        }
      }
    }
    CCTK_INFO ("-----------------------------------------\n");
    CCTK_TimerDestroyData (info);
  }
  else
  {
    CCTK_WARN (1, "Couldn't create timer info structure ! "
                  "No timing output available.");
  }
}