path: root/src/FlatBoundary.c

 /*@@
   @file      FlatBoundaryWrappers.c
   @date      
   @author    Gerd Lanfermann
   @desc 
   Flat boundary conditions in 1,2 and 3D
   @enddesc 
   @version $Header$
 @@*/

/*#define DEBUG_BOUND*/

#include <stdio.h>
#include <assert.h>
#include <stdlib.h>
#include <ctype.h>
#include <stdarg.h>
#include <string.h>

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

#include "BoundarySymmetries.h"
#include "Symmetry.h"


/* Internal routine prototypes */

static int BndApplyFlat3Di(cGH *GH, 
			   int  gdim, 
			   int *doBC, 
			   int *lssh,
			   int *stencil, 
			   CCTK_REAL *var);

static int BndApplyFlat2Di(cGH *GH, 
			   int  gdim, 
			   int *doBC, 
			   int *lssh,
			   int *stencil, 
			   CCTK_REAL *var);

static int BndApplyFlat1Di(cGH *GH, 
			   int  gdim, 
			   int *doBC, 
			   int *lssh,
			   int *stencil, 
			   CCTK_REAL *var);

static int ApplyBndFlat(cGH *GH, 
			int *stencil, 
			int first_var,
			int num_vars);

/* Local variables */



/********************************************************************
 ********************    External Routines   ************************
 ********************************************************************/

/*@@
   @routine    BndFlatGI
   @date       Thu Mar  2 11:11:40 2000
   @author     Gerd Lanfermann
   @desc 
   Apply flat boundary conditions by group index
   @enddesc 
   @calls     
   @calledby   
   @history 
 
   @endhistory 

@@*/

int BndFlatGI(cGH *GH, int *stencil, int gi) 
{
  int numvars, first_vi;

  first_vi = CCTK_FirstVarIndexI(gi);
  numvars  = CCTK_NumVarsInGroupI(gi);

  return ApplyBndFlat(GH,stencil,first_vi,numvars);
}

void CCTK_FCALL CCTK_FNAME(BndFlatGI)
     (int *ierr, cGH *GH, int *sw, int *gi)
{
  *ierr = BndFlatGI(GH, sw, *gi);
  return;
}




/*@@
   @routine    BndFlatGN
   @date       Thu Mar  2 11:11:40 2000
   @author     Gerd Lanfermann
   @desc 
   Apply flat boundary conditions by group  name
   @enddesc 
   @calls     
   @calledby   
   @history 
 
   @endhistory 

@@*/

int BndFlatGN(cGH *GH, int *sw, const char *gn)
{
  int retval;
  int gi;
  gi   = CCTK_GroupIndex(gn);
  
  if (gi>-1)
  {
    retval=BndFlatGI(GH, sw, gi);
  }
  else 
  {
    CCTK_VWarn(1,__LINE__,__FILE__,"Boundary",
	       "BndFlatGN: Grid variable %s not found",gn);
    retval=-1;
  }
  return retval;
}
 
void CCTK_FCALL CCTK_FNAME(BndFlatGN)
     (int *ierr, cGH *GH, int *sw, ONE_FORTSTRING_ARG)
{
  ONE_FORTSTRING_CREATE(gn)
  *ierr = BndFlatGN(GH, sw, gn);
  free(gn);
  return;
}



/*@@
   @routine    BndFlatVI
   @date       Thu Mar  2 11:11:40 2000
   @author     Gerd Lanfermann
   @desc 
   Apply flat boundary conditions by variable index
   @enddesc 
   @calls     
   @calledby   
   @history 
 
   @endhistory 

@@*/

int BndFlatVI(cGH *GH, int *stencil, int vi) 
{
  return ApplyBndFlat(GH,stencil,vi,1);
}

void CCTK_FCALL CCTK_FNAME(BndFlatVI)
     (int *ierr, cGH *GH, int *sw, int *vi)
{
  *ierr = BndFlatVI(GH, sw, *vi);
  return;
}



/*@@
   @routine    BndFlatVN
   @date       Thu Mar  2 11:11:40 2000
   @author     Gerd Lanfermann
   @desc 
   Apply flat boundary conditions by variable name
   @enddesc 
   @calls     
   @calledby   
   @history 
 
   @endhistory 

@@*/

int BndFlatVN(cGH *GH, 
	      int *sw,
	      const char *vn)
{
  int  vi;
  int  retval;

  vi   = CCTK_VarIndex(vn);
  
  if (vi>-1)
  {
    retval = BndFlatVI(GH, sw, vi);
  }
  else
  {
    retval = -1;
    CCTK_VWarn(1,__LINE__,__FILE__,"Boundary",
	       "BndFlatVN: Grid variable %s not found",vn);
  }
  return retval;
}

void CCTK_FCALL CCTK_FNAME(BndFlatVN)
     (int *ierr, cGH *GH, int *sw, ONE_FORTSTRING_ARG)
{  
  ONE_FORTSTRING_CREATE(vn)
  *ierr = BndFlatVN(GH, sw, vn);
  free(vn);
  return;
}
  


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

/*@@
   @routine    ApplyBndFlat
   @date       Jul 5 2000
   @author     Gabrielle Allen
   @desc 
   Wrapper for setting up and calling flat boundary conditions
   @enddesc 
   @history 
 
   @endhistory 

@@*/

int ApplyBndFlat(cGH *GH, 
		 int *stencil, 
		 int first_var,
		 int num_vars)
{
  int symmetry_handle;  /* handle for the optional symmetry structure */
  int vi, gi, dim;
  int idim;
  int berr,ierr;
  int *doBC, *dstag, *lssh;
  int timelevel;        /* timelevel that condition applied on */
  SymmetryGHex *sGHex;

  /* See if we have a symmetry array */
  symmetry_handle = CCTK_GHExtensionHandle("Symmetry");
  if (symmetry_handle < 0)
  {
    sGHex = NULL;
  }
  else
  {
    sGHex  = (SymmetryGHex*)GH->extensions[symmetry_handle];
  }

  /* Get group dimension */
  dim = CCTK_GroupDimFromVarI(first_var);

  /* allocate arrays */
  doBC   = (int *)malloc((2*dim)*sizeof(int)); 
  dstag  = (int *)malloc(dim*sizeof(int));
  lssh   = (int *)malloc(dim*sizeof(int));

  /* get the directional staggering of the group */
  gi   = CCTK_GroupIndexFromVarI(first_var);
  ierr = CCTK_GroupStaggerDirArrayGI(dstag, dim, gi);

  /* get the current timelevel */
  timelevel = CCTK_NumTimeLevelsFromVarI (first_var) - 1;
  if (timelevel > 0)
    timelevel--;

  for (vi=first_var; vi<first_var+num_vars; vi++) 
  {

    /* Apply condition if:
       + boundary is not a symmetry boundary 
       (no symmetry or unset(=unsed))
       + boundary is a physical boundary
       + have enough grid points
    */
    if (sGHex)
    {
      for (idim=0;idim<dim;idim++) 
      {
	doBC[idim*2]   = (((sGHex->GFSym[vi][idim*2]==GFSYM_NOSYM)||
			   (sGHex->GFSym[vi][idim*2]==GFSYM_UNSET)) &&
			  GH->cctk_lsh[idim]>1 && GH->cctk_bbox[idim*2]);
	doBC[idim*2+1] = (((sGHex->GFSym[vi][idim*2+1]==GFSYM_NOSYM)||
			   (sGHex->GFSym[vi][idim*2+1]==GFSYM_UNSET)) &&
			  GH->cctk_lsh[idim]>1 && GH->cctk_bbox[idim*2+1]);
	lssh[idim] = GH->cctk_lssh[CCTK_LSSH_IDX(dstag[idim],idim)];
      }
    }
    else
    {
      for (idim=0;idim<dim;idim++) 
      {
	doBC[idim*2]   = (GH->cctk_lsh[idim]>1 && GH->cctk_bbox[idim*2]);
	doBC[idim*2+1] = (GH->cctk_lsh[idim]>1 && GH->cctk_bbox[idim*2+1]);
	lssh[idim]     = GH->cctk_lssh[CCTK_LSSH_IDX(dstag[idim],idim)];
      }
    }

    switch (dim) 
    {
      case 1: berr = BndApplyFlat1Di(GH, 
				     dim,
				     doBC, 
				     lssh,
				     stencil,
				     GH->data[vi][timelevel]); break;
      case 2: berr = BndApplyFlat2Di(GH, 
				     dim,
				     doBC, 
				     lssh,
				     stencil,
				     GH->data[vi][timelevel]); break;
      case 3: berr = BndApplyFlat3Di(GH, 
				     dim,
				     doBC,
				     lssh,
				     stencil,
				     GH->data[vi][timelevel]); break;     
      default : berr = -1; CCTK_WARN(1, "ApplyBndFlat: No BC for dim>3");
    }
    berr = (berr>-1) ? 0 : -1;
  }
  
  free(dstag);
  free(doBC);
  free(lssh);
  
  return(ierr);
}


/*@@
   @routine    BndApplyFlat3Di
   @date       Thu Mar  2 11:14:08 2000
   @author     Gerd Lanfermann
   @desc 
     Apply flat boundary conditions to 3d variables
   @enddesc 
   @calls     
   @calledby   
   @history 
 
   @endhistory 

@@*/

int BndApplyFlat3Di(cGH *GH, 
		    int  gdim, 
		    int *doBC, 
		    int *lssh,
		    int *stencil, 
		    CCTK_REAL *var) 
{
  int i,j,k,sw;

  /* lower x */
  if (doBC[0] == 1)
  {
    for (k=0;k<GH->cctk_lsh[2];k++)
    {
      for (j=0;j<GH->cctk_lsh[1];j++)
      {
	for (sw=0;sw<stencil[0];sw++)
	{
	  var[CCTK_GFINDEX3D(GH,sw,j,k)]
	    = var[CCTK_GFINDEX3D(GH,stencil[0],j,k)];
	}
      }
    }
  }
  
  /* upper x */
  if (doBC[1] == 1)
  {
    for (k=0;k<GH->cctk_lsh[2];k++)
    {
      for (j=0;j<GH->cctk_lsh[1];j++)
      {
	for (sw=0;sw<stencil[0];sw++)
	{
	  var[CCTK_GFINDEX3D(GH,GH->cctk_lsh[0]-sw-1,j,k)]
	    = var[CCTK_GFINDEX3D(GH,GH->cctk_lsh[0]-stencil[0]-1,j,k)];
	}
      }
    }
  }
  
  /*lower y */
  if (doBC[2] == 1)
  {
    for (k=0;k<GH->cctk_lsh[2];k++)
    {
      for (i=0;i<GH->cctk_lsh[0];i++)
      {
	for (sw=0;sw<stencil[1];sw++)
	{
	  var[CCTK_GFINDEX3D(GH,i,sw,k)]
	    = var[CCTK_GFINDEX3D(GH,i,stencil[1],k)];
	}
      }
    }
  }
  
  /* upper y */
  if (doBC[3] == 1)
  {
    for (k=0;k<GH->cctk_lsh[2];k++)
    {
      for (i=0;i<GH->cctk_lsh[0];i++)
      {
	for (sw=0;sw<stencil[1];sw++)
	{
	  var[CCTK_GFINDEX3D(GH,i,GH->cctk_lsh[1]-sw-1,k)]
	    = var[CCTK_GFINDEX3D(GH,i,GH->cctk_lsh[1]-stencil[1]-1,k)];
	}
      }
    }
  }

  /* lower z */
  if (doBC[4] == 1)
  {
    for (j=0;j<GH->cctk_lsh[1];j++)
    {
      for (i=0;i<GH->cctk_lsh[0];i++)
      {
	for (sw=0;sw<stencil[2];sw++)
	{
	  var[CCTK_GFINDEX3D(GH,i,j,sw)]
	    = var[CCTK_GFINDEX3D(GH,i,j,stencil[2])];
	}
      }
    }
  }

  /* upper z */
  if (doBC[5] == 1)
  {
    for (j=0;j<GH->cctk_lsh[1];j++)
    {
      for (i=0;i<GH->cctk_lsh[0];i++)
      {
	for (sw=0;sw<stencil[2];sw++)
	{
	  var[CCTK_GFINDEX3D(GH,i,j,GH->cctk_lsh[2]-sw-1)]
	    = var[CCTK_GFINDEX3D(GH,i,j,GH->cctk_lsh[2]-stencil[2]-1)];
	}
      }
    }
  }
  return(0);
}



/*@@
   @routine    BndApplyFlat2Di
   @date       Thu Mar  2 11:14:08 2000
   @author     Gerd Lanfermann
   @desc 
     Apply flat boudnary conditions to 2d variables
   @enddesc 
   @calls     
   @calledby   
   @history 
 
   @endhistory 

@@*/
int BndApplyFlat2Di(cGH *GH, 
		    int  gdim, 
		    int *doBC, 
		    int *lssh,
		    int *stencil, 
		    CCTK_REAL *var) 
{
  int i,j,sw;
    
  /* lower x */
  if (doBC[0] == 1)
  {
    for (j=0;j<GH->cctk_lsh[1];j++)
    {
      for (sw=0;sw<stencil[0];sw++)
      {
	var[CCTK_GFINDEX2D(GH,sw,j)]
	  = var[CCTK_GFINDEX2D(GH,stencil[0],j)];
      }
    }
  }
  
  /* upper x */
  if (doBC[1] == 1)
  {
    for (j=0;j<GH->cctk_lsh[1];j++)
    {
      for (sw=0;sw<stencil[0];sw++)
      {
	var[CCTK_GFINDEX2D(GH,GH->cctk_lsh[0]-sw-1,j)]
	  = var[CCTK_GFINDEX2D(GH,GH->cctk_lsh[0]-stencil[0]-1,j)];
      }
    }
  }

  /* lower y */
  if (doBC[2] == 1)
  {
    for (i=0;i<GH->cctk_lsh[0];i++)
    {
      for (sw=0;sw<stencil[1];sw++)
      {
	var[CCTK_GFINDEX2D(GH,i,sw)]
	  = var[CCTK_GFINDEX2D(GH,i,stencil[1])];
      }
    }
  }
  
  /* upper y */
  if (doBC[3] == 1)
  {
    for (i=0;i<GH->cctk_lsh[0];i++)
    {
      for (sw=0;sw<stencil[1];sw++)
      {
	var[CCTK_GFINDEX2D(GH,i,GH->cctk_lsh[1]-sw-1)]
	  = var[CCTK_GFINDEX2D(GH,i,GH->cctk_lsh[1]-stencil[1]-1)];
      }
    }
  }

  return(0);
}



/*@@
   @routine    BndApplyFlat1Di
   @date       Thu Mar  2 11:14:08 2000
   @author     Gerd Lanfermann
   @desc 
     Apply flat boudnary conditions to 1d variables
   @enddesc 
   @calls     
   @calledby   
   @history 
 
   @endhistory 

@@*/

int BndApplyFlat1Di(cGH *GH, 
		    int  gdim, 
		    int *doBC, 
		    int *lssh,
		    int *stencil, 
		    CCTK_REAL *var) 
{ int sw;

  /* lower x */
  if (doBC[0] == 1)
  {
    for (sw=0;sw<stencil[0];sw++)
    {
      var[CCTK_GFINDEX1D(GH,sw)]
	= var[CCTK_GFINDEX1D(GH,stencil[0])];
    }
  }
  
  /* upper x */
  if (doBC[1] == 1)
  {
    for (sw=0;sw<stencil[0];sw++)
    {
      var[CCTK_GFINDEX1D(GH,GH->cctk_lsh[0]-sw-1)]
	= var[CCTK_GFINDEX1D(GH,GH->cctk_lsh[0]-stencil[0]-1)];
    }
  }
  return(0);
}



/********************************************************************
 *********************    DEPRECATED: Beta 8   **********************
 ********************************************************************/

int FlatBCGroup(cGH *GH, 
		int *stencil_size, 
		const char *impgrpname)
{
  CCTK_WARN(1,"FlatBCGroup: ROUTINE DEPRECATED, see Beta8 release notices");
  return(BndFlatGN(GH,stencil_size,impgrpname));
}

void CCTK_FCALL CCTK_FNAME(FlatBCGroup)(int *retval, 
					 cGH *GH, 
					 int *stencil_size, 
					 ONE_FORTSTRING_ARG) 
{
  ONE_FORTSTRING_CREATE(impgrpname)
  CCTK_WARN(1,"FlatBCGroup: ROUTINE DEPRECATED, see Beta8 release notices");
  *retval = BndFlatGN(GH,stencil_size,impgrpname);
  free(impgrpname);
}

int FlatBCVar(cGH *GH, 
	      int *stencil_size, 
	      const char *impvarname) 
{
  CCTK_WARN(1,"FlatBCVar: ROUTINE DEPRECATED, see Beta8 release notices");
  return(BndFlatVN(GH,stencil_size,impvarname));
}

void CCTK_FCALL CCTK_FNAME(FlatBCVar)(int *retval, 
				       cGH *GH, 
				       int *stencil_size, 
				       ONE_FORTSTRING_ARG) 
{
  ONE_FORTSTRING_CREATE(impvarname)
  CCTK_WARN(1,"FlatBCVar: ROUTINE DEPRECATED, see Beta8 release notices");
  *retval = BndFlatVN(GH, stencil_size, impvarname);
  free(impvarname);
}

int FlatBCVarI(cGH *GH, int *stencil_size, int vi) 
{
  CCTK_WARN(1,"FlatBCVarI: ROUTINE DEPRECATED, see Beta8 release notices");
  return(BndFlatVI(GH, stencil_size, vi));
}

void CCTK_FCALL CCTK_FNAME(FlatBCVarI)(int *retval, 
					cGH *GH, 
					int *stencil_size, 
					int *vi) 
{
  CCTK_WARN(1,"FlatBCVarI: ROUTINE DEPRECATED, see Beta8 release notices");
  *retval = BndFlatVI(GH, stencil_size, *vi);
}

int FlatBCGroupI(cGH *GH, 
		 int *stencil_size, 
		 int gi) 
{
  CCTK_WARN(1,"FlatBCGroupI: ROUTINE DEPRECATED, see Beta8 release notices");
  return(BndFlatGI(GH, stencil_size, gi));
}

void CCTK_FCALL CCTK_FNAME(FlatBCGroupI)(int *retval, 
					  cGH *GH, 
					  int *stencil_size, 
					  int *gi) 
{
  CCTK_WARN(1,"FlatBCGroupI: ROUTINE DEPRECATED, see Beta8 release notices");
  *retval=BndFlatGI(GH, stencil_size, *gi);
}