From 116ecf1f8ad5bbd3ad74f9367bf7907d8c43d0ce Mon Sep 17 00:00:00 2001 From: diener Date: Wed, 2 Jul 2003 09:48:57 +0000 Subject: 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 --- param.ccl | 6 ++++++ src/EHFinder_Generator_Sources.F90 | 17 ++++++++++------- src/EHFinder_Generator_Sources2.F90 | 17 ++++++++++------- src/EHFinder_Sources.F90 | 9 ++++++--- src/include/upwind_characteristic_second2.h | 9 ++++++--- 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 -- cgit v1.2.3