Added support for setting up generators on ellipsoidal initial data.

Added support for evolving the generators independently (used for testing).

Added support for detecting a peanut shaped excision region and take the
appropriate precautions.


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

4 files changed, 247 insertions, 39 deletions

diff --git a/param.ccl b/param.ccl
index 45e7cfb..e20eeeb 100644
--- a/param.ccl
+++ b/param.ccl
@@ -118,6 +118,7 @@ KEYWORD mode "Mode of operation"
   "normal"  :: "Find event horizons"
   "test_reparam" :: "Test reparametrization routines with MoL"
   "analysis" :: "Provide storage for f without evolving"
+  "generator" :: "Provide storage for f and initialize without evolving"
 } "normal"
 
 KEYWORD upwind_type "Type of upwinding used in evolving the ehfinder equations"
@@ -226,7 +227,7 @@ CCTK_INT number_of_generators "How many generators should be evolved"
 
 KEYWORD generator_distribution "What initial distribution should be used"
 {
-  "line"	:: "Put the generatos on a line on the x-axis"
+  "line"	:: "Put the generators on a line in the xz-plane"
 } "line"
 
 STRING surface_interpolator "What interpolator should be used to locate the surface"
diff --git a/schedule.ccl b/schedule.ccl
index 76ae893..e6b3fd4 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -503,16 +503,56 @@ if (CCTK_Equals(mode,"normal"))
   {
     LANG: Fortran
     SYNC: eh_mask_all
-  } "Finish modifying the mask"
+  } "Find excision boundaries"
 
-  schedule  EHFinder_ApplySymMask in EHFinder_SetMask after EHFinder_SetMask2
+  schedule  EHFinder_ApplySymMask as ASM1 in EHFinder_SetMask after EHFinder_SetMask2
   {
     LANG: Fortran
   } "Apply symmetry boundaries"
  
+  schedule EHFinder_SetMask3 in EHFinder_SetMask after EHFinder_ASM1
+  {
+    LANG: Fortran
+    SYNC: eh_mask_all
+  } "Check to see if the mask needs to be modified"
+
+  schedule  EHFinder_ApplySymMask as ASM2 in EHFinder_SetMask after EHFinder_SetMask3
+  {
+    LANG: Fortran
+  } "Apply symmetry boundaries"
+ 
+  schedule EHFinder_SetMask2 as SM2 in EHFinder_SetMask after ASM2
+  {
+    LANG: Fortran
+    SYNC: eh_mask_all
+  } "Find excision boundaries"
 
+  schedule  EHFinder_ApplySymMask as ASM3 in EHFinder_SetMask after SM2
+  {
+    LANG: Fortran
+  } "Apply symmetry boundaries"
+ 
 }
 
+  if ( evolve_generators)
+  {
+    if (CCTK_Equals(generator_tracking_method,"interpolate_before"))
+    {
+      schedule EHFinder_Generator_Sources as EGS in MoL_CalcRHS
+      {
+        LANG: Fortran
+      } "Calculate the source terms for the generator evolution"
+    }
+    if (CCTK_Equals(generator_tracking_method,"interpolate_after"))
+    {
+      schedule EHFinder_Generator_Sources2 as EGS2 in MoL_CalcRHS after EHFinder_Sources
+      {
+        LANG: Fortran
+      } "Calculate the source terms for the generator evolution"
+    }
+  }
+
+
 }
 
 #schedule  EHFinder_ApplySym at CCTK_POSTINITIAL after EHFinder_MaskInit
diff --git a/src/EHFinder_Init.F90 b/src/EHFinder_Init.F90
index 3036a1d..2bf3d7d 100644
--- a/src/EHFinder_Init.F90
+++ b/src/EHFinder_Init.F90
@@ -20,7 +20,7 @@ subroutine EHFinder_Init_F(CCTK_ARGUMENTS)
   CCTK_REAL, dimension(3,3) :: txyz
   CCTK_REAL :: cosa, sina, cosb, sinb, cosc, sinc
   CCTK_REAL :: last_time
-  CCTK_REAL :: theta, dtheta, thetamin, thetamax
+  CCTK_REAL :: theta, dtheta, thetamin, thetamax, r_el
   CCTK_INT, dimension(1) :: lsh, lbnd
 
   ! Initialize the current_iteration variable.
@@ -35,6 +35,46 @@ subroutine EHFinder_Init_F(CCTK_ARGUMENTS)
 
   cctk_time = last_time
 
+  if ( evolve_generators .gt. 0 ) then
+
+    call CCTK_GrouplbndGN ( status, cctkGH, 1, lbnd, "ehfinder::generators" )
+    if ( status .lt. 0 ) then
+      call CCTK_WARN ( 0, "cannot get lower bounds for generator arrays" )
+    end if
+    call CCTK_GrouplshGN ( status, cctkGH, 1, lsh, "ehfinder::generators" )
+    if ( status .lt. 0 ) then
+      call CCTK_WARN ( 0, "cannot get local size for generator arrays" )
+    end if
+
+    if ( CCTK_EQUALS( generator_distribution, 'line' ) ) then
+
+      if ( CCTK_EQUALS( domain, 'full' ) ) then
+        thetamin = zero; thetamax = pi
+      else if ( CCTK_EQUALS( domain, 'bitant') ) then
+        if ( CCTK_EQUALS( bitant_plane, 'xy' ) ) then
+          thetamin = zero; thetamax = half * pi
+        else 
+          thetamin = zero; thetamax = pi
+        end if
+      else if ( CCTK_EQUALS( domain, 'quadrant' ) ) then
+        if ( CCTK_EQUALS( quadrant_direction, 'x' ) .or. &
+             CCTK_EQUALS( quadrant_direction, 'y' ) ) then
+          thetamin = zero; thetamax = half * pi
+        else
+          thetamin = zero; thetamax = pi
+        end if
+      else if ( CCTK_EQUALS( domain, 'octant' ) ) then
+        thetamin = zero; thetamax = half * pi
+      end if
+                   
+      if ( number_of_generators .eq. 1 ) then
+        theta = half * ( thetamax - thetamin ) + thetamin
+      else
+        dtheta = ( thetamax - thetamin ) / ( number_of_generators - 1 )
+      end if
+    end if
+  end if
+
   ! If a sphere is requested...
   if ( CCTK_EQUALS( initial_f, 'sphere' ) ) then
 
@@ -44,42 +84,7 @@ subroutine EHFinder_Init_F(CCTK_ARGUMENTS)
               ( y - translate_y )**2 + &
               ( z - translate_z )**2 ) - initial_rad
     if ( evolve_generators .gt. 0 ) then
-
-      call CCTK_GrouplbndGN ( status, cctkGH, 1, lbnd, "ehfinder::generators" )
-      if ( status .lt. 0 ) then
-        call CCTK_WARN ( 0, "cannot get lower bounds for generator arrays" )
-      end if
-      call CCTK_GrouplshGN ( status, cctkGH, 1, lsh, "ehfinder::generators" )
-      if ( status .lt. 0 ) then
-        call CCTK_WARN ( 0, "cannot get local size for generator arrays" )
-      end if
-
       if ( CCTK_EQUALS( generator_distribution, 'line' ) ) then
-
-        if ( CCTK_EQUALS( domain, 'full' ) ) then
-          thetamin = zero; thetamax = pi
-        else if ( CCTK_EQUALS( domain, 'bitant') ) then
-          if ( CCTK_EQUALS( bitant_plane, 'xy' ) ) then
-            thetamin = zero; thetamax = half * pi
-          else 
-            thetamin = zero; thetamax = pi
-          end if
-        else if ( CCTK_EQUALS( domain, 'quadrant' ) ) then
-          if ( CCTK_EQUALS( quadrant_direction, 'x' ) .or. &
-               CCTK_EQUALS( quadrant_direction, 'y' ) ) then
-            thetamin = zero; thetamax = half * pi
-          else
-            thetamin = zero; thetamax = pi
-          end if
-        else if ( CCTK_EQUALS( domain, 'octant' ) ) then
-          thetamin = zero; thetamax = half * pi
-        end if
-                   
-        if ( number_of_generators .eq. 1 ) then
-          theta = half * ( thetamax - thetamin ) + thetamin
-        else
-          dtheta = ( thetamax - thetamin ) / ( number_of_generators - 1 )
-        end if
         do i = 1, lsh(1)
            theta = thetamin + dtheta * ( i + lbnd(1) - 1 )
            xg(i) = initial_rad * sin(theta) + translate_x
@@ -130,6 +135,23 @@ subroutine EHFinder_Init_F(CCTK_ARGUMENTS)
         end do
       end do
     end do
+
+    if ( evolve_generators .gt. 0 ) then
+      if ( CCTK_EQUALS( generator_distribution, 'line' ) ) then
+        do i = 1, lsh(1)
+          theta = thetamin + dtheta * ( i + lbnd(1) - 1 )
+          r_el = sqrt ( one / ( sin(theta)**2 / initial_a**2 + &
+                         cos(theta)**2 / initial_c**2 ) )
+          xp(1) = r_el * sin(theta) + translate_x
+          xp(2) = translate_y
+          xp(3) = r_el * cos(theta) + translate_z
+          xpt = matmul ( txyz, xp )
+          xg(i) = xpt(1)
+          yg(i) = xpt(2)
+          zg(i) = xpt(3)
+        end do
+      end if
+    end if
   end if
 
   ! if an ovaloid of Cassini is requested...
diff --git a/src/EHFinder_SetMask.F90 b/src/EHFinder_SetMask.F90
index 21d7fc9..9a64ca0 100644
--- a/src/EHFinder_SetMask.F90
+++ b/src/EHFinder_SetMask.F90
@@ -570,3 +570,148 @@ subroutine EHFinder_SetMask2(CCTK_ARGUMENTS)
   end if
 
 end subroutine EHFinder_SetMask2
+
+
+subroutine EHFinder_SetMask3(CCTK_ARGUMENTS)
+
+  use EHFinder_mod
+
+  implicit none
+
+  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_FUNCTIONS
+
+  CCTK_INT :: i, j, k
+  logical active, mask_modified
+
+  active = .false.
+
+! If re-parametrization has just been done, we are in the process of
+! resetting the mask so we set active=.true.
+  if ( mask_first ) active = .true.
+  if ( CCTK_EQUALS(re_param_method,'approx') .or. &
+       CCTK_EQUALS(re_param_method,'mixed') ) then
+    if ( reparametrize_every_approx .gt. 0 ) then
+      if ( mod(cctk_iteration,reparametrize_every_approx) .eq. 0 ) then
+        active = .true.
+      end if
+    end if
+  end if
+  if ( CCTK_EQUALS(re_param_method,'pde') .or. &
+       CCTK_EQUALS(re_param_method,'mixed') ) then
+    if ( reparametrize_every_pde .gt. 0 ) then
+      if ( mod(cctk_iteration,reparametrize_every_pde) .eq. 0 ) active = .true.
+    end if
+  end if
+
+! If the reparametrization was not undone...
+  if ( active .and. .not.reparam_undone ) then
+
+    mask_modified = .false.
+# include "include/physical_part.h"
+
+!   Make sure we are not on the physical outer boundary.
+    if ( ixl .eq. 1 ) ixl = 2
+    if ( ixr .eq. nx ) ixr = nx - 1
+    if ( jyl .eq. 1 ) jyl = 2
+    if ( jyr .eq. ny ) jyr = ny - 1
+    if ( kzl .eq. 1 ) kzl = 2
+    if ( kzr .eq. nz ) kzr = nz - 1
+
+    tm_mask = eh_mask
+
+    do k = kzl, kzr
+      do j = jyl, jyr
+        do i = ixl, ixr
+          if ( eh_mask(i,j,k) .ge. 0 ) then
+            ! If we are on an excision boundary in the -x-direction...
+            if ( btest(eh_mask(i,j,k),0) ) then
+              ! If any of the two closest neighbours in the +x-direction is
+              ! excised we have to excise this point
+              if ( ( eh_mask(i+1,j,k) .eq. -1 ) .or. &
+                   ( eh_mask(i+2,j,k) .eq. -1 ) ) then
+                tm_mask(i,j,k) = -1
+                mask_modified = .true.
+                print*, i, j, k, 'x-'
+                print*, eh_mask(i:i+2,j,k)
+              end if
+            end if
+
+            ! If we are on an excision boundary in the +x-direction...
+            if ( btest(eh_mask(i,j,k),1) ) then
+              ! If any of the two closest neighbours in the -x-direction is
+              ! excised we have to excise this point
+              if ( ( eh_mask(i-1,j,k) .eq. -1 ) .or. &
+                   ( eh_mask(i-2,j,k) .eq. -1 ) ) then
+                tm_mask(i,j,k) = -1
+                mask_modified = .true.
+                print*, i, j, k, 'x+'
+                print*, eh_mask(i:i-2,j,k)
+              end if
+            end if
+
+            ! If we are on an excision boundary in the -y-direction...
+            if ( btest(eh_mask(i,j,k),2) ) then
+              ! If any of the two closest neighbours in the +y-direction is
+              ! excised we have to excise this point
+              if ( ( eh_mask(i,j+1,k) .eq. -1 ) .or. &
+                   ( eh_mask(i,j+2,k) .eq. -1 ) ) then
+                tm_mask(i,j,k) = -1
+                mask_modified = .true.
+                print*, i, j, k, 'y-'
+                print*, eh_mask(i,j:j+2,k)
+              end if
+            end if
+
+            ! If we are on an excision boundary in the +y-direction...
+            if ( btest(eh_mask(i,j,k),3) ) then
+              ! If any of the two closest neighbours in the -y-direction is
+              ! excised we have to excise this point
+              if ( ( eh_mask(i,j-1,k) .eq. -1 ) .or. &
+                   ( eh_mask(i,j-2,k) .eq. -1 ) ) then
+                tm_mask(i,j,k) = -1
+                mask_modified = .true.
+                print*, i, j, k, 'y+'
+                print*, eh_mask(i,j:j-2,k)
+              end if
+            end if
+
+            ! If we are on an excision boundary in the -z-direction...
+            if ( btest(eh_mask(i,j,k),4) ) then
+              ! If any of the two closest neighbours in the +z-direction is
+              ! excised we have to excise this point
+              if ( ( eh_mask(i,j,k+1) .eq. -1 ) .or. &
+                   ( eh_mask(i,j,k+2) .eq. -1 ) ) then
+                tm_mask(i,j,k) = -1
+                mask_modified = .true.
+                print*, i, j, k, 'z-'
+                print*, eh_mask(i,j,k:k+2)
+              end if
+            end if
+
+            ! If we are on an excision boundary in the +z-direction...
+            if ( btest(eh_mask(i,j,k),5) ) then
+              ! If any of the two closest neighbours in the -z-direction is
+              ! excised we have to excise this point
+              if ( ( eh_mask(i,j,k-1) .eq. -1 ) .or. &
+                   ( eh_mask(i,j,k-2) .eq. -1 ) ) then
+                tm_mask(i,j,k) = -1
+                mask_modified = .true.
+                print*, i, j, k, 'z+'
+                print*, eh_mask(i,j,k:k-2)
+              end if
+            end if
+          end if
+
+        end do
+      end do
+    end do
+
+    if ( mask_modified ) then
+      call CCTK_WARN ( 1, "Mask modified because it was not convex" )
+    end if
+    eh_mask = tm_mask
+  end if
+
+end subroutine EHFinder_SetMask3