path: root/src/WaveToy.c

 /*@@
   @file      WaveToy.c
   @date      
   @author    Tom Goodale
   @desc 
              Evolution routines for the wave equation solver
   @enddesc 
   @version $Header$
 @@*/

#include <assert.h>

#if 0
#include <sys/time.h>
#include <sys/resource.h>
#include <unistd.h>
#endif

#include "cctk.h" 
#include "cctk_Arguments.h"
#include "cctk_Parameters.h"

#include "Symmetry.h"

static const char *rcsid = "$Header$";

CCTK_FILEVERSION(CactusWave_WaveToyC_WaveToy_c);

void WaveToyC_Boundaries(CCTK_ARGUMENTS);
void WaveToyC_Evolution(CCTK_ARGUMENTS);

#define DEBUG 0



static inline int min (int const a, int const b)
{
  return a<b ? a : b;
}

static inline int max (int const a, int const b)
{
  return a>b ? a : b;
}
    
#if 0
static inline long long timeval2longlong (const struct timeval * const tv)
{
  return tv->tv_sec * 1000000LL + tv->tv_usec;
}

static long long currenttime (void)
{
  struct rusage ru;
  getrusage (RUSAGE_SELF, &ru);
  return timeval2longlong (ru.ru_utime) + timeval2longlong (ru.ru_stime);
}
#endif

static long long currenttime (void)
{
  unsigned long eax, edx;
  asm volatile ("rdtsc" : "=a" (eax), "=d" (edx));
  return (unsigned long long)edx << 32 | (unsigned long long)eax;
}



 /*@@
   @routine    WaveToyC_Evolution
   @date       
   @author     Tom Goodale
   @desc 
               Evolution for the wave equation
   @enddesc 
   @calls      CCTK_SyncGroup, WaveToyC_Boundaries
   @calledby   
   @history 
   @endhistory 
@@*/

void  WaveToyC_Evolution(CCTK_ARGUMENTS)
{
  DECLARE_CCTK_ARGUMENTS;
  DECLARE_CCTK_PARAMETERS;
      
  int i,j,k;
  int index;
  int istart, jstart, kstart, iend, jend, kend;
  CCTK_REAL dx,dy,dz,dt,dx2,dy2,dz2,dt2;
  CCTK_REAL dt2dx2i,dt2dy2i,dt2dz2i;

  CCTK_REAL factor;

  /* Set up shorthands */
  dx = CCTK_DELTA_SPACE(0);
  dy = CCTK_DELTA_SPACE(1);
  dz = CCTK_DELTA_SPACE(2);
  dt = CCTK_DELTA_TIME;

  dx2 = dx*dx;
  dy2 = dy*dy;
  dz2 = dz*dz;
  dt2 = dt*dt;

  dt2dx2i = dt2/dx2;
  dt2dy2i = dt2/dy2;
  dt2dz2i = dt2/dz2;

  istart = 1;
  jstart = 1;
  kstart = 1;
  
  iend = cctk_lsh[0]-1;
  jend = cctk_lsh[1]-1;
  kend = cctk_lsh[2]-1;

  /* Do the evolution */
  factor = 2 * (1 - (dt2dx2i + dt2dy2i + dt2dz2i));

  {
    /* user-defined values */
    const char * const loopthorn = CCTK_THORNSTRING;
    const char * const loopfile = __FILE__;
    int          const loopline = __LINE__;
    const char * const loopname = "Evolution";
    
    int const dim = 3;
    
    int const imin = istart;
    int const jmin = jstart;
    int const kmin = kstart;
    int const imax = iend;
    int const jmax = jend;
    int const kmax = kend;
    int const ishape = cctk_lsh[0];
    int const jshape = cctk_lsh[1];
    int const kshape = cctk_lsh[2];
    
    /* tuning parameters */
    if (DEBUG) printf ("level1_stride = [%d,%d,%d]\n", level1_stride[0], level1_stride[1], level1_stride[2]);
    if (DEBUG) printf ("level2_stride = [%d,%d,%d]\n", level2_stride[0], level2_stride[1], level2_stride[2]);
    
    int const i1step = level1_stride[0]<0 ? imax-imin : level1_stride[0];
    int const j1step = level1_stride[1]<0 ? jmax-jmin : level1_stride[1];
    int const k1step = level1_stride[2]<0 ? kmax-kmin : level1_stride[2];
    
    int const i2step = level2_stride[0]<0 ? imax-imin : level2_stride[0];
    int const j2step = level2_stride[1]<0 ? jmax-jmin : level2_stride[1];
    int const k2step = level2_stride[2]<0 ? kmax-kmin : level2_stride[2];
    
    /* internal stuff */
    int rep;
    for (rep = 0; rep < repetitions; ++rep) {
    
    long long l0time  = 0, l1time  = 0, l2time  = 0;
    long long l0count = 0, l1count = 0, l2count = 0;
    
    /* level 2 loop */
    l2time = currenttime() - l2time;
    ++ l2count;
    
    int const i2min = imin;
    int const j2min = jmin;
    int const k2min = kmin;
    int const i2max = imax;
    int const j2max = jmax;
    int const k2max = kmax;
    
    int i2, j2, k2;
    
    for (k2 = k2min; k2 < k2max; k2 += k2step) {
    for (j2 = j2min; j2 < j2max; j2 += j2step) {
    for (i2 = i2min; i2 < i2max; i2 += i2step) {
    
    /* level 1 loop */
    l1time = currenttime() - l1time;
    ++ l1count;
    
    int const i1min = i2;
    int const j1min = j2;
    int const k1min = k2;
    int const i1max = min (imax, i2 + i2step);
    int const j1max = min (jmax, j2 + j2step);
    int const k1max = min (kmax, k2 + k2step);
    
    int i1, j1, k1;
    
    for (k1 = k1min; k1 < k1max; k1 += k1step) {
    for (j1 = j1min; j1 < j1max; j1 += j1step) {
    for (i1 = i1min; i1 < i1max; i1 += i1step) {
    
    /* level 0 loop */
    l0time = currenttime() - l0time;
    ++ l0count;
    
    int const i0min = i1;
    int const j0min = j1;
    int const k0min = k1;
    int const i0max = min (imax, i1 + i1step);
    int const j0max = min (jmax, j1 + j1step);
    int const k0max = min (kmax, k1 + k1step);
    
    int i0, j0, k0;
    
    for (k0 = k0min; k0 < k0max; ++k0) {
    for (j0 = j0min; j0 < j0max; ++j0) {
    for (i0 = i0min; i0 < i0max; ++i0) {
    
    int const i = i0;
    int const j = j0;
    int const k = k0;
    
    int const di = 1;
    int const dj = ishape;
    int const dk = ishape * jshape;
    
    int const index = i *di + j * dj + k * dk;
    
    {
        phi[index] =   factor * phi_p[index] - phi_p_p[index] 
                     + (phi_p[index+di] + phi_p[index-di]) * dt2dx2i
                     + (phi_p[index+dj] + phi_p[index-dj]) * dt2dy2i
                     + (phi_p[index+dk] + phi_p[index-dk]) * dt2dz2i;
    } 
    
    } } } /* end level 0 loop */
    l0time = currenttime() - l0time;
    
    } } } /* end level 1 loop */
    l1time = currenttime() - l1time;
    
    } } } /* end level 2 loop */
    l2time = currenttime() - l2time;
    
    if (DEBUG)
      CCTK_VInfo (CCTK_THORNSTRING,
                  "Timing information for loop %s::%s\n"
                  "   (file \"%s\", line %d):\n"
                  "   Level 0: shape %d %d %d, count %lld, time %lld\n"
                  "   Level 1: shape %d %d %d, count %lld, time %lld\n"
                  "   Level 2: shape %d %d %d, count %lld, time %lld",
                  loopthorn, loopname, loopfile, loopline,
                  i1step   , j1step   , k1step   , l0count, l0time,
                  i2step   , j2step   , k2step   , l1count, l1time,
                  imax-imin, jmax-jmin, kmax-kmin, l2count, l2time);
    
    } /* end repetition loop */
    
  }
    
  return;
}


 /*@@
   @routine    WaveToyC_Boundaries
   @date       
   @author     Tom Goodale
   @desc 
               Boundary conditions for the wave equation
   @enddesc 
   @history 
   @endhistory 
@@*/

void WaveToyC_Boundaries(CCTK_ARGUMENTS)
{
  DECLARE_CCTK_ARGUMENTS;
  DECLARE_CCTK_PARAMETERS;

  int ierr=0;

  ierr = CartSymGN(cctkGH,"wavetoy::scalarevolve");

  if (CCTK_EQUALS(bound,"flat") || CCTK_EQUALS(bound,"static") || 
      CCTK_EQUALS(bound,"radiation") || CCTK_EQUALS(bound,"robin") || 
      CCTK_EQUALS(bound,"zero") || CCTK_EQUALS(bound,"none"))
  {
    /* Uses all default arguments, so invalid table handle -1 can be passed */
    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, -1, "wavetoy::phi", 
				   bound);
  }

  if (ierr < 0) 
  {
    CCTK_VWarn(0,__LINE__,__FILE__,CCTK_THORNSTRING,
	       "WaveToyC_Boundaries: Error selecting boundary "
	       "condition %s",bound);
  }

  return;
}