diff options
| -rw-r--r-- | interface.ccl | 2 | ||||
| -rw-r--r-- | param.ccl | 4 | ||||
| -rw-r--r-- | schedule.ccl | 5 | ||||
| -rw-r--r-- | src/EHFinder_Sources.F90 | 166 | ||||
| -rw-r--r-- | src/make.code.defn | 1 |
5 files changed, 121 insertions, 57 deletions
diff --git a/interface.ccl b/interface.ccl index 796eff6..d543c91 100644 --- a/interface.ccl +++ b/interface.ccl @@ -2,7 +2,7 @@ # $Header$ implements: ehfinder -inherits: grid admbase coordgauge boundary +inherits: grid admbase coordgauge staticconformal boundary USES INCLUDE: Boundary.h USES INCLUDE: MoL.h @@ -93,3 +93,7 @@ shares: grid USES KEYWORD domain USES KEYWORD quadrant_direction USES KEYWORD bitant_plane + +shares: admbase + +USES KEYWORD metric_type diff --git a/schedule.ccl b/schedule.ccl index e728618..4465325 100644 --- a/schedule.ccl +++ b/schedule.ccl @@ -8,6 +8,11 @@ STORAGE: eh_mask_all STORAGE: rep_mask STORAGE: ftmp +schedule EHFinder_ParamCheck at CCTK_PARAMCHECK +{ + LANG: Fortran +} "Check parameters" + schedule EHFinder_Init at CCTK_POSTINITIAL { LANG: Fortran diff --git a/src/EHFinder_Sources.F90 b/src/EHFinder_Sources.F90 index 1ccfee2..6101d54 100644 --- a/src/EHFinder_Sources.F90 +++ b/src/EHFinder_Sources.F90 @@ -19,6 +19,7 @@ subroutine EHFinder_Sources(CCTK_ARGUMENTS) CCTK_REAL :: idx, idy, idz CCTK_REAL :: a, b, c CCTK_REAL :: g4tt, g4tx, g4ty, g4tz, g4xx, g4xy, g4xz, g4yy, g4yz, g4zz + CCTK_REAL :: gxxc, gxyc, gxzc, gyyc, gyzc, gzzc, psito4 CCTK_REAL :: idetg, alp2, tmp1, tmp2, tmp3 CCTK_REAL :: ratio CCTK_REAL, dimension(3) :: maxpos @@ -250,71 +251,124 @@ subroutine EHFinder_Sources(CCTK_ARGUMENTS) ! do j = 1+iy1, ny-iy2 ! do i = 1+ix1, nx-ix2 - 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 - g4tt = -one - g4tx = betax(i,j,k) - g4ty = betay(i,j,k) - g4tz = betaz(i,j,k) + 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 + g4tt = -one + g4tx = betax(i,j,k) + g4ty = betay(i,j,k) + g4tz = betaz(i,j,k) - 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) + 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 ) + idetg = one / ( gxx(i,j,k)*tmp1 + gxy(i,j,k)*tmp2 + gxz(i,j,k)*tmp3 ) + + g4xx = tmp1 * idetg + g4xy = tmp2 * idetg + g4xz = tmp3 * idetg + g4yy = ( gxx(i,j,k)*gzz(i,j,k) - gxz(i,j,k)**2 ) * idetg + g4yz = ( gxy(i,j,k)*gxz(i,j,k) - gxx(i,j,k)*gyz(i,j,k) ) * idetg + g4zz = ( gxx(i,j,k)*gyy(i,j,k) - gxy(i,j,k)**2 ) * idetg + + g4xx = alp2*g4xx - betax(i,j,k)**2 + g4xy = alp2*g4xy - betax(i,j,k)*betay(i,j,k) + g4xz = alp2*g4xz - betax(i,j,k)*betaz(i,j,k) + g4yy = alp2*g4yy - betay(i,j,k)**2 + g4yz = alp2*g4yz - betay(i,j,k)*betaz(i,j,k) + g4zz = alp2*g4zz - betaz(i,j,k)**2 + + a = g4tt + b = two * ( g4tx * dfx(i,j,k) + & + g4ty * dfy(i,j,k) + & + g4tz * dfz(i,j,k) ) + c = g4xx * dfx(i,j,k)**2 + & + g4yy * dfy(i,j,k)**2 + & + g4zz * dfz(i,j,k)**2 + & + two * ( g4xy * dfx(i,j,k) * dfy(i,j,k) + & + g4xz * dfx(i,j,k) * dfz(i,j,k) + & + g4yz * dfy(i,j,k) * dfz(i,j,k) ) + + if ( b .lt. zero ) then + sf(i,j,k) = ( -b + sqrt(b**2 - four*a*c) ) / ( two*a ) + else + sf(i,j,k) = two*c / ( -b - sqrt(b**2 - 4*a*c) ) + end if + else + sf(i,j,k) = zero + end if + end do + end do + end do + end if - g4xx = tmp1 * idetg - g4xy = tmp2 * idetg - g4xz = tmp3 * idetg - g4yy = ( gxx(i,j,k)*gzz(i,j,k) - gxz(i,j,k)**2 ) * idetg - g4yz = ( gxy(i,j,k)*gxz(i,j,k) - gxx(i,j,k)*gyz(i,j,k) ) * idetg - g4zz = ( gxx(i,j,k)*gyy(i,j,k) - gxy(i,j,k)**2 ) * idetg + 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 + g4tt = -one + g4tx = betax(i,j,k) + g4ty = betay(i,j,k) + g4tz = betaz(i,j,k) + + 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 + + idetg = one / ( gxxc*tmp1 + gxyc*tmp2 + gxzc*tmp3 ) + + g4xx = tmp1 * idetg + g4xy = tmp2 * idetg + g4xz = tmp3 * idetg + g4yy = ( gxxc*gzzc - gxzc**2 ) * idetg + g4yz = ( gxyc*gxzc - gxxc*gyzc ) * idetg + g4zz = ( gxxc*gyyc - gxyc**2 ) * idetg - g4xx = alp2*g4xx - betax(i,j,k)**2 - g4xy = alp2*g4xy - betax(i,j,k)*betay(i,j,k) - g4xz = alp2*g4xz - betax(i,j,k)*betaz(i,j,k) - g4yy = alp2*g4yy - betay(i,j,k)**2 - g4yz = alp2*g4yz - betay(i,j,k)*betaz(i,j,k) - g4zz = alp2*g4zz - betaz(i,j,k)**2 + g4xx = alp2*g4xx - betax(i,j,k)**2 + g4xy = alp2*g4xy - betax(i,j,k)*betay(i,j,k) + g4xz = alp2*g4xz - betax(i,j,k)*betaz(i,j,k) + g4yy = alp2*g4yy - betay(i,j,k)**2 + g4yz = alp2*g4yz - betay(i,j,k)*betaz(i,j,k) + g4zz = alp2*g4zz - betaz(i,j,k)**2 - a = g4tt - b = two * ( g4tx * dfx(i,j,k) + & - g4ty * dfy(i,j,k) + & - g4tz * dfz(i,j,k) ) - c = g4xx * dfx(i,j,k)**2 + & - g4yy * dfy(i,j,k)**2 + & - g4zz * dfz(i,j,k)**2 + & - two * ( g4xy * dfx(i,j,k) * dfy(i,j,k) + & - g4xz * dfx(i,j,k) * dfz(i,j,k) + & - g4yz * dfy(i,j,k) * dfz(i,j,k) ) + a = g4tt + b = two * ( g4tx * dfx(i,j,k) + & + g4ty * dfy(i,j,k) + & + g4tz * dfz(i,j,k) ) + c = g4xx * dfx(i,j,k)**2 + & + g4yy * dfy(i,j,k)**2 + & + g4zz * dfz(i,j,k)**2 + & + two * ( g4xy * dfx(i,j,k) * dfy(i,j,k) + & + g4xz * dfx(i,j,k) * dfz(i,j,k) + & + g4yz * dfy(i,j,k) * dfz(i,j,k) ) - if ( b .lt. zero ) then - sf(i,j,k) = ( -b + sqrt(b**2 - four*a*c) ) / ( two*a ) + if ( b .lt. zero ) then + sf(i,j,k) = ( -b + sqrt(b**2 - four*a*c) ) / ( two*a ) + else + sf(i,j,k) = two*c / ( -b - sqrt(b**2 - 4*a*c) ) + end if + sf(i,j,k) = sf(i,j,k)*sign(one,alp(i,j,k)) else - sf(i,j,k) = two*c / ( -b - sqrt(b**2 - 4*a*c) ) + sf(i,j,k) = zero end if - else - sf(i,j,k) = zero - end if -! if ( CCTK_MyProc(cctkGH) .eq. 0 ) then -! if ( abs(b) .gt. 0.00000001 ) then -! print*,i,j,k -! print*,dfx(i,j,k),dfy(i,j,k),dfz(i,j,k) -! print*,g4tx,g4ty,g4tz -! print*,idetg -! print*,g4xx,g4yy,g4zz -! print*,g4xy,g4xz,g4yz -! print*,a,b,c -! print*,sf(i,j,k) -! pause -! endif -! end if + end do end do - end do - end do + end do + end if ! print*,maxval(sf) ! print* diff --git a/src/make.code.defn b/src/make.code.defn index 47d924a..c2cafaa 100644 --- a/src/make.code.defn +++ b/src/make.code.defn @@ -3,6 +3,7 @@ # Source files in this directory SRCS = EHFinder_mod.F90 \ + EHFinder_ParamCheck.F90 \ EHFinder_Init.F90 \ MoL_Init.F90 \ EHFinder_Sources.F90 \ |