path: root/src/RobinBoundary.c

 /*@@
   @file      RobinBoundary.c
   @date      July 6th 2000
   @author    Gabrielle Allen, Gerd Lanfermann
   @desc
              Routines for Robin boundary conditions
   @enddesc
   @history
   @hdate     Tue 10 Apr 2001
   @hauthor   Thomas Radke
   @hdesc     BC routines generalized for applying to arbitrary CCTK data types
   @endhistory
   @version   $Id$
 @@*/

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

#include "cctk.h"
#include "cctk_Parameters.h"
#include "cctk_FortranString.h"

#include "Symmetry.h"
#include "Boundary.h"

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


/********************************************************************
 ********************    External Routines   ************************
 ********************************************************************/
/* prototypes for external C routines are declared in header Boundary.h
   here only follow the fortran wrapper prototypes */
void CCTK_FCALL CCTK_FNAME (BndRobinGI)
                           (int *ierr,
                            cGH *GH,
                            int stencil[],
                            const CCTK_REAL *finf,
                            const int *npow,
                            const int *gi);
void CCTK_FCALL CCTK_FNAME (BndRobinGN)
                           (int *ierr,
                            cGH *GH,
                            int stencil[],
                            const CCTK_REAL *finf,
                            const int *npow,
                            ONE_FORTSTRING_ARG);
void CCTK_FCALL CCTK_FNAME (BndRobinVI)
                           (int *ierr,
                            cGH *GH,
                            int stencil[],
                            const CCTK_REAL  *finf,
                            const int *npow,
                            const int *vi);
void CCTK_FCALL CCTK_FNAME (BndRobinVN)
                           (int *ierr,
                            cGH *GH,
                            int stencil[],
                            const CCTK_REAL *finf,
                            const int *npow,
                            ONE_FORTSTRING_ARG);


/********************************************************************
 ********************    Internal Routines   ************************
 ********************************************************************/
static int ApplyBndRobin (cGH *GH,
                          int stencil[],
                          CCTK_REAL finf,
                          int npow,
                          int first_var,
                          int num_vars);


/*@@
   @routine    BndRobinGI
   @date       Tue Jul 18 18:08:28 2000
   @author     Gerd Lanfermann
   @desc
               Apply Robin boundary conditions by group index
   @enddesc
   @calls      ApplyBndRobin

   @var        GH
   @vdesc      Pointer to CCTK grid hierarchy
   @vtype      cGH *
   @vio        in
   @endvar
   @var        stencil
   @vdesc      stencil width array
   @vtype      int [ dimension of group ]
   @vio        in
   @endvar
   @var        finf
   @vdesc      value of f at infimum
   @vtype      CCTK_REAL
   @vio        in
   @endvar
   @var        npow
   @vdesc      power of decay rate
   @vtype      int
   @vio        in
   @endvar
   @var        gi
   @vdesc      index of group to apply boundary conditions to
   @vtype      int
   @vio        in
   @endvar

   @returntype int
   @returndesc
               return code of @seeroutine ApplyBndRobin <BR>
               -1 if invalid group index was given
   @endreturndesc
@@*/
int BndRobinGI (cGH *GH,
                int stencil[],
                CCTK_REAL finf,
                int npow,
                int gi)
{
  int first_vi, retval;


  first_vi = CCTK_FirstVarIndexI (gi);
  if (first_vi >= 0)
  {
    retval = ApplyBndRobin (GH, stencil, finf, npow, first_vi,
                            CCTK_NumVarsInGroupI (gi));
  }
  else
  {
    CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING,
                "Invalid group index %d in BndFlatGI", gi);
    retval = -1;
  }

  return (retval);
}

void CCTK_FCALL CCTK_FNAME (BndRobinGI)
                           (int *ierr,
                            cGH *GH,
                            int stencil[],
                            const CCTK_REAL *finf,
                            const int *npow,
                            const int *gi)
{
  *ierr = BndRobinGI (GH, stencil, *finf, *npow, *gi);
}


/*@@
   @routine    BndRobinGN
   @date       Tue Jul 18 18:08:28 2000
   @author     Gerd Lanfermann
   @desc
               Apply Robin boundary conditions by group name
   @enddesc
   @calls      BndRobinGI

   @var        GH
   @vdesc      Pointer to CCTK grid hierarchy
   @vtype      cGH *
   @vio        in
   @endvar
   @var        stencil
   @vdesc      stencil width array
   @vtype      int [ dimension of group ]
   @vio        in
   @endvar
   @var        finf
   @vdesc      value of f at infimum
   @vtype      CCTK_REAL
   @vio        in
   @endvar
   @var        npow
   @vdesc      power of decay rate
   @vtype      int
   @vio        in
   @endvar
   @var        gn
   @vdesc      name of group to apply boundary conditions to
   @vtype      char []
   @vio        in
   @endvar

   @returntype int
   @returndesc
               return code of @seeroutine BndRobinGI <BR>
               -1 if invalid group index was given
   @endreturndesc
@@*/
int BndRobinGN (cGH *GH,
                int stencil[],
                CCTK_REAL finf,
                int npow,
                char gn[])
{
  int gi, retval;


  gi = CCTK_GroupIndex (gn);
  if (gi >= 0)
  {
    retval = BndRobinGI (GH, stencil, finf, npow, gi);
  }
  else
  {
    CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING,
                "Invalid group name '%s' in BndRobinGN", gn);
    retval = -1;
  }

  return (retval);
}

void CCTK_FCALL CCTK_FNAME (BndRobinGN)
                           (int *ierr,
                            cGH *GH,
                            int stencil[],
                            const CCTK_REAL *finf,
                            const int *npow,
                            ONE_FORTSTRING_ARG)
{
  ONE_FORTSTRING_CREATE (gn)
  *ierr = BndRobinVN (GH, stencil, *finf, *npow, gn);
  free (gn);
}


/*@@
   @routine    BndRobinVI
   @date       Tue Jul 18 18:08:28 2000
   @author     Gerd Lanfermann
   @desc
               Apply Robin boundary conditions by variable index
   @enddesc
   @calls      ApplyBndRobin

   @var        GH
   @vdesc      Pointer to CCTK grid hierarchy
   @vtype      cGH *
   @vio        in
   @endvar
   @var        stencil
   @vdesc      stencil width array
   @vtype      int [ dimension of variable ]
   @vio        in
   @endvar
   @var        finf
   @vdesc      value of f at infimum
   @vtype      CCTK_REAL
   @vio        in
   @endvar
   @var        npow
   @vdesc      power of decay rate
   @vtype      int
   @vio        in
   @endvar
   @var        vi
   @vdesc      index of variable to apply boundary conditions to
   @vtype      char []
   @vio        in
   @endvar

   @returntype int
   @returndesc
               return code of @seeroutine ApplyBndRobin <BR>
               -1 if invalid group index was given
   @endreturndesc
@@*/
int BndRobinVI (cGH *GH,
                int stencil[],
                CCTK_REAL finf,
                int npow,
                int vi)
{
  int retval;


  if (vi >= 0 && vi < CCTK_NumVars ())
  {
    retval = ApplyBndRobin (GH, stencil, finf, npow, vi, 1);
  }
  else
  {
    CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING,
                "BndRobinVI: Invalid variable index %d", vi);
    retval = -1;
  }

  return (retval);
}

void CCTK_FCALL CCTK_FNAME (BndRobinVI)
                           (int *ierr,
                            cGH *GH,
                            int stencil[],
                            const CCTK_REAL  *finf,
                            const int *npow,
                            const int *vi)
{
  *ierr = BndRobinVI (GH, stencil, *finf, *npow, *vi);
}


/*@@
   @routine    BndRobinVN
   @date       Tue Jul 18 18:08:28 2000
   @author     Gerd Lanfermann
   @desc
               Apply Robin boundary conditions by variable name
   @enddesc
   @calls      BndRobinVI

   @var        GH
   @vdesc      Pointer to CCTK grid hierarchy
   @vtype      cGH *
   @vio        in
   @endvar
   @var        stencil
   @vdesc      stencil width array
   @vtype      int [ dimension of variable ]
   @vio        in
   @endvar
   @var        finf
   @vdesc      value of f at infimum
   @vtype      CCTK_REAL
   @vio        in
   @endvar
   @var        npow
   @vdesc      power of decay rate
   @vtype      int
   @vio        in
   @endvar
   @var        vn
   @vdesc      name of variable to apply boundary conditions to
   @vtype      char []
   @vio        in
   @endvar

   @returntype int
   @returndesc
               return code of @seeroutine BndRobinVI <BR>
               -1 if invalid group index was given
   @endreturndesc
@@*/
int BndRobinVN (cGH *GH,
                int stencil[],
                CCTK_REAL finf,
                int npow,
                char vn[])
{
  int vi, retval;


  vi = CCTK_VarIndex (vn);
  if (vi >= 0)
  {
    retval = BndRobinVI (GH, stencil, finf, npow, vi);
  }
  else
  {
    CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING,
                "Invalid variable name '%s' in BndRobinVN", vn);
    retval = -1;
  }

  return (retval);
}

void CCTK_FCALL CCTK_FNAME (BndRobinVN)
                           (int *ierr,
                            cGH *GH,
                            int stencil[],
                            const CCTK_REAL *finf,
                            const int *npow,
                            ONE_FORTSTRING_ARG)
{
  ONE_FORTSTRING_CREATE (vn)
  *ierr = BndRobinVN (GH, stencil, *finf, *npow, vn);
  free (vn);
}


/********************************************************************
 *********************     Local Routines   *************************
 ********************************************************************/

/* maximum dimension we can deal with */
#define MAXDIM  3

/* macro to compute the linear index of a 3D point */
#define INDEX_3D(lsh, i, j, k)      ((i) + (lsh)[0]*((j) + (lsh)[1]*(k)))

/*@@
   @routine    ROBIN_BOUNDARY_3D_TYPED
   @date       Tue 10 April 2001
   @author     Thomas Radke
   @desc
               Macro to apply Robin boundary conditions to a 3D variable
               of given datatype
   @enddesc

   @var        doBC
   @vdesc      flag telling whether to apply boundary conditions or not
   @vtype      int
   @vio        in
   @endvar
   @var        iend, jend
   @vdesc      upper ranges for the loopers
   @vtype      int
   @vio        in
   @endvar
   @var        ii, jj, kk
   @vdesc      indices of the current grid point
   @vtype      int
   @vio        in
   @endvar
   @var        offset
   @vdesc      offset to the next neighbouring grid point
   @vtype      int
   @vio        in
   @endvar
   @var        cctk_type
   @vdesc      CCTK datatype of the variable
   @vtype      <cctk_type>
   @vio        in
   @endvar
@@*/
#define ROBIN_BOUNDARY_3D_TYPED(doBC,                                         \
                                iend, jend,                                   \
                                ii, jj, kk,                                   \
                                offset,                                       \
                                cctk_type)                                    \
{                                                                             \
  if (doBC)                                                                   \
  {                                                                           \
    for (j = 0; j < jend; j++)                                                \
    {                                                                         \
      for (i = 0; i < iend; i++)                                              \
      {                                                                       \
        int _i = CCTK_GFINDEX3D (GH, ii, jj, kk);                             \
        cctk_type *_var = (cctk_type *) (GH->data[var][timelvl]) + _i;        \
        const CCTK_REAL *_r = radius + _i;                                    \
        /* the native C compiler on DEC Alphas has problems with automatic    \
           type-casting of integers to doubles for inlined math routines */   \
        const double dpow = (double) npow;                                    \
                                                                              \
                                                                              \
        _var[0] = (cctk_type) (                                               \
                   (  1 * pow (_r[3*offset], dpow) * (_var[3*offset] - finf)  \
                    - 3 * pow (_r[2*offset], dpow) * (_var[2*offset] - finf)  \
                    + 3 * pow (_r[1*offset], dpow) * (_var[1*offset] - finf)  \
                   ) /    pow (_r[0], dpow)                                   \
                  + finf);                                                    \
      }                                                                       \
    }                                                                         \
  }                                                                           \
}


/*@@
   @routine    ROBIN_BOUNDARY
   @date       Sat 20 Jan 2001
   @author     Thomas Radke
   @desc
               Macro to apply Robin boundary conditions to a variable
               of a given datatype in all directions
               Currently it is limited up to 3D variables only.
   @enddesc

   @var        cctk_type
   @vdesc      CCTK datatype of the variable
   @vtype      <cctk_type>
   @vio        in
   @endvar
@@*/
#define ROBIN_BOUNDARY(cctk_type)                                             \
{                                                                             \
  /* check the dimensionality */                                              \
  if (gdim != 3)                                                              \
  {                                                                           \
    CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,                      \
                "ApplyBndRobin: variable dimension of %d not supported",      \
                gdim);                                                        \
    return (-5);                                                              \
  }                                                                           \
                                                                              \
  /* now set the boundaries face by face */                                   \
  /* lower x */                                                               \
  ROBIN_BOUNDARY_3D_TYPED (doBC[0], lsh[1], lsh[2], 0, i, j,                  \
                           offset[0], cctk_type);                             \
  /* upper x */                                                               \
  ROBIN_BOUNDARY_3D_TYPED (doBC[1], lsh[1], lsh[2], lsh[0]-1, i, j,           \
                           -offset[0], cctk_type);                            \
  /* lower y */                                                               \
  ROBIN_BOUNDARY_3D_TYPED (doBC[2], lsh[0], lsh[2], i, 0, j,                  \
                           offset[1], cctk_type);                             \
  /* upper y */                                                               \
  ROBIN_BOUNDARY_3D_TYPED (doBC[3], lsh[0], lsh[2], i, lsh[1]-1, j,           \
                           -offset[1], cctk_type);                            \
  /* lower z */                                                               \
  ROBIN_BOUNDARY_3D_TYPED (doBC[4], lsh[0], lsh[1], i, j, 0,                  \
                           offset[2], cctk_type);                             \
  /* upper z */                                                               \
  ROBIN_BOUNDARY_3D_TYPED (doBC[5], lsh[0], lsh[1], i, j, lsh[2]-1,           \
                           -offset[2], cctk_type);                            \
}


/*@@
   @routine    ApplyBndRobin
   @date       Tue Jul 18 18:08:28 2000
   @author     Gerd Lanfermann
   @desc
               Apply Robin boundary conditions to a group of grid functions
               given by their indices
               This routine is called by the various BndRobinXXX wrappers.

               Although it is currently limited to handle 1D, 2D, or 3D
               variables only it can easily be extended for higher dimensions
               by adapting the appropriate macros.
   @enddesc

   @var        GH
   @vdesc      Pointer to CCTK grid hierarchy
   @vtype      cGH *
   @vio        in
   @endvar
   @var        stencil
   @vdesc      stencil width array
   @vtype      int [ dimension of variable ]
   @vio        in
   @endvar
   @var        finf
   @vdesc      value of f at infimum
   @vtype      CCTK_REAL
   @vio        in
   @endvar
   @var        npow
   @vdesc      power of decay rate
   @vtype      int
   @vio        in
   @endvar
   @var        first_var
   @vdesc      index of first variable to apply boundary conditions to
   @vtype      int
   @vio        in
   @endvar
   @var        num_vars
   @vdesc      number of variables
   @vtype      int
   @vio        in
   @endvar

   @calls      CCTK_VarTypeI
               CCTK_GroupDimFromVarI
               CCTK_NumTimeLevelsFromVarI
               ROBIN_BOUNDARY
   @history
   @hdate      Tue 10 Apr 2001
   @hauthor    Thomas Radke
   @hdesc      Merged separate routines for 1D, 2D, and 3D
               into a single generic routine
   @endhistory

   @returntype int
   @returndesc
                0 for success
               -1 if variable dimension is not supported
               -2 if NULL pointer passed as stencil width array
               -3 if stencil width is other than 1
               -4 if variable type is not supported
               -5 if variable dimension is other than 3D
   @endreturndesc
@@*/
static int ApplyBndRobin (cGH *GH,
                          int stencil[],
                          CCTK_REAL finf,
                          int npow,
                          int first_var,
                          int num_vars)
{
  int i, j;
  int var, vtype, gindex, gdim, timelvl;
  int doBC[2*MAXDIM], lsh[MAXDIM], offset[MAXDIM];
  SymmetryGHex *sGHex;
  const CCTK_REAL *radius;
  char coord_system_name[20];


  /* get the group index of the variables */
  gindex = CCTK_GroupIndexFromVarI (first_var);

  /* get the number of dimensions and the variables' type */
  gdim  = CCTK_GroupDimI (gindex);
  vtype = CCTK_VarTypeI (first_var);

  /* make sure we can deal with this number of dimensions */
  if (gdim > MAXDIM)
  {
    CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
                "ApplyBndRobin: Variable dimension of %d not supported", gdim);
    return (-1);
  }

  /* check the stencil width */
  if (! stencil)
  {
    CCTK_WARN (1, "ApplyBndRobin: NULL pointer passed for stencil width array");
    return (-2);
  }

  for (i = 0; i < gdim; i++)
  {
    if (stencil[i] != 1)
    {
      CCTK_WARN (1, "ApplyBndRobin: Stencil width must be 1 for Robin boundary conditions");
      return (-3);
    }
    offset[i] = i == 0 ? 1 : offset[i-1] * GH->cctk_lsh[i-1];
  }

  /* initialize array for variables with less dimensions than MAXDIM
     so that we can use the INDEX_3D macro later on */
  memset (lsh, 0, sizeof (lsh));

  /* get the current timelevel */
  timelvl = CCTK_NumTimeLevelsFromVarI (first_var) - 1;

  /* Robin boundaries need the underlying radius coordinate */
  sprintf (coord_system_name, "spher%dd", gdim);
  radius = GH->data[CCTK_CoordIndex (-1, "r", coord_system_name)][0];

  /* see if we have a symmetry array */
  sGHex = (SymmetryGHex *) CCTK_GHExtension (GH, "Symmetry");

  /* now loop over all variables */
  for (var = first_var; var < first_var + num_vars; var++)
  {
    /* Apply condition if:
       + boundary is not a symmetry boundary (no symmetry or unset(=unsed))
       + boundary is a physical boundary
       + have enough grid points
    */
    memset (doBC, 1, sizeof (doBC));
    if (sGHex)
    {
      for (i = 0; i < 2 * gdim; i++)
      {
        doBC[i] = sGHex->GFSym[var][i] == GFSYM_NOSYM ||
                  sGHex->GFSym[var][i] == GFSYM_UNSET;
      }
    }
    for (i = 0; i < gdim; i++)
    {
      lsh[i]       = GH->cctk_lsh[i];
      doBC[i*2]   &= GH->cctk_lsh[i] > 1 && GH->cctk_bbox[i*2];
      doBC[i*2+1] &= GH->cctk_lsh[i] > 1 && GH->cctk_bbox[i*2+1];
    }

    switch (CCTK_VarTypeI (var))
    {
      case CCTK_VARIABLE_CHAR:
        ROBIN_BOUNDARY (CCTK_CHAR); break;

      case CCTK_VARIABLE_INT:
        ROBIN_BOUNDARY (CCTK_INT); break;

      case CCTK_VARIABLE_REAL:
        ROBIN_BOUNDARY (CCTK_REAL); break;

#ifdef CCTK_INT2
      case CCTK_VARIABLE_INT2:
        ROBIN_BOUNDARY (CCTK_INT2); break;
#endif

#ifdef CCTK_INT4
      case CCTK_VARIABLE_INT4:
        ROBIN_BOUNDARY (CCTK_INT4); break;
#endif

#ifdef CCTK_INT8
      case CCTK_VARIABLE_INT8:
        ROBIN_BOUNDARY (CCTK_INT8); break;
#endif

#ifdef CCTK_REAL4
      case CCTK_VARIABLE_REAL4:
        ROBIN_BOUNDARY (CCTK_REAL4); break;
#endif

#ifdef CCTK_REAL8
      case CCTK_VARIABLE_REAL8:
        ROBIN_BOUNDARY (CCTK_REAL8); break;
#endif

#ifdef CCTK_REAL16
      case CCTK_VARIABLE_REAL16:
        ROBIN_BOUNDARY (CCTK_REAL16); break;
#endif

      default:
        CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
                    "ApplyBndRobin: Unsupported variable type %d for variable '%s'",
                    CCTK_VarTypeI (var), CCTK_VarName (var));
        return (-4);
    }
  }

  return(0);
}