Check whether the symmetry boundary is indeed at the origin (or a

user-specified location).

Offer to poison the symmetry boundary points before the symmetry is
applied, to check whether the symmetry boundary condision is indeed
applied to all boundary points.

Add support for complex grid functions.

Correct a small memory leak.


git-svn-id: http://svn.cactuscode.org/arrangements/CactusNumerical/RotatingSymmetry180/trunk@33 20f44201-0f4f-0410-9130-e5fc2714a787

3 files changed, 179 insertions, 7 deletions

diff --git a/interface.ccl b/interface.ccl
index 3b61ff1..a252cfb 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -71,3 +71,16 @@ CCTK_INT FUNCTION                            \
          CCTK_INT ARRAY OUT table_handles,   \
          CCTK_STRING IN  bc_name)
 REQUIRES FUNCTION Boundary_SelectedGVs
+
+
+
+CCTK_INT FUNCTION GetDomainSpecification \
+  (CCTK_INT IN size, \
+   CCTK_REAL OUT ARRAY physical_min, \
+   CCTK_REAL OUT ARRAY physical_max, \
+   CCTK_REAL OUT ARRAY interior_min, \
+   CCTK_REAL OUT ARRAY interior_max, \
+   CCTK_REAL OUT ARRAY exterior_min, \
+   CCTK_REAL OUT ARRAY exterior_max, \
+   CCTK_REAL OUT ARRAY spacing)
+USES FUNCTION GetDomainSpecification
diff --git a/param.ccl b/param.ccl
index 00f26e8..c67ca31 100644
--- a/param.ccl
+++ b/param.ccl
@@ -4,3 +4,30 @@
 BOOLEAN verbose "Produce screen output while applying boundary conditions"
 {
 } "no"
+
+
+
+BOOLEAN use_coordbase "Get symmetry boundary location from CoordBase (and ignore symmetry_boundary_* parameters)"
+{
+} "no"
+
+CCTK_REAL symmetry_boundary_x "x-coordinate of symmetry boundary (ignored if use_coordbase is set)"
+{
+  *:* :: ""
+} 0.0
+
+CCTK_REAL symmetry_boundary_y "y-coordinate of symmetry boundary (ignored if use_coordbase is set)"
+{
+  *:* :: ""
+} 0.0
+
+
+
+BOOLEAN poison_boundaries "Fill the symmetry boundary with a poison value before the symmetry is applied, and check afterwards whether it has been overwritten"
+{
+} "no"
+
+CCTK_REAL poison_value "Poison value (must not occur naturally)"
+{
+  *:* :: ""
+} -424242.0
diff --git a/src/rotatingsymmetry180.c b/src/rotatingsymmetry180.c
index 7404770..fdb69e9 100644
--- a/src/rotatingsymmetry180.c
+++ b/src/rotatingsymmetry180.c
@@ -39,6 +39,7 @@ int BndRot180VI (cGH const * restrict const cctkGH,
   int fake_bbox[6];
   
   CCTK_REAL x0[3], dx[3];
+  CCTK_REAL symbnd[3];          /* location of symmetry boundary */
   CCTK_REAL origin[3], dorigin[3];
   int avoid_origin[3], iorigin[3];
   int offset[3];                /* offset 0..1 due to avoid_origin */
@@ -274,11 +275,36 @@ int BndRot180VI (cGH const * restrict const cctkGH,
   alldirs[0] = 0;
   alldirs[1] = 1;
   
+  /* Find location of symmetry boundary */
+  if (use_coordbase) {
+    CCTK_INT const size = 3;
+    CCTK_REAL physical_min[3];
+    CCTK_REAL physical_max[3];
+    CCTK_REAL interior_min[3];
+    CCTK_REAL interior_max[3];
+    CCTK_REAL exterior_min[3];
+    CCTK_REAL exterior_max[3];
+    CCTK_REAL spacing;
+    GetDomainSpecification
+      (size,
+       physical_min, physical_max,
+       interior_min, interior_max,
+       exterior_min, exterior_max,
+       & spacing);
+    symbnd[0] = physical_min[0];
+    symbnd[1] = physical_min[1];
+    symbnd[2] = physical_min[2];
+  } else {
+    symbnd[0] = symmetry_boundary_x;
+    symbnd[1] = symmetry_boundary_y;
+    symbnd[2] = 0;              /* unused */
+  }
+  
   /* Find grid point that corresponds to the origin */
   for (q=0; q<2; ++q) {
     d = alldirs[q];
-    /* x0 + dx * origin == 0 */
-    origin[d] = - x0[d] / dx[d];
+    /* x0 + dx * origin == symbnd */
+    origin[d] = (symbnd[d] - x0[d]) / dx[d];
     dorigin[d] = origin[d] - floor(origin[d]);
     if (fabs(dorigin[d]) < 1.0e-6 || fabs(dorigin[d] - 1.0) < 1.0e-6) {
       avoid_origin[d] = 0;
@@ -307,6 +333,15 @@ int BndRot180VI (cGH const * restrict const cctkGH,
                   "xyz"[otherdir]);
     }
     
+    if (iorigin[dir] != data.nghostzones[dir]) {
+      CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
+                  "The group \"%s\" has in the %c-direction %d symmetry points (grid points outside of the symmetry boundary).  "
+                  "This is not equal to the number of ghost zones, which is %d.  "
+                  "The number of symmetry points must be equal to the number of ghost zones.",
+                  CCTK_GroupNameFromVarI(vis[var]), "xyz"[dir],
+                  iorigin[dir], data.nghostzones[dir]);
+    }
+    
     if (data.gsh[dir] < 2*data.nghostzones[dir] + offset[dir]) {
       assert (nvars > 0);
       CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
@@ -318,6 +353,47 @@ int BndRot180VI (cGH const * restrict const cctkGH,
                   2*data.nghostzones[dir] + offset[dir]);
     }
     
+    if (poison_boundaries) {
+      /* poison destination grid points */
+      if (cctkGH->cctk_bbox[2*dir]) {
+        
+        int imin[3], imax[3];
+        int i, j, k;
+        for (d=0; d<3; ++d) {
+          imin[d] = xferinfo[d].dst.off;
+          imax[d] = xferinfo[d].dst.off + xferinfo[d].dst.len;
+          imin[d] -= cctk_lbnd[d];
+          imax[d] -= cctk_lbnd[d];
+          if (imin[d] < 0) imin[d] = 0;
+          if (imax[d] >= cctk_lsh[d]) imax[d] = cctk_lsh[d];
+        }
+        
+        var = 0;
+        
+        assert (group.dim == 3);
+        switch (group.vartype) {
+        case CCTK_VARIABLE_INT:
+          /* do nothing */
+          break;
+        case CCTK_VARIABLE_REAL:
+          for (k=imin[0]; k<imax[2]; ++k) {
+            for (j=imin[1]; j<imax[1]; ++j) {
+              for (i=imin[2]; i<imax[0]; ++i) {
+                const int ind = CCTK_GFINDEX3D(cctkGH, i, j, k);
+                ((CCTK_REAL *) varptrs[var]) [ind] = poison_value;
+              }
+            }
+          }
+          break;
+        case CCTK_VARIABLE_COMPLEX:
+          /* do nothing */
+          break;
+        default:
+          assert (0);
+        } /* switch grouptype */
+      } /* if bbox */
+    } /* if poison_boundaries */
+    
     /* Allocate slab transfer description */
     xferinfo = malloc(group.dim * sizeof *xferinfo);
     assert (xferinfo);
@@ -361,6 +437,48 @@ int BndRot180VI (cGH const * restrict const cctkGH,
        nvars, vartypes, varptrs, vartypes, varptrs);
     assert (!ierr);
     
+    if (poison_boundaries) {
+      /* check destination grid points for poison */
+      if (cctkGH->cctk_bbox[2*dir]) {
+        
+        int imin[3], imax[3];
+        int i, j, k;
+        for (d=0; d<3; ++d) {
+          imin[d] = xferinfo[d].dst.off;
+          imax[d] = xferinfo[d].dst.off + xferinfo[d].dst.len;
+          imin[d] -= cctk_lbnd[d];
+          imax[d] -= cctk_lbnd[d];
+          if (imin[d] < 0) imin[d] = 0;
+          if (imax[d] >= cctk_lsh[d]) imax[d] = cctk_lsh[d];
+        }
+        
+        var = 0;
+        
+        assert (group.dim == 3);
+        switch (group.vartype) {
+        case CCTK_VARIABLE_INT:
+          /* do nothing */
+          break;
+        case CCTK_VARIABLE_REAL:
+          for (k=imin[0]; k<imax[2]; ++k) {
+            for (j=imin[1]; j<imax[1]; ++j) {
+              for (i=imin[2]; i<imax[0]; ++i) {
+                const int ind = CCTK_GFINDEX3D(cctkGH, i, j, k);
+                CCTK_REAL const val = ((CCTK_REAL const *) varptrs[var]) [ind];
+                assert (val != poison_value);
+              }
+            }
+          }
+          break;
+        case CCTK_VARIABLE_COMPLEX:
+          /* do nothing */
+          break;
+        default:
+          assert (0);
+        } /* switch grouptype */
+      } /* if bbox */
+    } /* if poison_boundaries */
+    
     /* take parity into account */
     if (cctkGH->cctk_bbox[2*dir]) {
       for (var=0; var<nvars; ++var) {
@@ -401,17 +519,31 @@ int BndRot180VI (cGH const * restrict const cctkGH,
               }
             }
             break;
+          case CCTK_VARIABLE_COMPLEX:
+            {
+              CCTK_COMPLEX const czero = CCTK_Cmplx (0.0, 0.0);
+              for (k=imin[0]; k<imax[2]; ++k) {
+                for (j=imin[1]; j<imax[1]; ++j) {
+                  for (i=imin[2]; i<imax[0]; ++i) {
+                    const int ind = CCTK_GFINDEX3D(cctkGH, i, j, k);
+                    CCTK_COMPLEX * restrict const ptr =
+                      & ((CCTK_COMPLEX *) varptrs[var]) [ind];
+                    * ptr = CCTK_CmplxSub (czero, * ptr);
+                  }
+                }
+              }
+            }
+            break;
           default:
             assert (0);
           } /* switch grouptype */
         }
       } /* for var */
-      
-      free (xferinfo);
-      
-    } /* if global_bbox */
+    } /* if bbox */
+    
+    free (xferinfo);
     
-  } /* if bbox */
+  } /* if global_bbox */
   
   free (gis);
   free (varptrs);