Simplified the implementation of the evolution equations.

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/EHFinder/trunk@114 2a26948c-0e4f-0410-aee8-f1d3e353619c

1 files 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