From d1afeabbb8ff8f04de1aeda143a5c33b427a526d Mon Sep 17 00:00:00 2001 From: diener Date: Tue, 3 Jun 2003 11:44:15 +0000 Subject: Simplified the implementation of the evolution equations. git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/EHFinder/trunk@114 2a26948c-0e4f-0410-aee8-f1d3e353619c --- src/EHFinder_Sources.F90 | 134 ++++++++++++++--------------------------------- 1 file changed, 40 insertions(+), 94 deletions(-) diff --git a/src/EHFinder_Sources.F90 b/src/EHFinder_Sources.F90 index c20f25e..b629783 100644 --- a/src/EHFinder_Sources.F90 +++ b/src/EHFinder_Sources.F90 @@ -18,7 +18,6 @@ subroutine EHFinder_Sources(CCTK_ARGUMENTS) CCTK_INT :: i, j, k CCTK_REAL :: idx, idy, idz, mdelta CCTK_REAL :: a, b, c - CCTK_REAL :: g4tt, g4tx, g4ty, g4tz, g4xx, g4xy, g4xz, g4yy, g4yz, g4zz CCTK_REAL :: g3xx, g3xy, g3xz, g3yy, g3yz, g3zz CCTK_REAL :: gxxc, gxyc, gxzc, gyyc, gyzc, gzzc, psito4 CCTK_REAL :: idetg, alp2, tmp1, tmp2, tmp3 @@ -79,10 +78,6 @@ subroutine EHFinder_Sources(CCTK_ARGUMENTS) 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) @@ -90,35 +85,27 @@ subroutine EHFinder_Sources(CCTK_ARGUMENTS) 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 ) + 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 - sqrt ( alp2 * tmp2 ) else - sf(i,j,k) = two*c / ( -b - sqrt(b**2 - 4*a*c) ) + call CCTK_WARN ( 0, "3-metric not positive definite: Stopping" ) end if else sf(i,j,k) = zero @@ -134,10 +121,6 @@ subroutine EHFinder_Sources(CCTK_ARGUMENTS) 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 @@ -153,65 +136,28 @@ subroutine EHFinder_Sources(CCTK_ARGUMENTS) 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 - - 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 .le. zero ) then - sf(i,j,k) = ( -b + sqrt(b**2 - four*a*c) ) / ( two*a ) + 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 - sqrt ( alp2 * tmp2 ) else - sf(i,j,k) = two*c / ( -b - sqrt(b**2 - 4*a*c) ) + call CCTK_WARN ( 0, "3-metric not positive definite: Stopping" ) end if - sf(i,j,k) = sf(i,j,k)*sign(one,alp(i,j,k)) - -! if ( i .eq. 41 .and. j .eq. 20 .and. k .eq. 25 ) then -! print*,i,j,k -! print*,f(i,j,k),eh_mask(i,j,k) -! print*,f(i-1,j,k),f(i+1,j,k),eh_mask(i-1,j,k),eh_mask(i+1,j,k) -! print*,f(i,j-1,k),f(i,j+1,k),eh_mask(i,j-1,k),eh_mask(i,j+1,k) -! print*,f(i,j,k-1),f(i,j,k+1),eh_mask(i,j,k-1),eh_mask(i,j,k+1) -! print*,dfx(i,j,k),dfy(i,j,k),dfz(i,j,k) -! print* -! print*,x(i,j,k),y(i,j,k),z(i,j,k) -! print*,alp(i,j,k) -! print*,betax(i,j,k),betay(i,j,k),betaz(i,j,k) -! print*,gxx(i,j,k),gxy(i,j,k),gxz(i,j,k) -! print*,gyy(i,j,k),gyz(i,j,k),gzz(i,j,k) -! print*,psi(i,j,k) -! print* -! print*,tmp1,tmp2,tmp3 -! print*,gxxc*tmp1 + gxyc*tmp2 + gxzc*tmp3 -! print*,idetg -! print* -! print*,g4xx,g4xy,g4xz -! print*,g4yy,g4yz,g4zz -! print* -! print*,a,b,c -! print*,sf(i,j,k) -! print* -! print* -! end if else sf(i,j,k) = zero end if -- cgit v1.2.3