path: root/src/utils.cc

#include <stdio.h>
#include <string.h>
#include <math.h>
#include <string>
#include <map>
#include <sys/stat.h>
#include <errno.h>
#include <iostream>

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

#ifdef HAVE_CAPABILITY_HDF5
// We currently support the HDF5 1.6 API (and when using 1.8 the
// compatibility mode introduced by H5_USE_16_API).  Several machines
// in SimFactory use HDF5 1.6, so we cannot drop support for it.  It
// seems it is hard to support both the 1.6 and 1.8 API
// simultaneously; for example H5Fopen takes a different number of
// arguments in the two versions.
#define H5_USE_16_API
#include <hdf5.h>
#endif

#include "utils.hh"
#include "integrate.hh"

using namespace std;

////////////////////////////////////////////////////////////////////////
// File manipulation
////////////////////////////////////////////////////////////////////////

FILE *Multipole_OpenOutputFile(CCTK_ARGUMENTS, const string &name)
{
  DECLARE_CCTK_ARGUMENTS;
  DECLARE_CCTK_PARAMETERS;
  
  bool first_time = cctk_iteration == 0;
  const char *mode = first_time ? "w" : "a";

  string output_name(string(out_dir) + "/" + name);

  FILE *fp = fopen(output_name.c_str(), mode);

  if (fp == 0)
  {
    CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, (string("Could not open output file ") + output_name).c_str());
  }

  return fp;
}

////////////////////////////////////////////////////////////////////////
// Unused
////////////////////////////////////////////////////////////////////////

void Multipole_OutputArray(CCTK_ARGUMENTS, FILE *f, int array_size, 
                           CCTK_REAL const th[], CCTK_REAL const ph[], 
                           CCTK_REAL const xs[], CCTK_REAL const ys[], CCTK_REAL const zs[],
                           CCTK_REAL const data[])
{
  DECLARE_CCTK_PARAMETERS;
  DECLARE_CCTK_ARGUMENTS;

  CCTK_REAL last_ph = ph[0];

  for (int i = 0; i < array_size; i++)
  {
    if (ph[i] != last_ph) // Separate blocks for gnuplot
      fprintf(f, "\n");
    fprintf(f, "%f %f %f %f %f %f %.19g\n", cctk_time, th[i], ph[i], xs[i], ys[i], zs[i], data[i]);
    last_ph = ph[i];
  }
}


void Multipole_OutputArrayToFile(CCTK_ARGUMENTS, const string &name, int array_size,
                                 CCTK_REAL const th[], CCTK_REAL const ph[],
                                 CCTK_REAL const xs[], CCTK_REAL const ys[], CCTK_REAL const zs[],
                                 CCTK_REAL const data[])
{
  DECLARE_CCTK_ARGUMENTS;

  if (FILE *fp = Multipole_OpenOutputFile(CCTK_PASS_CTOC, name))
  {
    Multipole_OutputArray(CCTK_PASS_CTOC, fp, array_size, th, ph, xs, ys, zs, data);
    fclose(fp);
  }
}

////////////////////////////////////////////////////////////////////////
// Misc
////////////////////////////////////////////////////////////////////////

void Multipole_Output1D(CCTK_ARGUMENTS, const string &name, int array_size,
                        CCTK_REAL const th[], CCTK_REAL const ph[], mp_coord coord,
                        CCTK_REAL const data[])
{
  DECLARE_CCTK_ARGUMENTS;
  DECLARE_CCTK_PARAMETERS;

  if (FILE *f = Multipole_OpenOutputFile(CCTK_PASS_CTOC, name))
  {
    fprintf(f, "\"Time = %.19g\n", cctk_time);

    if (coord == mp_theta)
    {
      for (int i = 0; i <= ntheta; i++)
      {
        int idx = Multipole_Index(i, 0, ntheta);
        fprintf(f, "%f %.19g\n", th[idx], data[idx]);
      }
    }
    else if (coord == mp_phi)
    {
      for (int i = 0; i <= nphi; i++)
      {
        int idx = Multipole_Index(ntheta / 4, i, ntheta);
        fprintf(f, "%f %.19g\n", ph[idx], data[idx]);
      }
    }
    fprintf(f, "\n\n");
    fclose(f);
  }
}

////////////////////////////////////////////////////////////////////////
// Complex IO
////////////////////////////////////////////////////////////////////////

void Multipole_OutputComplex(CCTK_ARGUMENTS, FILE *fp, CCTK_REAL redata, CCTK_REAL imdata)
{
  DECLARE_CCTK_PARAMETERS;
  DECLARE_CCTK_ARGUMENTS;
  fprintf(fp, "%f %.19g %.19g\n", cctk_time, redata, imdata);
}

void Multipole_OutputComplexToFile(CCTK_ARGUMENTS, const string &name, CCTK_REAL redata, CCTK_REAL imdata)
{
  DECLARE_CCTK_ARGUMENTS;

  if (FILE *fp = Multipole_OpenOutputFile(CCTK_PASS_CTOC, name))
  {
    Multipole_OutputComplex(CCTK_PASS_CTOC, fp, redata, imdata);
    fclose(fp);
  }
}

////////////////////////////////////////////////////////////////////////
// HDF5 complex output
////////////////////////////////////////////////////////////////////////

#ifdef HAVE_CAPABILITY_HDF5

static bool file_exists(const string &name)
{
  struct stat sts;
  return !(stat(name.c_str(), &sts) == -1 && errno == ENOENT);
}

bool dataset_exists(hid_t file, const string &dataset_name)
{
  // To test whether a dataset exists, the recommended way in API 1.6
  // is to use H5Gget_objinfo, but this prints an error to stderr if
  // the dataset does not exist.  We explicitly avoid this by wrapping
  // the call in H5E_BEGIN_TRY/H5E_END_TRY statements.  In 1.8,
  // H5Gget_objinfo is deprecated, and H5Lexists does the job.  See
  // http://www.mail-archive.com/hdf-forum@hdfgroup.org/msg00125.html

  #if 1
  bool exists;
  H5E_BEGIN_TRY
  {
    exists = H5Gget_objinfo(file, dataset_name.c_str(), 1, NULL) >= 0;
  }
  H5E_END_TRY;
  return exists;
  #else
  return H5Lexists(file, dataset_name.c_str(), H5P_DEFAULT);
  #endif
}

void Multipole_OutputComplexToH5File(CCTK_ARGUMENTS, const string &basename, const string &datasetname,
                                     CCTK_REAL redata, CCTK_REAL imdata)
{
  DECLARE_CCTK_ARGUMENTS;
  DECLARE_CCTK_PARAMETERS;

  string output_name = out_dir + string("/") + basename;
  int status = 0;
  static map<string,bool> checked; // Has the given file been checked
                                   // for truncation? map<*,bool>
                                   // defaults to false

  hid_t file;

  if (!file_exists(output_name) || (!checked[output_name] && IO_TruncateOutputFiles(cctkGH)))
  {
    file = H5Fcreate(output_name.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
  }
  else
  {
    file = H5Fopen(output_name.c_str(), H5F_ACC_RDWR, H5P_DEFAULT);
  }

  checked[output_name] = true;

  hid_t dataset = -1;

  if (dataset_exists(file, datasetname))
  {
    dataset = H5Dopen(file, datasetname.c_str());
  }
  else
  {
    hsize_t dims[2] = {0,3};
    hsize_t maxdims[2] = {H5S_UNLIMITED, 3};
    hid_t dataspace = H5Screate_simple(2, dims, maxdims);

    hid_t cparms = -1;
    hsize_t chunk_dims[2] = {hdf5_chunk_size,3};
    cparms = H5Pcreate (H5P_DATASET_CREATE);
    status = H5Pset_chunk(cparms, 2, chunk_dims);

    dataset = H5Dcreate(file, datasetname.c_str(), H5T_NATIVE_DOUBLE, dataspace, cparms);
    H5Pclose(cparms);
  }

  hid_t filespace = H5Dget_space (dataset);
  hsize_t filedims[2];
  hsize_t maxdims[2];
  status = H5Sget_simple_extent_dims(filespace, filedims, maxdims);

  filedims[0] += 1;
  hsize_t size[2] = {filedims[0], filedims[1]};
  status = H5Dextend (dataset, size);
  H5Sclose(filespace);

  /* Select a hyperslab  */
  hsize_t offset[2] = {filedims[0]-1, 0};
  hsize_t dims2[2] = {1,3};
  filespace = H5Dget_space (dataset);
  status = H5Sselect_hyperslab (filespace, H5S_SELECT_SET, offset, NULL,
                                dims2, NULL);

  CCTK_REAL data[] = {cctk_time, redata, imdata};

  hid_t memdataspace = H5Screate_simple(2, dims2, NULL);

  /* Write the data to the hyperslab  */
  status = H5Dwrite (dataset, H5T_NATIVE_DOUBLE, memdataspace, filespace,
                     H5P_DEFAULT, data);

  H5Dclose(dataset);
  H5Sclose(filespace);
  H5Sclose(memdataspace);

  status = H5Fclose(file);
}

#else

void Multipole_OutputComplexToH5File(CCTK_ARGUMENTS, const string &basename, const string &datasetname,
                                     CCTK_REAL redata, CCTK_REAL imdata)
{
  CCTK_WARN(0,"HDF5 output has been requested but Cactus has been compiled without HDF5 support");
}

#endif

////////////////////////////////////////////////////////////////////////
// Coordinates
////////////////////////////////////////////////////////////////////////

void Multipole_CoordSetup(int ntheta, int nphi, 
                          CCTK_REAL xhat[], CCTK_REAL yhat[], 
                          CCTK_REAL zhat[], CCTK_REAL th[], 
                          CCTK_REAL ph[])
{
  const CCTK_REAL PI = acos(-1.0);

  for (int it = 0; it <= ntheta; it++)
  {
    for (int ip = 0; ip <= nphi; ip++)
    {
      const int i = Multipole_Index(it, ip, ntheta);

      th[i] = it * PI / (ntheta);
      ph[i] = ip * 2 * PI / nphi;
      xhat[i] = cos(ph[i])*sin(th[i]);
      yhat[i] = sin(ph[i])*sin(th[i]);
      zhat[i] = cos(th[i]);
    }
  }
}

void Multipole_ScaleCartesian(int ntheta, int nphi, CCTK_REAL r,
                              CCTK_REAL const xhat[], CCTK_REAL const yhat[], CCTK_REAL const zhat[],
                              CCTK_REAL x[], CCTK_REAL y[], CCTK_REAL z[])
{
  for (int it = 0; it <= ntheta; it++)
  {
    for (int ip = 0; ip <= nphi; ip++)
    {
      const int i = Multipole_Index(it, ip, ntheta);

      x[i] = r * xhat[i];
      y[i] = r * yhat[i];
      z[i] = r * zhat[i];
    }
  }
}

////////////////////////////////////////////////////////////////////////
// Integration
////////////////////////////////////////////////////////////////////////

void Multipole_Integrate(int array_size, int nthetap,
    CCTK_REAL const array1r[], CCTK_REAL const array1i[],
    CCTK_REAL const array2r[], CCTK_REAL const array2i[],
    CCTK_REAL const th[], CCTK_REAL const ph[], 
    CCTK_REAL *outre, CCTK_REAL *outim)
{
  DECLARE_CCTK_PARAMETERS

  int il = Multipole_Index(0,0,ntheta);
  int iu = Multipole_Index(1,0,ntheta);
  CCTK_REAL dth = th[iu] - th[il];
  iu = Multipole_Index(0,1,ntheta);
  CCTK_REAL dph = ph[iu] - ph[il];

  if (CCTK_Equals(integration_method, "midpoint"))
  {
    CCTK_REAL tempr = 0.0;
    CCTK_REAL tempi = 0.0;
      
    for (int i=0; i < array_size; i++) {
      // the below calculations take the integral of conj(array1)*array2
      tempr += ( array1r[i]*array2r[i] + array1i[i]*array2i[i] )
        *sin(th[i])*dth*dph;
      tempi += ( array1r[i]*array2i[i] - array1i[i]*array2r[i] )
        *sin(th[i])*dth*dph;
      
      *outre = tempr;
      *outim = tempi;
    } 
  }
  else if (CCTK_Equals(integration_method, "simpson"))
  {
    static CCTK_REAL *fr = 0;
    static CCTK_REAL *fi = 0;
    static bool allocated_memory = false;

    // Construct an array for the real integrand
    if (!allocated_memory)
    {
      fr = new CCTK_REAL[array_size];
      fi = new CCTK_REAL[array_size];
      allocated_memory = true;
    }
    
    for (int i = 0; i < array_size; i++)
    {
      fr[i] = (array1r[i] * array2r[i] + 
               array1i[i] * array2i[i] ) * sin(th[i]);
      fi[i] = (array1r[i] * array2i[i] - 
               array1i[i] * array2r[i] ) * sin(th[i]);
    }
    
    *outre = Simpson2DIntegral(fr, ntheta, nphi, dth, dph);
    *outim = Simpson2DIntegral(fi, ntheta, nphi, dth, dph);
  }
}