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

1 files changed, 139 insertions, 7 deletions

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);