Patch from Erik Schnetter:

the following patch (which is only as large because a part of CartGrid3D is
now within an if block, changing indentation) should make CartGrid3D work
with mesh refinement.  I ensures that CartGrid3D calculates grid spacing,
registers coordinates etc. only once, namely for the coarsest level.  On
finer levels, it reuses the information from the coarser levels.  This patch
should not make any differences to unigrid codes.  It only changes behaviour
when cctk_levfac[] is different from one.


git-svn-id: http://svn.cactuscode.org/arrangements/CactusBase/CartGrid3D/trunk@151 c78560ca-4b45-4335-b268-5f3340f3cb52

1 files changed, 307 insertions, 255 deletions

diff --git a/src/CartGrid3D.c b/src/CartGrid3D.c
index 384e41a..5a5a2d0 100644
--- a/src/CartGrid3D.c
+++ b/src/CartGrid3D.c
@@ -41,272 +41,292 @@ void CartGrid3D(CCTK_ARGUMENTS)
   CCTK_REAL lowerx,upperx,lowery,uppery,lowerz,upperz;
   int lowerxi,upperxi,loweryi,upperyi,lowerzi,upperzi;
   char infoline[120];
+  int ierr;
 
   int domainsym[6];
   int cntstag[3];
 
-
-  /* Some compilers warn about variables which might be uninitialized
-     because they are assigned in some if/else statements only */
-  x_origin = y_origin = z_origin = 0.0;
-  lowerxi = loweryi = lowerzi = 0;
-  upperxi = cctk_gsh[0]-1;
-  upperyi = cctk_gsh[1]-1;
-  upperzi = cctk_gsh[2]-1;
-  this_dx = this_dy = this_dz = 0.0;
-
-  /* Avoid origin?  Default is yes */
-  cntstag[0] = no_origin && no_originx && avoid_origin && avoid_originx;
-  cntstag[1] = no_origin && no_originy && avoid_origin && avoid_originy;
-  cntstag[2] = no_origin && no_originz && avoid_origin && avoid_originz;
-
-  /* Determine symmetries of domain */
-  DecodeSymParameters3D(domainsym);
-  
-  /* Calculate physical indices, using symmetries and periodicity */
-  if (domainsym[0]) 
-  {
-    lowerxi = cctk_nghostzones[0];
-  }  
-  else
+  if (cctk_levfac[0] == 1 && cctk_levfac[1] == 1 && cctk_levfac[2] == 1) 
   {
-    lowerxi = 0;
-  }
-  if (domainsym[2]) 
-  {
-    loweryi = cctk_nghostzones[1];
-  }  
-  else
-  {
-    loweryi = 0;
-  }
-  if (domainsym[4]) 
-  {
-    lowerzi = cctk_nghostzones[2];
-  }  
-  else
-  {
-    lowerzi = 0;
-  }
 
-  if (periodic)
-  {
-    if (periodic_x)
+    /* Calculate the coordinate ranges only for the coarsest level */
+    
+    /* Some compilers warn about variables which might be uninitialized
+       because they are assigned in some if/else statements only */
+    x_origin = y_origin = z_origin = 0.0;
+    lowerxi = loweryi = lowerzi = 0;
+    upperxi = cctk_gsh[0]-1;
+    upperyi = cctk_gsh[1]-1;
+    upperzi = cctk_gsh[2]-1;
+    this_dx = this_dy = this_dz = 0.0;
+
+    /* Avoid origin?  Default is yes */
+    cntstag[0] = no_origin && no_originx && avoid_origin && avoid_originx;
+    cntstag[1] = no_origin && no_originy && avoid_origin && avoid_originy;
+    cntstag[2] = no_origin && no_originz && avoid_origin && avoid_originz;
+
+    /* Determine symmetries of domain */
+    DecodeSymParameters3D(domainsym);
+  
+    /* Calculate physical indices, using symmetries and periodicity */
+    if (domainsym[0]) 
     {
-      lowerxi = cctk_nghostzones[0]-1;
-      upperxi = cctk_gsh[0]-1-cctk_nghostzones[0];
-    }
-    if (periodic_y)
+      lowerxi = cctk_nghostzones[0];
+    }  
+    else
     {
-      loweryi = cctk_nghostzones[1]-1;
-      upperyi = cctk_gsh[1]-1-cctk_nghostzones[1];
+      lowerxi = 0;
     }
-    if (periodic_z)
+    if (domainsym[2]) 
     {
-      lowerzi = cctk_nghostzones[2]-1;
-      upperzi = cctk_gsh[2]-1-cctk_nghostzones[2];
-    }
-  }
-
-  dconv = pow(2, cctk_convlevel);
-  iconv = (int)dconv;
-
-  /****************************************************************
-   *
-   *  BYRANGE
-   *
-   *  User gives: minimum and maximum values of coordinates and
-   *              the number of gridpoints on the coarse grid
-   *
-   ***************************************************************/
-  /**************************************************************
-   *
-   *   BOX (-0.5 to 0.5)
-   *
-   *   User gives: number of gridpoints on the coarse grid
-   *
-   **************************************************************/
-
-  if (CCTK_Equals(type,"byrange") || CCTK_Equals(type,"box")) 
-  { 
-    
-    if (CCTK_Equals(type,"box"))
+      loweryi = cctk_nghostzones[1];
+    }  
+    else
     {
-      
-      /*  Coordinates are all -0.5 to 0.5 */
-      xmax1 = 0.5;
-      ymax1 = 0.5;
-      zmax1 = 0.5;
-      
-      xmin1 = -0.5;
-      ymin1 = -0.5;
-      zmin1 = -0.5;
-      
+      loweryi = 0;
     }
+    if (domainsym[4]) 
+    {
+      lowerzi = cctk_nghostzones[2];
+    }  
     else
     {
+      lowerzi = 0;
+    }
 
-      if (xyzmax != -424242)
+    if (periodic)
+    {
+      if (periodic_x)
       {
-          xmax1 = xyzmax;
-          ymax1 = xyzmax;
-          zmax1 = xyzmax;
+	lowerxi = cctk_nghostzones[0];
+	upperxi = cctk_gsh[0]-1-cctk_nghostzones[0];
       }
-      else
+      if (periodic_y)
       {
-          xmax1 = xmax;
-          ymax1 = ymax;
-          zmax1 = zmax;
+	loweryi = cctk_nghostzones[1];
+	upperyi = cctk_gsh[1]-1-cctk_nghostzones[1];
       }
-  
-      if (xyzmin != -424242)
+      if (periodic_z)
       {
-          xmin1 = xyzmin;
-          ymin1 = xyzmin;
-          zmin1 = xyzmin;
+	lowerzi = cctk_nghostzones[2];
+	upperzi = cctk_gsh[2]-1-cctk_nghostzones[2];
       }
-      else
+    }
+
+    dconv = pow(2, cctk_convlevel);
+    iconv = (int)dconv;
+
+    /****************************************************************
+     *
+     *  BYRANGE
+     *
+     *  User gives: minimum and maximum values of coordinates and
+     *              the number of gridpoints on the coarse grid
+     *
+     ***************************************************************/
+    /**************************************************************
+     *
+     *   BOX (-0.5 to 0.5)
+     *
+     *   User gives: number of gridpoints on the coarse grid
+     *
+     **************************************************************/
+
+    if (CCTK_Equals(type,"byrange") || CCTK_Equals(type,"box")) 
+    { 
+      
+      if (CCTK_Equals(type,"box"))
       {
-          xmin1 = xmin;
-          ymin1 = ymin;
-          zmin1 = zmin;
+      
+	/*  Coordinates are all -0.5 to 0.5 */
+	xmax1 = 0.5;
+	ymax1 = 0.5;
+	zmax1 = 0.5;
+	
+	xmin1 = -0.5;
+	ymin1 = -0.5;
+	zmin1 = -0.5;
+	
       }
+      else
+      {
+	
+	if (xyzmax != -424242)
+	{
+	  xmax1 = xyzmax;
+	  ymax1 = xyzmax;
+	  zmax1 = xyzmax;
+	}
+	else
+	{
+	  xmax1 = xmax;
+	  ymax1 = ymax;
+	  zmax1 = zmax;
+	}
+	
+	if (xyzmin != -424242)
+	{
+	  xmin1 = xyzmin;
+	  ymin1 = xyzmin;
+	  zmin1 = xyzmin;
+	}
+	else
+	{
+	  xmin1 = xmin;
+	  ymin1 = ymin;
+	  zmin1 = zmin;
+	}
 
-    }
-
+      }
 
-    
-    /* Grid spacing on coarsest grid */
-    /* TODO: Put the coordinates into arrays, and loop over the
-       dimensions.  That gets ride of all the triplicated code. */
-    if (domainsym[0]) 
-    {
-      if (cntstag[0]) 
+      /* Grid spacing on coarsest grid */
+      /* TODO: Put the coordinates into arrays, and loop over the
+	 dimensions.  That gets ride of all the triplicated code. */
+      if (domainsym[0]) 
       {
-        *coarse_dx = xmax1 / (cctk_gsh[0] - cctk_nghostzones[0] - 0.5);
-        x_origin = - (cctk_nghostzones[0] - 0.5) * *coarse_dx * iconv;
+	if (cntstag[0]) 
+	{
+	  *coarse_dx = xmax1 / (cctk_gsh[0] - cctk_nghostzones[0] - 0.5);
+	  x_origin = - (cctk_nghostzones[0] - 0.5) * *coarse_dx * iconv;
+	} 
+	else 
+	{
+	  *coarse_dx = xmax1 / (cctk_gsh[0] - cctk_nghostzones[0] - 1);
+	  x_origin = -  cctk_nghostzones[0]        * *coarse_dx * iconv;
+	}
       } 
       else 
       {
-        *coarse_dx = xmax1 / (cctk_gsh[0] - cctk_nghostzones[0] - 1);
-        x_origin = -  cctk_nghostzones[0]        * *coarse_dx * iconv;
+	*coarse_dx = (xmax1 - xmin1) / max(cctk_gsh[0] - 1, 1);
+	x_origin = xmin1;
       }
-    } 
-    else 
-    {
-      *coarse_dx = (xmax1 - xmin1) / max(cctk_gsh[0] - 1, 1);
-      x_origin = xmin1;
-    }
 
-    if (domainsym[2]) 
-    {
-      if (cntstag[1]) 
+      if (domainsym[2]) 
       {
-        *coarse_dy = ymax1 / (cctk_gsh[1] - cctk_nghostzones[1] - 0.5);
-        y_origin = - (cctk_nghostzones[1] - 0.5) * *coarse_dy * iconv;
+	if (cntstag[1]) 
+	{
+	  *coarse_dy = ymax1 / (cctk_gsh[1] - cctk_nghostzones[1] - 0.5);
+	  y_origin = - (cctk_nghostzones[1] - 0.5) * *coarse_dy * iconv;
+	} 
+	else 
+	{
+	  *coarse_dy = ymax1 / (cctk_gsh[1] - cctk_nghostzones[1] - 1);
+	  y_origin = -  cctk_nghostzones[1]      * *coarse_dy * iconv;
+	}
       } 
       else 
       {
-        *coarse_dy = ymax1 / (cctk_gsh[1] - cctk_nghostzones[1] - 1);
-        y_origin = -  cctk_nghostzones[1]      * *coarse_dy * iconv;
+	*coarse_dy = (ymax1 - ymin1) / max(cctk_gsh[1] - 1, 1);
+	y_origin = ymin1;
       }
-    } 
-    else 
-    {
-      *coarse_dy = (ymax1 - ymin1) / max(cctk_gsh[1] - 1, 1);
-      y_origin = ymin1;
-    }
     
-    if (domainsym[4]) 
-    {
-      if (cntstag[2]) 
+      if (domainsym[4]) 
       {
-        *coarse_dz = zmax1 / (cctk_gsh[2] - cctk_nghostzones[2] - 0.5);
-        z_origin = - (cctk_nghostzones[2] - 0.5) * *coarse_dz * iconv;
+	if (cntstag[2]) 
+	{
+	  *coarse_dz = zmax1 / (cctk_gsh[2] - cctk_nghostzones[2] - 0.5);
+	  z_origin = - (cctk_nghostzones[2] - 0.5) * *coarse_dz * iconv;
+	} 
+	else 
+	{
+	  *coarse_dz = zmax1 / (cctk_gsh[2] - cctk_nghostzones[2] - 1);
+	  z_origin = -  cctk_nghostzones[2]      * *coarse_dz * iconv;
+	}
       } 
       else 
       {
-        *coarse_dz = zmax1 / (cctk_gsh[2] - cctk_nghostzones[2] - 1);
-        z_origin = -  cctk_nghostzones[2]      * *coarse_dz * iconv;
+	*coarse_dz = (zmax1 - zmin1) / max(cctk_gsh[2] - 1, 1);
+	z_origin = zmin1;
       }
-    } 
-    else 
-    {
-      *coarse_dz = (zmax1 - zmin1) / max(cctk_gsh[2] - 1, 1);
-      z_origin = zmin1;
+      
+      /* dx,dy,dz on the grid we are on */
+      this_dx = *coarse_dx * iconv;
+      this_dy = *coarse_dy * iconv;
+      this_dz = *coarse_dz * iconv;
+      
     }
-    
-    /* dx,dy,dz on the grid we are on */
-    this_dx = *coarse_dx * iconv;
-    this_dy = *coarse_dy * iconv;
-    this_dz = *coarse_dz * iconv;
 
-  }
-
-  /**************************************************************
-   *  BYSPACING
-   *
-   *  User gives: grid spacing on the coarsest GH and
-   *              the number of gridpoints on the coarsest GH
-   *
-   **************************************************************/
-
-  else if (CCTK_Equals(type,"byspacing")) 
-  {
+    /**************************************************************
+     *  BYSPACING
+     *
+     *  User gives: grid spacing on the coarsest GH and
+     *              the number of gridpoints on the coarsest GH
+     *
+     **************************************************************/
 
-    /* Dx, Dy, Dx on the coarsest grid */
-
-    if (dxyz > 0) 
-    {
-      *coarse_dx = dxyz;
-      *coarse_dy = dxyz;
-      *coarse_dz = dxyz;
-    }
-    else
+    else if (CCTK_Equals(type,"byspacing")) 
     {
-      *coarse_dx = dx;
-      *coarse_dy = dy;
-      *coarse_dz = dz;
-    }
+      
+      /* Dx, Dy, Dx on the coarsest grid */
+      
+      if (dxyz > 0) 
+      {
+	*coarse_dx = dxyz;
+	*coarse_dy = dxyz;
+	*coarse_dz = dxyz;
+      }
+      else
+      {
+	*coarse_dx = dx;
+	*coarse_dy = dy;
+	*coarse_dz = dz;
+      }
 
-    /* dx, dy, dz on the grid we are on */
-    this_dx = *coarse_dx * iconv;
-    this_dy = *coarse_dy * iconv;
-    this_dz = *coarse_dz * iconv;
+      /* dx, dy, dz on the grid we are on */
+      this_dx = *coarse_dx * iconv;
+      this_dy = *coarse_dy * iconv;
+      this_dz = *coarse_dz * iconv;
     
+      /* Set minimum values of coordinates */
+      if (domainsym[0]) 
+      {
+	x_origin = (- cctk_nghostzones[0] + cntstag[0] * 0.5) * this_dx;
+      } 
+      else 
+      {
+	x_origin = - 0.5 * (cctk_gsh[0]-1 - cntstag[0] * (cctk_gsh[0])%2)
+	  * this_dx;
+      }
+      
+      if (domainsym[2]) 
+      {
+	y_origin = (- cctk_nghostzones[1] + cntstag[1] * 0.5) * this_dy;
+      } 
+      else 
+      {
+	y_origin = - 0.5 * (cctk_gsh[1]-1 - cntstag[1] * (cctk_gsh[1])%2)
+	  * this_dy;
+      }
+      
+      if (domainsym[4]) 
+      {
+	z_origin = (- cctk_nghostzones[2] + cntstag[2] * 0.5) * this_dz;
+      } 
+      else 
+      {
+	z_origin = - 0.5 * (cctk_gsh[2]-1 - cntstag[2] * (cctk_gsh[2])%2)
+	  * this_dz;
+      }
+    }
+  } 
+  else 
+  { 
+    /* if (not coarsest refinement level) */
     
+    /* Use the already calculated coordinate ranges for all but the
+       coarsest levels */
     
-    /* Set minimum values of coordinates */
-    if (domainsym[0]) 
-    {
-      x_origin = (- cctk_nghostzones[0] + cntstag[0] * 0.5) * this_dx;
-    } 
-    else 
-    {
-      x_origin = - 0.5 * (cctk_gsh[0]-1 - cntstag[0] * (cctk_gsh[0])%2)
-        * this_dx;
-    }
+    this_dx = cctk_delta_space[0] / cctk_levfac[0];
+    this_dy = cctk_delta_space[1] / cctk_levfac[1];
+    this_dz = cctk_delta_space[2] / cctk_levfac[2];
     
-    if (domainsym[2]) 
-    {
-      y_origin = (- cctk_nghostzones[1] + cntstag[1] * 0.5) * this_dy;
-    } 
-    else 
-    {
-      y_origin = - 0.5 * (cctk_gsh[1]-1 - cntstag[1] * (cctk_gsh[1])%2)
-        * this_dy;
-    }
+    x_origin = cctk_origin_space[0];
+    y_origin = cctk_origin_space[1];
+    z_origin = cctk_origin_space[2];
     
-    if (domainsym[4]) 
-    {
-      z_origin = (- cctk_nghostzones[2] + cntstag[2] * 0.5) * this_dz;
-    } 
-    else 
-    {
-      z_origin = - 0.5 * (cctk_gsh[2]-1 - cntstag[2] * (cctk_gsh[2])%2)
-       * this_dz;
-    }
-  }
+  } /* if (not coarsest refinement level) */
+
+
 
   /* Set spatial coordinates */
 
@@ -326,51 +346,83 @@ void CartGrid3D(CCTK_ARGUMENTS)
     }
   }
   
-  cctk_delta_space[0] = this_dx;
-  cctk_delta_space[1] = this_dy;
-  cctk_delta_space[2] = this_dz;
+  if (cctk_levfac[0] == 1 && cctk_levfac[1] == 1 && cctk_levfac[2] == 1) 
+  {
+
+    /* Register the coordinate ranges only for the coarsest level */
+    
+    cctk_delta_space[0] = this_dx;
+    cctk_delta_space[1] = this_dy;
+    cctk_delta_space[2] = this_dz;
+    
+    cctk_origin_space[0] = x_origin;
+    cctk_origin_space[1] = y_origin;
+    cctk_origin_space[2] = z_origin;
+    
+    lowerx = x_origin;
+    upperx = x_origin+this_dx*(cctk_gsh[0]-1); 
+    ierr = CCTK_CoordRegisterRange(cctkGH,lowerx,upperx,-1,"x","cart3d");
+    if (ierr < 0)
+    {
+      CCTK_WARN(0,"Failed to register x-coordinate computational range");
+    }
+
+    ierr = CCTK_CoordRegisterRangePhysIndex(cctkGH,lowerxi,upperxi,-1,"x","cart3d");
+    if (ierr < 0)
+    {
+      CCTK_WARN(0,"Failed to register x-coordinate physical range");
+    }
+
+    lowery = y_origin;
+    uppery = y_origin+this_dy*(cctk_gsh[1]-1); 
+    ierr = CCTK_CoordRegisterRange(cctkGH,lowery,uppery,-1,"y","cart3d");
+    if (ierr < 0)
+    {
+      CCTK_WARN(0,"Failed to register y-coordinate computational range");
+    }
+    ierr = CCTK_CoordRegisterRangePhysIndex(cctkGH,loweryi,upperyi,-1,"y","cart3d");
+    if (ierr < 0)
+    {
+      CCTK_WARN(0,"Failed to register y-coordinate physical range");
+    }
+
+    lowerz = z_origin;
+    upperz = z_origin+this_dz*(cctk_gsh[2]-1); 
+    ierr = CCTK_CoordRegisterRange(cctkGH,lowerz,upperz,-1,"z","cart3d");
+    if (ierr < 0)
+    {
+      CCTK_WARN(0,"Failed to register z-coordinate computational range");
+    }
+    ierr = CCTK_CoordRegisterRangePhysIndex(cctkGH,lowerzi,upperzi,-1,"z","cart3d");
+    if (ierr < 0)
+    {
+      CCTK_WARN(0,"Failed to register z-coordinate physical range");
+    }
   
-  cctk_origin_space[0] = x_origin;
-  cctk_origin_space[1] = y_origin;
-  cctk_origin_space[2] = z_origin;
+    CCTK_INFO("Grid Spacings:");
+    sprintf(infoline," %s%12.7e  %s%12.7e  %s%12.7e  ",  
+	    "dx=>",cctk_delta_space[0],
+	    "dy=>",cctk_delta_space[1],
+	    "dz=>",cctk_delta_space[2]);
+    CCTK_INFO(infoline);
   
-  lowerx = x_origin;
-  upperx = x_origin+this_dx*(cctk_gsh[0]-1); 
-  CCTK_CoordRegisterRange(cctkGH,lowerx,upperx,-1,"x","cart3d");
-  CCTK_CoordRegisterRangePhysIndex(cctkGH,lowerxi,upperxi,-1,"x","cart3d");
-
-  lowery = y_origin;
-  uppery = y_origin+this_dy*(cctk_gsh[1]-1); 
-  CCTK_CoordRegisterRange(cctkGH,lowery,uppery,-1,"y","cart3d");
-  CCTK_CoordRegisterRangePhysIndex(cctkGH,loweryi,upperyi,-1,"y","cart3d");
-
-  lowerz = z_origin;
-  upperz = z_origin+this_dz*(cctk_gsh[2]-1); 
-  CCTK_CoordRegisterRange(cctkGH,lowerz,upperz,-1,"z","cart3d");
-  CCTK_CoordRegisterRangePhysIndex(cctkGH,lowerzi,upperzi,-1,"z","cart3d");
+    CCTK_INFO("Computational Coordinates:");
+    sprintf(infoline," %s[%6.3f,%6.3f]  %s[%6.3f,%6.3f]  %s[%6.3f,%6.3f]  ", 
+	    "x=>",lowerx,upperx,
+	    "y=>",lowery,uppery,
+	    "z=>",lowerz,upperz);
+    CCTK_INFO(infoline);
+
+    CCTK_INFO("Indices of Physical Coordinates:");
+    sprintf(infoline," %s[%d,%d]  %s[%d,%d]  %s[%d,%d]  ", 
+	    "x=>",lowerxi,upperxi,
+	    "y=>",loweryi,upperyi,
+	    "z=>",lowerzi,upperzi);
+    CCTK_INFO(infoline);
+  
+  } /* if (coarsest refinement level) */
   
-  CCTK_INFO("Grid Spacings:");
-  sprintf(infoline," %s%12.7e  %s%12.7e  %s%12.7e  ",  
-                   "dx=>",cctk_delta_space[0],
-                   "dy=>",cctk_delta_space[1],
-                   "dz=>",cctk_delta_space[2]);
-  CCTK_INFO(infoline);
   
-  CCTK_INFO("Computational Coordinates:");
-  sprintf(infoline," %s[%6.3f,%6.3f]  %s[%6.3f,%6.3f]  %s[%6.3f,%6.3f]  ", 
-                   "x=>",lowerx,upperx,
-                   "y=>",lowery,uppery,
-                   "z=>",lowerz,upperz);
-  CCTK_INFO(infoline);
-
-  CCTK_INFO("Indices of Physical Coordinates:");
-  sprintf(infoline," %s[%d,%d]  %s[%d,%d]  %s[%d,%d]  ", 
-                   "x=>",lowerxi,upperxi,
- 	           "y=>",loweryi,upperyi,
-	           "z=>",lowerzi,upperzi);
-  CCTK_INFO(infoline);
-
-
 #ifdef CCTK_DEBUG
   printf("\n");
   printf("CartGrid3D\n");