Changed the way I handle boundaries on the local grid. Should now

work for full, bitant, quadrant and octant mode for the event horizon
evolution and for re-parametrization with euler pde. Still need testing though.


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/EHFinder/trunk@13 2a26948c-0e4f-0410-aee8-f1d3e353619c

7 files changed, 128 insertions, 16 deletions

diff --git a/param.ccl b/param.ccl
index 8196afb..02d3961 100644
--- a/param.ccl
+++ b/param.ccl
@@ -71,6 +71,11 @@ KEYWORD re_param_int_method "Integration method in pde re-parametrization"
   "rk2"    :: "Second order Runge-Kutta scheme"
 } "rk2"
 
+INT re_param_max_iter "maximum number of iteration in the re-parametrization"
+{
+(0:        :: "Positive please"
+} 800
+
 KEYWORD pde_differences "Type of finite diffencing used in pde re-parametrization"
 {
   "centered" :: "Use 2nd order centered differences except at the boundaries"
@@ -86,3 +91,5 @@ INT reparametrize_every "Re-parametrize every"
 shares: grid
 
 USES KEYWORD domain
+USES KEYWORD quadrant_direction
+USES KEYWORD bitant_plane
diff --git a/schedule.ccl b/schedule.ccl
index e19e0cc..e728618 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -29,6 +29,12 @@ schedule EHFinder_SetMask at CCTK_POSTINITIAL after EHFinder_MaskInit
   LANG: Fortran
 } "Set the mask"
 
+schedule  EHFinder_ApplySym at CCTK_POSTINITIAL after EHFinder_MaskInit
+{
+  LANG: Fortran
+  SYNC: level_set
+} "apply symmetry boundaries"
+
 #schedule GROUP EHFinder_ReParamInit at CCTK_POSTINITIAL after EHFinder_MaskInit
 #{
 #  LANG: Fortran
@@ -85,9 +91,7 @@ if (CCTK_Equals(mode,"test_reparam"))
 
 schedule GROUP EHFinder_PostStep in MoL_PostStep
 {
-  LANG: Fortran
-#  SYNC: level_set
-} "Schedule syncing of level set"
+} "Schedule group for symmetry boundaries and syncing"
 
 schedule EHFinder_ApplySym in EHFinder_PostStep
 {
@@ -111,11 +115,21 @@ if (!CCTK_Equals(mode,"test_reparam"))
 
     if (CCTK_Equals(re_param_int_method,"euler"))
     {
-      schedule EHFinder_ReParametrize4 in EHFinder_ReParametrize AFTER EHFinder_ReParametrize1 WHILE ehfinder::re_param_control
+      schedule GROUP Euler_ReParametrize in EHFinder_ReParametrize AFTER EHFinder_ReParametrize1 WHILE ehfinder::re_param_control
+      {
+      } "Schedule group for Euler re-parametrization"
+
+      schedule EHFinder_ReParametrize4 in Euler_ReParametrize 
       {
         LANG: Fortran
         SYNC: level_set
       } "Euler scheme"
+
+      schedule  EHFinder_ApplySym in Euler_ReParametrize after EHFinder_ReParametrize4
+      {
+        LANG: Fortran
+        SYNC: level_set
+      } "apply symmetry boundaries"
     }
 
     if (CCTK_Equals(re_param_int_method,"rk2"))
diff --git a/src/EHFinder_ReParametrize.F90 b/src/EHFinder_ReParametrize.F90
index 39e0ee5..149d624 100644
--- a/src/EHFinder_ReParametrize.F90
+++ b/src/EHFinder_ReParametrize.F90
@@ -251,9 +251,9 @@ subroutine EHFinder_ReParametrize4(CCTK_ARGUMENTS)
 
   maxdfloc = zero
   mindfloc = 1.0d23
-  do k = 1, nz
-    do j = 1, ny
-      do i =1, nx
+  do k = kzl, kzr
+    do j = jyl, jyr
+      do i = ixl, ixr
         if ( eh_mask(i,j,k) .ge. 0 ) then
           maxdfloc = max ( maxdfloc, abs(dfsq(i,j,k)) )
           mindfloc = min ( mindfloc, abs(dfsq(i,j,k)) )
diff --git a/src/EHFinder_Sources.F90 b/src/EHFinder_Sources.F90
index 6e6b912..1ccfee2 100644
--- a/src/EHFinder_Sources.F90
+++ b/src/EHFinder_Sources.F90
@@ -26,6 +26,8 @@ subroutine EHFinder_Sources(CCTK_ARGUMENTS)
 !  print*,'Entered Sources'
 ! calculate 1/(2*delta) in each direction
 
+#include "include/physical_part.h"
+
   idx = half / cctk_delta_space(1)
   idy = half / cctk_delta_space(2)
   idz = half / cctk_delta_space(3)
@@ -70,9 +72,9 @@ subroutine EHFinder_Sources(CCTK_ARGUMENTS)
 !    end where
 !  end where
 
-  do k = 1, nz
-    do j = 1, ny
-      do i = 1, nx
+  do k = kzl, kzr
+    do j = jyl, jyr
+      do i = ixl, ixr
         if ( eh_mask(i,j,k) .ge. 0 ) then
           if ( ( .not. btest(eh_mask(i,j,k),0) ) .and. &
                ( .not. btest(eh_mask(i,j,k),1) ) ) then
@@ -247,9 +249,10 @@ subroutine EHFinder_Sources(CCTK_ARGUMENTS)
 !  do k = 1+iz1, nz-iz2
 !    do j = 1+iy1, ny-iy2
 !      do i = 1+ix1, nx-ix2
-  do k = 1, nz
-    do j = 1, ny
-      do i = 1, nx
+
+  do k = kzl, kzr
+    do j = jyl, jyr
+      do i = ixl, ixr
         if ( eh_mask(i,j,k) .ge. 0 ) then
           alp2 = alp(i,j,k)**2
           g4tt = -one
diff --git a/src/EHFinder_mod.F90 b/src/EHFinder_mod.F90
index 2c88ee5..3979003 100644
--- a/src/EHFinder_mod.F90
+++ b/src/EHFinder_mod.F90
@@ -8,6 +8,7 @@ module EHFinder_mod
   CCTK_INT, dimension(0:5), parameter :: jy = (/ 0, 0, -1, 1, 0, 0 /)
   CCTK_INT, dimension(0:5), parameter :: kz = (/ 0, 0, 0, 0, -1, 1 /)
   CCTK_INT, dimension(0:5), parameter :: ll = (/ 1, 2, 4, 8, 16, 32 /)
+  CCTK_INT :: ixl, ixr, jyl, jyr, kzl, kzr
 !  CCTK_INT, parameter :: ixl = 1
 !  CCTK_INT, parameter :: ixr = 2
 !  CCTK_INT, parameter :: iyl = 4
diff --git a/src/include/physical_part.h b/src/include/physical_part.h
new file mode 100644
index 0000000..44166ea
--- /dev/null
+++ b/src/include/physical_part.h
@@ -0,0 +1,78 @@
+! Here I figure out the extension of the physical part of a chunk of a grid.
+
+! nx, ny and nz contains the full size of the chunk.
+
+nx = cctk_lsh(1)
+ny = cctk_lsh(2)
+nz = cctk_lsh(3)
+
+! ngx, ngy, ngz contains the size of any ghostzones.
+
+ngx = cctk_nghostzones(1)
+ngy = cctk_nghostzones(2)
+ngz = cctk_nghostzones(3)
+
+! ixl and ixr is set to the minimum and maximum local cell index in the
+! x-direction. The corresponding is done for the y- and z-directions.
+
+ixl = 1
+ixr = nx
+jyl = 1
+jyr = ny
+kzl = 1
+kzr = nz
+
+! I check to see if we have any processor boundaries. If we have, adjust the
+! minimum and maximum cell index appropriately with the size of the ghostzone.
+
+if ( cctk_bbox(1) .eq. 0 ) ixl = ixl + ngx
+if ( cctk_bbox(2) .eq. 0 ) ixr = ixr - ngx
+
+if ( cctk_bbox(3) .eq. 0 ) jyl = jyl + ngy
+if ( cctk_bbox(4) .eq. 0 ) jyr = jyr - ngy
+
+if ( cctk_bbox(5) .eq. 0 ) kzl = kzl + ngz
+if ( cctk_bbox(6) .eq. 0 ) kzr = kzr - ngz
+
+! In octant I check to see if the chunk of grid contains any symmetry
+! boundaries. This is done by checking if the lowest index of the local grid
+! as seen on the global grid is equal to zero. If this is the case adjust the
+! minimum and maximum cell index appropriately with the size of the ghostzone.
+
+if ( CCTK_Equals(domain,'octant') ) then
+  if (cctk_lbnd(1) .eq. 0 ) ixl = ixl + ngx
+  if (cctk_lbnd(2) .eq. 0 ) jyl = jyl + ngy
+  if (cctk_lbnd(3) .eq. 0 ) kzl = kzl + ngz
+end if
+
+! In quadrant I do the same, except here I have to check in which direction
+! the quadrant is directed.
+
+if ( CCTK_Equals(domain,'quadrant') ) then
+  if ( CCTK_Equals(quadrant_direction,'x') ) then
+    if (cctk_lbnd(2) .eq. 0 ) jyl = jyl + ngy
+    if (cctk_lbnd(3) .eq. 0 ) kzl = kzl + ngz
+  end if
+  if ( CCTK_Equals(quadrant_direction,'y') ) then
+    if (cctk_lbnd(1) .eq. 0 ) ixl = ixl + ngx
+    if (cctk_lbnd(3) .eq. 0 ) kzl = kzl + ngz
+  end if
+  if ( CCTK_Equals(quadrant_direction,'z') ) then
+    if (cctk_lbnd(1) .eq. 0 ) ixl = ixl + ngx
+    if (cctk_lbnd(2) .eq. 0 ) jyl = jyl + ngy
+  end if
+end if
+
+! Ditto for bitant. Here I have to check which plane is the symmetry-plane.
+
+if ( CCTK_Equals(domain,'bitant') ) then
+  if ( CCTK_Equals(bitant_plane,'xy') ) then
+    if (cctk_lbnd(3) .eq. 0 ) kzl = kzl + ngz
+  endif
+  if ( CCTK_Equals(bitant_plane,'xz') ) then
+    if (cctk_lbnd(2) .eq. 0 ) jyl = jyl + ngy
+  endif
+  if ( CCTK_Equals(bitant_plane,'yz') ) then
+    if (cctk_lbnd(1) .eq. 0 ) ixl = ixl + ngx
+  endif
+end if
diff --git a/src/include/upwind_second.h b/src/include/upwind_second.h
index ff13ed0..84d41da 100644
--- a/src/include/upwind_second.h
+++ b/src/include/upwind_second.h
@@ -1,10 +1,19 @@
+# include "physical_part.h"
+
 idx = half / cctk_delta_space(1)
 idy = half / cctk_delta_space(2)
 idz = half / cctk_delta_space(3)
 
-do k = 1, nz
-  do j = 1, ny
-    do i = 1, nx
+if ( cctk_bbox(1) .eq. 0 ) eh_mask(1:ngx,:,:) = -2
+if ( cctk_bbox(2) .eq. 0 ) eh_mask(nx-ngx+1:nx,:,:) = -2
+if ( cctk_bbox(3) .eq. 0 ) eh_mask(:,1:ngy,:) = -2
+if ( cctk_bbox(4) .eq. 0 ) eh_mask(:,ny-ngy+1:ny,:) = -2
+if ( cctk_bbox(5) .eq. 0 ) eh_mask(:,:,1:ngz) = -2
+if ( cctk_bbox(6) .eq. 0 ) eh_mask(:,:,nz-ngz+1:nz) = -2
+
+do k = kzl, kzr
+  do j = jyl, jyr
+    do i = ixl, ixr
       if ( eh_mask(i,j,k) .ge. 0 ) then
         if ( ( .not. btest(eh_mask(i,j,k),0) ) .and. &
              ( .not. btest(eh_mask(i,j,k),1) ) ) then