aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--param.ccl6
-rw-r--r--src/EHFinder_Generator_Sources.F9017
-rw-r--r--src/EHFinder_Generator_Sources2.F9017
-rw-r--r--src/EHFinder_Sources.F909
-rw-r--r--src/include/upwind_characteristic_second2.h9
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
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-03 14:49:29 +0000