diff options
Diffstat (limited to 'src/EHFinder_Sources.F90')
| -rw-r--r-- | src/EHFinder_Sources.F90 | 192 |
1 files changed, 101 insertions, 91 deletions
diff --git a/src/EHFinder_Sources.F90 b/src/EHFinder_Sources.F90 index cdf2095..65e0488 100644 --- a/src/EHFinder_Sources.F90 +++ b/src/EHFinder_Sources.F90 @@ -15,7 +15,7 @@ subroutine EHFinder_Sources(CCTK_ARGUMENTS) DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_FUNCTIONS - CCTK_INT :: i, j, k + CCTK_INT :: i, j, k, l CCTK_REAL :: idx, idy, idz, mdelta CCTK_REAL :: a, b, c CCTK_REAL :: g3xx, g3xy, g3xz, g3yy, g3yz, g3zz @@ -41,132 +41,142 @@ subroutine EHFinder_Sources(CCTK_ARGUMENTS) if ( CCTK_EQUALS ( surface_direction, 'inward' ) ) ssign = -one if ( CCTK_EQUALS ( upwind_type, 'intrinsic' ) ) then - do k = kzl, kzr - do j = jyl, jyr - do i = ixl, ixr - if ( f(i,j,k) .gt. fmin_bound - three*mdelta ) then + do l = 1, eh_number_level_sets + do k = kzl, kzr + do j = jyl, jyr + do i = ixl, ixr + if ( f(i,j,k,l) .gt. fmin_bound - three*mdelta ) then #include "include/centered_second2.h" - else + else #include "include/upwind_second2.h" - end if + end if + end do end do end do end do end if if ( CCTK_EQUALS ( upwind_type, 'shift' ) ) then - do k = kzl, kzr - do j = jyl, jyr - do i = ixl, ixr + do l = 1, eh_number_level_sets + do k = kzl, kzr + do j = jyl, jyr + do i = ixl, ixr #include "include/upwind_shift_second2.h" + end do end do end do end do end if if ( CCTK_EQUALS ( upwind_type, 'characteristic' ) ) then - do k = kzl, kzr - do j = jyl, jyr - do i = ixl, ixr + do l = 1, eh_number_level_sets + do k = kzl, kzr + do j = jyl, jyr + do i = ixl, ixr #include "include/metric.h" #include "include/upwind_characteristic_second2.h" + end do end do end do end do end if if ( CCTK_EQUALS ( metric_type, 'physical' ) ) then - 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 - - tmp1 = gyy(i,j,k)*gzz(i,j,k) - gyz(i,j,k)**2 - tmp2 = gxz(i,j,k)*gyz(i,j,k) - gxy(i,j,k)*gzz(i,j,k) - tmp3 = gxy(i,j,k)*gyz(i,j,k) - gxz(i,j,k)*gyy(i,j,k) - - idetg = one / ( gxx(i,j,k)*tmp1 + gxy(i,j,k)*tmp2 + gxz(i,j,k)*tmp3 ) + do l = 1, eh_number_level_sets + do k = kzl, kzr + do j = jyl, jyr + do i = ixl, ixr + if ( eh_mask(i,j,k,l) .ge. 0 ) then + alp2 = alp(i,j,k)**2 - g3xx = tmp1 * idetg - g3xy = tmp2 * idetg - g3xz = tmp3 * idetg - g3yy = ( gxx(i,j,k)*gzz(i,j,k) - gxz(i,j,k)**2 ) * idetg - g3yz = ( gxy(i,j,k)*gxz(i,j,k) - gxx(i,j,k)*gyz(i,j,k) ) * idetg - g3zz = ( gxx(i,j,k)*gyy(i,j,k) - gxy(i,j,k)**2 ) * idetg - - tmp1 = betax(i,j,k) * dfx(i,j,k) + & - betay(i,j,k) * dfy(i,j,k) + & - betaz(i,j,k) * dfz(i,j,k) - - tmp2 = g3xx * dfx(i,j,k)**2 + & - g3yy * dfy(i,j,k)**2 + & - g3zz * dfz(i,j,k)**2 + & - two * ( g3xy * dfx(i,j,k) * dfy(i,j,k) + & - 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 - ssign * sqrt ( alp2 * tmp2 ) + tmp1 = gyy(i,j,k)*gzz(i,j,k) - gyz(i,j,k)**2 + tmp2 = gxz(i,j,k)*gyz(i,j,k) - gxy(i,j,k)*gzz(i,j,k) + tmp3 = gxy(i,j,k)*gyz(i,j,k) - gxz(i,j,k)*gyy(i,j,k) + + idetg = one / ( gxx(i,j,k)*tmp1 + gxy(i,j,k)*tmp2 + gxz(i,j,k)*tmp3 ) + + g3xx = tmp1 * idetg + g3xy = tmp2 * idetg + g3xz = tmp3 * idetg + g3yy = ( gxx(i,j,k)*gzz(i,j,k) - gxz(i,j,k)**2 ) * idetg + g3yz = ( gxy(i,j,k)*gxz(i,j,k) - gxx(i,j,k)*gyz(i,j,k) ) * idetg + g3zz = ( gxx(i,j,k)*gyy(i,j,k) - gxy(i,j,k)**2 ) * idetg + + tmp1 = betax(i,j,k) * dfx(i,j,k,l) + & + betay(i,j,k) * dfy(i,j,k,l) + & + betaz(i,j,k) * dfz(i,j,k,l) + + tmp2 = g3xx * dfx(i,j,k,l)**2 + & + g3yy * dfy(i,j,k,l)**2 + & + g3zz * dfz(i,j,k,l)**2 + & + two * ( g3xy * dfx(i,j,k,l) * dfy(i,j,k,l) + & + g3xz * dfx(i,j,k,l) * dfz(i,j,k,l) + & + g3yz * dfy(i,j,k,l) * dfz(i,j,k,l) ) + if ( tmp2 .ge. zero ) then + sf(i,j,k,l) = tmp1 - ssign * sqrt ( alp2 * tmp2 ) + else + call CCTK_WARN ( 0, "3-metric not positive definite: Stopping" ) + end if else - call CCTK_WARN ( 0, "3-metric not positive definite: Stopping" ) + sf(i,j,k,l) = zero end if - else - sf(i,j,k) = zero - end if + end do end do - end do - end do + end do + end do end if if ( CCTK_EQUALS ( metric_type, 'static conformal' ) ) then - 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 + do l = 1, eh_number_level_sets + do k = kzl, kzr + do j = jyl, jyr + do i = ixl, ixr + if ( eh_mask(i,j,k,l) .ge. 0 ) then + alp2 = alp(i,j,k)**2 - psito4 = psi(i,j,k)**4 - gxxc = gxx(i,j,k) * psito4 - gxyc = gxy(i,j,k) * psito4 - gxzc = gxz(i,j,k) * psito4 - gyyc = gyy(i,j,k) * psito4 - gyzc = gyz(i,j,k) * psito4 - gzzc = gzz(i,j,k) * psito4 + psito4 = psi(i,j,k)**4 + gxxc = gxx(i,j,k) * psito4 + gxyc = gxy(i,j,k) * psito4 + gxzc = gxz(i,j,k) * psito4 + gyyc = gyy(i,j,k) * psito4 + gyzc = gyz(i,j,k) * psito4 + gzzc = gzz(i,j,k) * psito4 - tmp1 = gyyc*gzzc - gyzc**2 - tmp2 = gxzc*gyzc - gxyc*gzzc - tmp3 = gxyc*gyzc - gxzc*gyyc + tmp1 = gyyc*gzzc - gyzc**2 + tmp2 = gxzc*gyzc - gxyc*gzzc + tmp3 = gxyc*gyzc - gxzc*gyyc - idetg = one / ( gxxc*tmp1 + gxyc*tmp2 + gxzc*tmp3 ) - - g3xx = tmp1 * idetg - g3xy = tmp2 * idetg - g3xz = tmp3 * idetg - g3yy = ( gxxc*gzzc - gxzc**2 ) * idetg - g3yz = ( gxyc*gxzc - gxxc*gyzc ) * idetg - g3zz = ( gxxc*gyyc - gxyc**2 ) * idetg - - tmp1 = betax(i,j,k) * dfx(i,j,k) + & - betay(i,j,k) * dfy(i,j,k) + & - betaz(i,j,k) * dfz(i,j,k) - - tmp2 = g3xx * dfx(i,j,k)**2 + & - g3yy * dfy(i,j,k)**2 + & - g3zz * dfz(i,j,k)**2 + & - two * ( g3xy * dfx(i,j,k) * dfy(i,j,k) + & - 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 - ssign * sqrt ( alp2 * tmp2 ) + idetg = one / ( gxxc*tmp1 + gxyc*tmp2 + gxzc*tmp3 ) + + g3xx = tmp1 * idetg + g3xy = tmp2 * idetg + g3xz = tmp3 * idetg + g3yy = ( gxxc*gzzc - gxzc**2 ) * idetg + g3yz = ( gxyc*gxzc - gxxc*gyzc ) * idetg + g3zz = ( gxxc*gyyc - gxyc**2 ) * idetg + + tmp1 = betax(i,j,k) * dfx(i,j,k,l) + & + betay(i,j,k) * dfy(i,j,k,l) + & + betaz(i,j,k) * dfz(i,j,k,l) + + tmp2 = g3xx * dfx(i,j,k,l)**2 + & + g3yy * dfy(i,j,k,l)**2 + & + g3zz * dfz(i,j,k,l)**2 + & + two * ( g3xy * dfx(i,j,k,l) * dfy(i,j,k,l) + & + g3xz * dfx(i,j,k,l) * dfz(i,j,k,l) + & + g3yz * dfy(i,j,k,l) * dfz(i,j,k,l) ) + if ( tmp2 .ge. zero ) then + sf(i,j,k,l) = tmp1 - ssign * sqrt ( alp2 * tmp2 ) + else + call CCTK_WARN ( 0, "3-metric not positive definite: Stopping" ) + end if else - call CCTK_WARN ( 0, "3-metric not positive definite: Stopping" ) + sf(i,j,k,l) = zero end if - else - sf(i,j,k) = zero - end if + end do end do - end do - end do + end do + end do end if return |