Added a parameter to schoose inward going surfaces instead of outward

going surfaces. If tracking inward going surfaces it is not an event horizon
finder anymore, but this feature might be useful for other things.


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

5 files changed, 38 insertions, 20 deletions

diff --git a/param.ccl b/param.ccl
index e20eeeb..8f62563 100644
--- a/param.ccl
+++ b/param.ccl
@@ -128,6 +128,12 @@ KEYWORD upwind_type "Type of upwinding used in evolving the ehfinder equations"
   "characteristic" :: "Use characteristic information"
 } "intrinsic"
 
+KEYWORD surface_direction "Should we track outward or inward moving surfaces"
+{
+  "outward" :: "Track outward moving surfaces. Use for event horizon finding."
+  "inward"  :: "Track inward moving surfaces."
+} "outward"
+
 BOOLEAN reparam_undo "Should re-parametrization be undone at pinch-off"
 {
 } "no"
diff --git a/src/EHFinder_Generator_Sources.F90 b/src/EHFinder_Generator_Sources.F90
index 95cb22e..a89a548 100644
--- a/src/EHFinder_Generator_Sources.F90
+++ b/src/EHFinder_Generator_Sources.F90
@@ -33,9 +33,12 @@ subroutine EHFinder_Generator_Sources(CCTK_ARGUMENTS)
                                                       0, 0, 0, 0, 0, &
                                                       1, 2, 3, 0 /)
   CCTK_INT, dimension(14) :: out_types
-  CCTK_REAL :: alp2, psi4, dfux, dfuy, dfuz, factor
+  CCTK_REAL :: alp2, psi4, dfux, dfuy, dfuz, factor, ssign
   CCTK_REAL :: idetg, guxx, guxy, guxz, guyy, guyz, guzz
 
+  if ( CCTK_EQUALS ( surface_direction, 'outward' ) ) ssign = one
+  if ( CCTK_EQUALS ( surface_direction, 'inward' ) ) ssign = -one
+
   out_types = CCTK_VARIABLE_REAL
 
   ! Convert the generator_interpolator string parameter to a Fortran string.
@@ -172,9 +175,9 @@ subroutine EHFinder_Generator_Sources(CCTK_ARGUMENTS)
                                       dfuz * dfzg(i) ) )
 
       ! Finally obtain dx^i/dt.
-      dxg(i) = - betaxg(i) + factor * dfux
-      dyg(i) = - betayg(i) + factor * dfuy
-      dzg(i) = - betazg(i) + factor * dfuz
+      dxg(i) = - betaxg(i) + ssign * factor * dfux
+      dyg(i) = - betayg(i) + ssign * factor * dfuy
+      dzg(i) = - betazg(i) + ssign * factor * dfuz
     end do
   else if ( CCTK_EQUALS ( metric_type, 'static conformal' ) ) then
     call Util_TableSetIntArray ( status, table_handle, 14, &
@@ -230,9 +233,9 @@ subroutine EHFinder_Generator_Sources(CCTK_ARGUMENTS)
       factor = alp2 / sqrt ( alp2 * ( dfux * dfxg(i) + &
                                       dfuy * dfyg(i) + &
                                       dfuz * dfzg(i) ) )
-      dxg(i) = - betaxg(i) + factor * dfux
-      dyg(i) = - betayg(i) + factor * dfuy
-      dzg(i) = - betazg(i) + factor * dfuz
+      dxg(i) = - betaxg(i) + ssign * factor * dfux
+      dyg(i) = - betayg(i) + ssign * factor * dfuy
+      dzg(i) = - betazg(i) + ssign * factor * dfuz
     end do
   end if
   return
diff --git a/src/EHFinder_Generator_Sources2.F90 b/src/EHFinder_Generator_Sources2.F90
index aab47e1..b50f9f1 100644
--- a/src/EHFinder_Generator_Sources2.F90
+++ b/src/EHFinder_Generator_Sources2.F90
@@ -29,9 +29,12 @@ subroutine EHFinder_Generator_Sources2(CCTK_ARGUMENTS)
   CCTK_INT, dimension(3), parameter :: op_indices = (/ 0, 1, 2 /), &
                                         op_codes = (/ 0, 0, 0 /)
   CCTK_INT, dimension(3) :: out_types
-  CCTK_REAL :: alp2, psi4, dfux, dfuy, dfuz, factor
+  CCTK_REAL :: alp2, psi4, dfux, dfuy, dfuz, factor, ssign
   CCTK_REAL :: idetg, guxx, guxy, guxz, guyy, guyz, guzz
 
+  if ( CCTK_EQUALS ( surface_direction, 'outward' ) ) ssign = one
+  if ( CCTK_EQUALS ( surface_direction, 'inward' ) ) ssign = -one
+
   out_types = CCTK_VARIABLE_REAL
 
   ! Convert the generator_interpolator string parameter to a Fortran string.
@@ -140,9 +143,9 @@ subroutine EHFinder_Generator_Sources2(CCTK_ARGUMENTS)
                                             dfuz * dfz(i,j,k) ) )
 
             ! Finally obtain dx^i/dt.
-            xgf(i,j,k) = - betax(i,j,k) + factor * dfux
-            ygf(i,j,k) = - betay(i,j,k) + factor * dfuy
-            zgf(i,j,k) = - betaz(i,j,k) + factor * dfuz
+            xgf(i,j,k) = - betax(i,j,k) + ssign * factor * dfux
+            ygf(i,j,k) = - betay(i,j,k) + ssign * factor * dfuy
+            zgf(i,j,k) = - betaz(i,j,k) + ssign * factor * dfuz
           else
             xgf(i,j,k) = zero
             ygf(i,j,k) = zero
@@ -185,9 +188,9 @@ subroutine EHFinder_Generator_Sources2(CCTK_ARGUMENTS)
       factor = alp2 / sqrt ( alp2 * ( dfux * dfxg(i) + &
                                       dfuy * dfyg(i) + &
                                       dfuz * dfzg(i) ) )
-      dxg(i) = - betaxg(i) + factor * dfux
-      dyg(i) = - betayg(i) + factor * dfuy
-      dzg(i) = - betazg(i) + factor * dfuz
+      dxg(i) = - betaxg(i) + ssign * factor * dfux
+      dyg(i) = - betayg(i) + ssign * factor * dfuy
+      dzg(i) = - betazg(i) + ssign * factor * dfuz
     end do
   end if
 
diff --git a/src/EHFinder_Sources.F90 b/src/EHFinder_Sources.F90
index b629783..cdf2095 100644
--- a/src/EHFinder_Sources.F90
+++ b/src/EHFinder_Sources.F90
@@ -21,7 +21,7 @@ subroutine EHFinder_Sources(CCTK_ARGUMENTS)
   CCTK_REAL :: g3xx, g3xy, g3xz, g3yy, g3yz, g3zz
   CCTK_REAL :: gxxc, gxyc, gxzc, gyyc, gyzc, gzzc, psito4
   CCTK_REAL :: idetg, alp2, tmp1, tmp2, tmp3
-  CCTK_REAL :: ratio, cfactor
+  CCTK_REAL :: ratio, cfactor, ssign
   CCTK_REAL, dimension(3) :: maxpos, cdx, dfup
   CCTK_REAL :: al, ar, bl, br, cl, cr
   CCTK_REAL :: alminus, alplus, blminus, blplus, clminus, clplus
@@ -37,6 +37,9 @@ subroutine EHFinder_Sources(CCTK_ARGUMENTS)
 
   mdelta = maxval ( cctk_delta_space )
 
+  if ( CCTK_EQUALS ( surface_direction, 'outward' ) ) ssign = one
+  if ( CCTK_EQUALS ( surface_direction, 'inward' ) ) ssign = -one
+
   if ( CCTK_EQUALS ( upwind_type, 'intrinsic' ) ) then
     do k = kzl, kzr
       do j = jyl, jyr
@@ -103,7 +106,7 @@ subroutine EHFinder_Sources(CCTK_ARGUMENTS)
                            g3xz * dfx(i,j,k) * dfz(i,j,k) + &
                            g3yz * dfy(i,j,k) * dfz(i,j,k) )
             if ( tmp2 .ge. zero ) then
-              sf(i,j,k) = tmp1 - sqrt ( alp2 * tmp2 )
+              sf(i,j,k) = tmp1 - ssign * sqrt ( alp2 * tmp2 )
             else
               call CCTK_WARN ( 0, "3-metric not positive definite: Stopping" )
             end if
@@ -154,7 +157,7 @@ subroutine EHFinder_Sources(CCTK_ARGUMENTS)
                            g3xz * dfx(i,j,k) * dfz(i,j,k) + &
                            g3yz * dfy(i,j,k) * dfz(i,j,k) )
             if ( tmp2 .ge. zero ) then
-              sf(i,j,k) = tmp1 - sqrt ( alp2 * tmp2 )
+              sf(i,j,k) = tmp1 - ssign * sqrt ( alp2 * tmp2 )
             else
               call CCTK_WARN ( 0, "3-metric not positive definite: Stopping" )
             end if
diff --git a/src/include/upwind_characteristic_second2.h b/src/include/upwind_characteristic_second2.h
index 25c317d..5dacaa4 100644
--- a/src/include/upwind_characteristic_second2.h
+++ b/src/include/upwind_characteristic_second2.h
@@ -48,9 +48,12 @@ if ( eh_mask(i,j,k) .ge. 0 ) then
                                   dfup(2) * dfy(i,j,k) + &
                                   dfup(3) * dfz(i,j,k) ) )
 
-  cdx(1) = ( -betax(i,j,k) + alp2 * dfup(1) * cfactor ) * cctk_delta_time
-  cdx(2) = ( -betay(i,j,k) + alp2 * dfup(2) * cfactor ) * cctk_delta_time
-  cdx(3) = ( -betaz(i,j,k) + alp2 * dfup(3) * cfactor ) * cctk_delta_time
+  cdx(1) = ( -betax(i,j,k) + ssign * alp2 * dfup(1) * cfactor ) * &
+                                                      cctk_delta_time
+  cdx(2) = ( -betay(i,j,k) + ssign * alp2 * dfup(2) * cfactor ) * &
+                                                      cctk_delta_time
+  cdx(3) = ( -betaz(i,j,k) + ssign * alp2 * dfup(3) * cfactor ) * &
+                                                      cctk_delta_time
 
   if ( ( .not. btest(eh_mask(i,j,k),0) ) .and. &
        ( .not. btest(eh_mask(i,j,k),1) ) ) then