Adding directional boundary conditions for everything apart from Robin.

See thorn documentation for details

Next to come, boundary conditions for grid functions which aren't CCTK_REALs.


git-svn-id: http://svn.cactuscode.org/arrangements/CactusBase/Boundary/trunk@131 6a38eb6e-646e-4a02-a296-d141613ad6c4

1 files changed, 392 insertions, 113 deletions

diff --git a/src/FlatBoundary.c b/src/FlatBoundary.c
index 2cc8780..616f2da 100644
--- a/src/FlatBoundary.c
+++ b/src/FlatBoundary.c
@@ -8,7 +8,7 @@
    @version $Header$
  @@*/
 
-/*#define DEBUG_BOUND*/
+/*#define DEBUG_BOUNDARY*/
 
 #include <stdio.h>
 #include <assert.h>
@@ -28,30 +28,29 @@
 /* Internal routine prototypes */
 
 static int BndApplyFlat3Di(cGH *GH, 
-			   int  gdim, 
-			   int *doBC, 
+			   int doBC, 
 			   int *lssh,
-			   int *stencil, 
+			   int *stencilarray, 
 			   CCTK_REAL *var);
 
 static int BndApplyFlat2Di(cGH *GH, 
-			   int  gdim, 
-			   int *doBC, 
+			   int doBC, 
 			   int *lssh,
-			   int *stencil, 
+			   int *stencil_array, 
 			   CCTK_REAL *var);
 
 static int BndApplyFlat1Di(cGH *GH, 
-			   int  gdim, 
-			   int *doBC, 
+			   int doBC, 
 			   int *lssh,
-			   int *stencil, 
+			   int *stencil_array, 
 			   CCTK_REAL *var);
 
 static int ApplyBndFlat(cGH *GH, 
-			int *stencil, 
+			int stencil,
+			int *stencil_array, 
 			int first_var,
-			int num_vars);
+			int num_vars,
+			int dir);
 
 /* Local variables */
 
@@ -62,6 +61,39 @@ static int ApplyBndFlat(cGH *GH,
  ********************************************************************/
 
 /*@@
+   @routine    BndFlatDirGI
+   @date       Sun Jan 21 2001
+   @author     Gabrielle Allen
+   @desc 
+   Apply flat boundary conditions by group index in given direction
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+int BndFlatDirGI(cGH *GH, int stencil, int dir, int gi) 
+{
+  int numvars, first_vi;
+
+  first_vi = CCTK_FirstVarIndexI(gi);
+  numvars  = CCTK_NumVarsInGroupI(gi);
+
+  return ApplyBndFlat(GH,stencil,NULL,first_vi,numvars,dir);
+}
+
+void CCTK_FCALL CCTK_FNAME(BndFlatDirGI)
+     (int *ierr, cGH *GH, int *stencil, int *dir, int *gi)
+{
+  *ierr = BndFlatDirGI(GH, *stencil, *dir, *gi);
+  return;
+}
+
+
+/*@@
    @routine    BndFlatGI
    @date       Thu Mar  2 11:11:40 2000
    @author     Gerd Lanfermann
@@ -83,7 +115,7 @@ int BndFlatGI(cGH *GH, int *stencil, int gi)
   first_vi = CCTK_FirstVarIndexI(gi);
   numvars  = CCTK_NumVarsInGroupI(gi);
 
-  return ApplyBndFlat(GH,stencil,first_vi,numvars);
+  return ApplyBndFlat(GH,-1,stencil,first_vi,numvars,0);
 }
 
 void CCTK_FCALL CCTK_FNAME(BndFlatGI)
@@ -94,6 +126,51 @@ void CCTK_FCALL CCTK_FNAME(BndFlatGI)
 }
 
 
+ /* ********************************************************************/
+
+
+/*@@
+   @routine    BndFlatDirGN
+   @date       Sun Jan 21 2001
+   @author     Gabrielle Allen
+   @desc 
+   Apply flat boundary conditions by group name in given direction
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+int BndFlatDirGN(cGH *GH, int stencil, int dir, const char *gn)
+{
+  int retval;
+  int gi;
+  gi   = CCTK_GroupIndex(gn);
+  
+  if (gi>-1)
+  {
+    retval=BndFlatDirGI(GH, stencil, dir, gi);
+  }
+  else 
+  {
+    CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING,
+	       "BndFlatDirGN: Grid variable %s not found",gn);
+    retval=-1;
+  }
+  return retval;
+}
+ 
+void CCTK_FCALL CCTK_FNAME(BndFlatDirGN)
+     (int *ierr, cGH *GH, int *stencil, int *dir, ONE_FORTSTRING_ARG)
+{
+  ONE_FORTSTRING_CREATE(gn)
+  *ierr = BndFlatDirGN(GH, *stencil, *dir, gn);
+  free(gn);
+  return;
+}
 
 
 /*@@
@@ -140,6 +217,37 @@ void CCTK_FCALL CCTK_FNAME(BndFlatGN)
 }
 
 
+ /* ********************************************************************/
+
+
+/*@@
+   @routine    BndFlatDirVI
+   @date       Sun Jan 21 2001
+   @author     Gabrielle Allen
+   @desc 
+   Apply flat boundary conditions by variable index in given direction
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+int BndFlatDirVI(cGH *GH, int stencil, int dir, int vi) 
+{
+  return ApplyBndFlat(GH,stencil,NULL,vi,1,dir);
+}
+
+void CCTK_FCALL CCTK_FNAME(BndFlatDirVI)
+     (int *ierr, cGH *GH, int *stencil, int *dir, int *vi)
+{
+  *ierr = BndFlatDirVI(GH, *stencil, *dir, *vi);
+  return;
+}
+
+
 
 /*@@
    @routine    BndFlatVI
@@ -158,7 +266,7 @@ void CCTK_FCALL CCTK_FNAME(BndFlatGN)
 
 int BndFlatVI(cGH *GH, int *stencil, int vi) 
 {
-  return ApplyBndFlat(GH,stencil,vi,1);
+  return ApplyBndFlat(GH,-1,stencil,vi,1,0);
 }
 
 void CCTK_FCALL CCTK_FNAME(BndFlatVI)
@@ -169,6 +277,54 @@ void CCTK_FCALL CCTK_FNAME(BndFlatVI)
 }
 
 
+ /* ********************************************************************/
+
+
+/*@@
+   @routine    BndFlatDirVN
+   @date       Sun Jan 21 2001
+   @author     Gabrielle Allen
+   @desc 
+   Apply flat boundary conditions by variable name in given direction
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+int BndFlatDirVN(cGH *GH, int stencil, int dir, const char *vn)
+{
+  int  vi;
+  int  retval;
+
+  vi   = CCTK_VarIndex(vn);
+  
+  if (vi>-1)
+  {
+    retval = BndFlatDirVI(GH, stencil, dir, vi);
+  }
+  else
+  {
+    retval = -1;
+    CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING,
+	       "BndFlatDirVN: Grid variable %s not found",vn);
+  }
+  return retval;
+}
+
+void CCTK_FCALL CCTK_FNAME(BndFlatDirVN)
+     (int *ierr, cGH *GH, int *stencil, int *dir, ONE_FORTSTRING_ARG)
+{  
+  ONE_FORTSTRING_CREATE(vn)
+  *ierr = BndFlatDirVN(GH, *stencil, *dir, vn);
+  free(vn);
+  return;
+}
+  
+
 
 /*@@
    @routine    BndFlatVN
@@ -236,118 +392,208 @@ void CCTK_FCALL CCTK_FNAME(BndFlatVN)
 @@*/
 
 int ApplyBndFlat(cGH *GH, 
-		 int *stencil, 
+		 int stencil,
+		 int *stencil_array, 
 		 int first_var,
-		 int num_vars)
+		 int num_vars,
+		 int dir)
 {
   int symmetry_handle;  /* handle for the optional symmetry structure */
   int vi, gi, dim;
   int idim; 
-  int berr=0,ierr=0; 
-  int *doBC, *dstag, *lssh;
+  int berr=0,ierr=0;
+  int retval=0;
+  int dirloop1,dirloop2;
+  int doBC;
+  int count;
+  int *dstag;
+  int *sw;
+  int *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)
+#ifdef DEBUG_BOUNDARY
+  printf("Input arguments for ApplyBndFlat:\n");
+  printf("GH = %x\n",GH);
+  printf("stencil = %d\n",stencil);
+  printf("stencil_array = %x\n",stencil_array);
+  printf("first_var = %d\n",first_var);
+  printf("num_vars = %d\n",num_vars); 
+  printf("dir = %d\n",dir);
+  printf("-----------------------------------\n");
+#endif
+
+   /* Get group dimension */
+  dim = CCTK_GroupDimFromVarI(first_var);
+
+  if (dim<1 || dim>3)
   {
-    sGHex = NULL;
+    CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING,
+	       "ApplyBndScalar: Invalid grid function dimension %d",
+	       dim);
+    retval=-1;
   }
   else
   {
-    sGHex  = (SymmetryGHex*)GH->extensions[symmetry_handle];
-  }
+    /* 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);
+    /* Set up stencil width array if needed */
+    if (stencil_array && dir == 0)
+    {
+      sw = stencil_array;
+    }
+    else if (dir==0)
+    {
+      CCTK_WARN(0,"ApplyBndScalar: Direction 0 invalid for directional BC");
+    }
+    else if (dir>dim||-dir>dim)
+    {
+      CCTK_VWarn(0,__LINE__,__FILE__,CCTK_THORNSTRING,
+		 "ApplyBndScalar: Direction %d greater than dimension %d",
+		 dir,dim);
+    }
+    else
+    {
+      int i;
+      /* Should be single direction */
+      sw  = (int *)malloc(dim*sizeof(int));
+      if (sw)
+      {
+	for (i=0;i<dim;i++)
+	{
+	  sw[i] = 0;
+	}
+	if (dir > 0)
+	{
+	  sw[dir-1] = stencil;
+	}
+	else
+	{
+	  sw[-dir-1] = stencil;
+	}
+      }
+      else
+      {
+	CCTK_WARN(0,"ApplyBndScalar: Malloc sw failed");
+      }
+    }
 
-  /* 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);
-  if (ierr<0) CCTK_VWarn(1,__LINE__,__FILE__,"Boundary",
-	"Cannot get stagger information for group '%s' (GF '$s')\n",
-	 CCTK_GroupName(gi), CCTK_VarName(first_var));
-  
+    /* allocate arrays */
+    dstag  = (int *)malloc(dim*sizeof(int));
+    lssh   = (int *)malloc(dim*sizeof(int));
 
-  /* get the current timelevel */
-  timelevel = CCTK_NumTimeLevelsFromVarI (first_var) - 1;
-  /*  if (timelevel > 0)
-      timelevel--;*/
+    /* get the directional staggering of the group */
+    gi   = CCTK_GroupIndexFromVarI(first_var);
+    ierr = CCTK_GroupStaggerDirArrayGI(dstag, dim, gi);
 
-  for (vi=first_var; vi<first_var+num_vars; vi++) 
-  {
+    /* get the current timelevel */
+    timelevel = CCTK_NumTimeLevelsFromVarI (first_var) - 1;
 
-    /* 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)
+    /* Treat all boundaries or just one? */
+    if (dir>0)
     {
-      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)];
-      }
+      dirloop1 = 2*(dir-1)+1;
+      dirloop2 = dirloop1+1;
+    }
+    else if (dir<0)
+    {
+      dirloop1 = 2*(-dir-1);
+      dirloop2 = dirloop1+1;
     }
     else
     {
-      for (idim=0;idim<dim;idim++) 
+      dirloop1 = 0;
+      dirloop2 = 2*dim;
+    }
+  
+    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          */
+
+      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)];
+	lssh[idim] = (int)GH->cctk_lssh[CCTK_LSSH_IDX(dstag[idim],idim)];
       }
+
+      for (count=dirloop1;count<dirloop2;count++)
+      {
+	if (sGHex)
+	{
+	  doBC = (
+		  ( (sGHex->GFSym[vi][count]==GFSYM_NOSYM) ||
+		    (sGHex->GFSym[vi][count]==GFSYM_UNSET) ) &&
+		  GH->cctk_lsh[count/2]>1 && GH->cctk_bbox[count]) ?
+	    count : -1;
+	}
+	else
+	{
+	  doBC  = (GH->cctk_lsh[count/2]>1 && GH->cctk_bbox[count]) ?
+	    count : -1 ;
+	}
+	
+	if (doBC >=0)
+	{
+	  switch (dim)
+	  {
+	  case 1: 	
+	    berr = BndApplyFlat1Di(GH, 
+				   doBC, 
+				   lssh,
+				   sw,
+				   GH->data[vi][timelevel]); 
+	    break;
+	  case 2:
+	    berr = BndApplyFlat2Di(GH, 
+				   doBC, 
+				   lssh,
+				   sw,
+				   GH->data[vi][timelevel]); 
+	    break;
+	  case 3:
+	    berr = BndApplyFlat3Di(GH, 
+				   doBC, 
+				   lssh,
+				   sw,
+				   GH->data[vi][timelevel]); 
+	    break;
+	  default : 
+	    berr = -1; 
+	    CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING, 
+		       "ApplyBndFlat: No BC for dim>3: grid variable %s",
+		       CCTK_VarName(vi));
+	  }
+	}
+      }
+      /*	  retval = (berr>-1) ? 0 : -1;*/
     }
+  
+    free(dstag);
+    free(lssh);
 
-    switch (dim) 
+    if (!stencil_array)
     {
-      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_VWarn(1,__LINE__,__FILE__,"Boundary", 
-		   "ApplyBndFlat: No BC for dim>3: grid variable '%s'",
-		   CCTK_VarName(vi));
+      free(sw);
     }
-    berr = (berr>-1) ? 0 : -1;
   }
-  
-  free(dstag);
-  free(doBC);
-  free(lssh);
-  
-  return(berr);
+
+  return retval;
 }
 
 
+
 /*@@
    @routine    BndApplyFlat3Di
    @date       Thu Mar  2 11:14:08 2000
@@ -364,8 +610,7 @@ int ApplyBndFlat(cGH *GH,
 @@*/
 
 int BndApplyFlat3Di(cGH *GH, 
-		    int  gdim, 
-		    int *doBC, 
+		    int doBC, 
 		    int *lssh,
 		    int *stencil, 
 		    CCTK_REAL *var) 
@@ -373,8 +618,11 @@ int BndApplyFlat3Di(cGH *GH,
   int i,j,k,sw;
 
   /* lower x */
-  if (doBC[0] == 1)
+  if (doBC == 0)
   {
+#ifdef DEBUG_BOUNDARY
+    printf("Applying flat boundary (lower x)\n");
+#endif
     for (k=0;k<GH->cctk_lsh[2];k++)
     {
       for (j=0;j<GH->cctk_lsh[1];j++)
@@ -389,8 +637,11 @@ int BndApplyFlat3Di(cGH *GH,
   }
   
   /* upper x */
-  if (doBC[1] == 1)
+  if (doBC == 1)
   {
+#ifdef DEBUG_BOUNDARY
+    printf("Applying flat boundary (upper x)\n");
+#endif
     for (k=0;k<GH->cctk_lsh[2];k++)
     {
       for (j=0;j<GH->cctk_lsh[1];j++)
@@ -405,8 +656,11 @@ int BndApplyFlat3Di(cGH *GH,
   }
   
   /*lower y */
-  if (doBC[2] == 1)
+  if (doBC == 2)
   {
+#ifdef DEBUG_BOUNDARY
+    printf("Applying flat boundary (lower y)\n");
+#endif
     for (k=0;k<GH->cctk_lsh[2];k++)
     {
       for (i=0;i<GH->cctk_lsh[0];i++)
@@ -421,8 +675,11 @@ int BndApplyFlat3Di(cGH *GH,
   }
   
   /* upper y */
-  if (doBC[3] == 1)
+  if (doBC == 3)
   {
+#ifdef DEBUG_BOUNDARY
+    printf("Applying flat boundary (upper y)\n");
+#endif
     for (k=0;k<GH->cctk_lsh[2];k++)
     {
       for (i=0;i<GH->cctk_lsh[0];i++)
@@ -437,8 +694,11 @@ int BndApplyFlat3Di(cGH *GH,
   }
 
   /* lower z */
-  if (doBC[4] == 1)
+  if (doBC == 4)
   {
+#ifdef DEBUG_BOUNDARY
+    printf("Applying flat boundary (lower z)\n");
+#endif
     for (j=0;j<GH->cctk_lsh[1];j++)
     {
       for (i=0;i<GH->cctk_lsh[0];i++)
@@ -453,8 +713,11 @@ int BndApplyFlat3Di(cGH *GH,
   }
 
   /* upper z */
-  if (doBC[5] == 1)
+  if (doBC == 5)
   {
+#ifdef DEBUG_BOUNDARY
+    printf("Applying flat boundary (upper z)\n");
+#endif
     for (j=0;j<GH->cctk_lsh[1];j++)
     {
       for (i=0;i<GH->cctk_lsh[0];i++)
@@ -487,8 +750,7 @@ int BndApplyFlat3Di(cGH *GH,
 
 @@*/
 int BndApplyFlat2Di(cGH *GH, 
-		    int  gdim, 
-		    int *doBC, 
+		    int doBC, 
 		    int *lssh,
 		    int *stencil, 
 		    CCTK_REAL *var) 
@@ -496,8 +758,11 @@ int BndApplyFlat2Di(cGH *GH,
   int i,j,sw;
     
   /* lower x */
-  if (doBC[0] == 1)
+  if (doBC == 0)
   {
+#ifdef DEBUG_BOUNDARY
+    printf("Applying flat boundary (lower x)\n");
+#endif
     for (j=0;j<GH->cctk_lsh[1];j++)
     {
       for (sw=0;sw<stencil[0];sw++)
@@ -509,8 +774,11 @@ int BndApplyFlat2Di(cGH *GH,
   }
   
   /* upper x */
-  if (doBC[1] == 1)
+  if (doBC == 1)
   {
+#ifdef DEBUG_BOUNDARY
+    printf("Applying flat boundary (upper x)\n");
+#endif
     for (j=0;j<GH->cctk_lsh[1];j++)
     {
       for (sw=0;sw<stencil[0];sw++)
@@ -522,8 +790,11 @@ int BndApplyFlat2Di(cGH *GH,
   }
 
   /* lower y */
-  if (doBC[2] == 1)
+  if (doBC == 2)
   {
+#ifdef DEBUG_BOUNDARY
+    printf("Applying flat boundary (lower y)\n");
+#endif
     for (i=0;i<GH->cctk_lsh[0];i++)
     {
       for (sw=0;sw<stencil[1];sw++)
@@ -535,8 +806,11 @@ int BndApplyFlat2Di(cGH *GH,
   }
   
   /* upper y */
-  if (doBC[3] == 1)
+  if (doBC == 3)
   {
+#ifdef DEBUG_BOUNDARY
+    printf("Applying flat boundary (upper y)\n");
+#endif
     for (i=0;i<GH->cctk_lsh[0];i++)
     {
       for (sw=0;sw<stencil[1];sw++)
@@ -568,16 +842,18 @@ int BndApplyFlat2Di(cGH *GH,
 @@*/
 
 int BndApplyFlat1Di(cGH *GH, 
-		    int  gdim, 
-		    int *doBC, 
+		    int doBC, 
 		    int *lssh,
 		    int *stencil, 
 		    CCTK_REAL *var) 
 { int sw;
 
   /* lower x */
-  if (doBC[0] == 1)
+  if (doBC == 0)
   {
+#ifdef DEBUG_BOUNDARY
+    printf("Applying flat boundary (lower x)\n");
+#endif
     for (sw=0;sw<stencil[0];sw++)
     {
       var[CCTK_GFINDEX1D(GH,sw)]
@@ -586,8 +862,11 @@ int BndApplyFlat1Di(cGH *GH,
   }
   
   /* upper x */
-  if (doBC[1] == 1)
+  if (doBC == 1)
   {
+#ifdef DEBUG_BOUNDARY
+    printf("Applying flat boundary (upper x)\n");
+#endif
     for (sw=0;sw<stencil[0];sw++)
     {
       var[CCTK_GFINDEX1D(GH,GH->cctk_lsh[0]-sw-1)]