From 24fccc7f8f07b8891cb883165fad85e3ee6d18df Mon Sep 17 00:00:00 2001 From: diener Date: Mon, 1 Sep 2003 12:29:47 +0000 Subject: Remove this experimental routine, that turned out to be a dead end. git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/EHFinder/trunk@131 2a26948c-0e4f-0410-aee8-f1d3e353619c --- src/EHFinder_Integrate.F90 | 340 --------------------------------------------- 1 file changed, 340 deletions(-) delete mode 100644 src/EHFinder_Integrate.F90 diff --git a/src/EHFinder_Integrate.F90 b/src/EHFinder_Integrate.F90 deleted file mode 100644 index 04a2a46..0000000 --- a/src/EHFinder_Integrate.F90 +++ /dev/null @@ -1,340 +0,0 @@ -! Integration over the surface -! $Header$ - -#include "cctk.h" -#include "cctk_Parameters.h" -#include "cctk_Arguments.h" - -subroutine EHFinder_Integrate(CCTK_ARGUMENTS) - - use EHFinder_mod - - implicit none - - DECLARE_CCTK_PARAMETERS - DECLARE_CCTK_ARGUMENTS - DECLARE_CCTK_FUNCTIONS - - CCTK_INT :: i, j, k - CCTK_REAL :: idx, idy, idz - CCTK_REAL :: al, ar, bl, br, cl, cr - CCTK_REAL :: alminus, alplus, blminus, blplus, clminus, clplus - CCTK_REAL :: arminus, arplus, brminus, brplus, crminus, crplus - CCTK_REAL :: dv, local_area, area, eps, eps_inv - CCTK_REAL, dimension(3,3) :: ginv - CCTK_REAL :: g2xx, g2xy, g2xz, g2yy, g2yz, g2zz - CCTK_REAL :: g22xx, g22xy, g22yy - CCTK_REAL, dimension(3) :: norm - CCTK_REAL :: dfsquare, dfsquare_i - CCTK_REAL, dimension(3) :: vec1, vec2 - CCTK_REAL :: tmp1, tmp2, tmp3 - CCTK_REAL :: detg, idetg, sqrtdg2 - CCTK_REAL :: a, b, c, r1, r2, r3, q - - interface - CCTK_REAL function dirac ( x, eps, eps_inv ) - CCTK_REAL, intent(in) :: x, eps, eps_inv - end function dirac - end interface - -# include "include/physical_part.h" - - idx = half / cctk_delta_space(1) - idy = half / cctk_delta_space(2) - idz = half / cctk_delta_space(3) - - dv = cctk_delta_space(1)*cctk_delta_space(2)*cctk_delta_space(3) - eps = minval(cctk_delta_space) - eps_inv = one / eps - - local_area = zero - - if ( CCTK_EQUALS ( metric_type, 'physical' ) ) then - do k = kzl, kzr - do j = jyl, jyr - do i = ixl, ixr - if ( abs(f(i,j,k)) .lt. eps ) then - -#include "include/upwind_second2.h" - -! Calculate the inverse of the determinant of the spatial 3-metric - - 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) - - detg = gxx(i,j,k)*tmp1 + gxy(i,j,k)*tmp2 + gxz(i,j,k)*tmp3 - idetg = one / detg - -! Invert the spatial 3-metric - - ginv(1,1) = tmp1 * idetg - ginv(1,2) = tmp2 * idetg - ginv(1,3) = tmp3 * idetg - ginv(2,1) = ginv(1,2) - ginv(2,2) = ( gxx(i,j,k)*gzz(i,j,k) - gxz(i,j,k)**2 ) * idetg - ginv(2,3) = ( gxy(i,j,k)*gxz(i,j,k) - gxx(i,j,k)*gyz(i,j,k) )*idetg - ginv(3,1) = ginv(1,3) - ginv(3,2) = ginv(2,3) - ginv(3,3) = ( gxx(i,j,k)*gyy(i,j,k) - gxy(i,j,k)**2 ) * idetg - - dfsquare = sqrt ( ginv(1,1) * dfx(i,j,k)**2 + & - two * ginv(1,2) * dfx(i,j,k) * dfy(i,j,k) + & - two * ginv(1,3) * dfx(i,j,k) * dfz(i,j,k) + & - ginv(2,2) * dfy(i,j,k)**2 + & - two * ginv(2,3) * dfy(i,j,k)*dfz(i,j,k) + & - ginv(3,3) * dfz(i,j,k)**2 ) -! print*,dfsquare - - dfsquare_i = one / dfsquare -! print*,dfsquare_i -! dfsquare_i = one / sqrt ( dfx(i,j,k)**2 + & -! dfy(i,j,k)**2 + & -! dfz(i,j,k)**2 ) - - norm(1) = dfx(i,j,k) * dfsquare_i - norm(2) = dfy(i,j,k) * dfsquare_i - norm(3) = dfz(i,j,k) * dfsquare_i -! print*,norm - -! norm(1) = dfx(i,j,k) -! norm(2) = dfy(i,j,k) -! norm(3) = dfz(i,j,k) - - call find_orthogonal ( norm, vec1, vec2, ginv, idetg ) -! print*,norm(1),norm(2),norm(3) -! print*,vec1 -! print*,vec2 -! print*,norm(1)*vec1(1)+norm(2)*vec1(2)+norm(3)*vec1(3) -! print*,norm(1)*vec2(1)+norm(2)*vec2(2)+norm(3)*vec2(3) -! print*,gxx(i,j,k)*vec1(1)*vec2(1)+gxy(i,j,k)*vec1(1)*vec2(2)+& -! gxz(i,j,k)*vec1(1)*vec2(3)+gxy(i,j,k)*vec1(2)*vec2(1)+& -! gyy(i,j,k)*vec1(2)*vec2(2)+gyz(i,j,k)*vec1(2)*vec2(3)+& -! gxz(i,j,k)*vec1(3)*vec2(1)+gyz(i,j,k)*vec1(3)*vec2(2)+& -! gzz(i,j,k)*vec1(3)*vec2(3) -! print*,gxx(i,j,k)*vec1(1)*vec1(1)+gxy(i,j,k)*vec1(1)*vec1(2)+& -! gxz(i,j,k)*vec1(1)*vec1(3)+gxy(i,j,k)*vec1(2)*vec1(1)+& -! gyy(i,j,k)*vec1(2)*vec1(2)+gyz(i,j,k)*vec1(2)*vec1(3)+& -! gxz(i,j,k)*vec1(3)*vec1(1)+gyz(i,j,k)*vec1(3)*vec1(2)+& -! gzz(i,j,k)*vec1(3)*vec1(3) -! print*,gxx(i,j,k)*vec2(1)*vec2(1)+gxy(i,j,k)*vec2(1)*vec2(2)+& -! gxz(i,j,k)*vec2(1)*vec2(3)+gxy(i,j,k)*vec2(2)*vec2(1)+& -! gyy(i,j,k)*vec2(2)*vec2(2)+gyz(i,j,k)*vec2(2)*vec2(3)+& -! gxz(i,j,k)*vec2(3)*vec2(1)+gyz(i,j,k)*vec2(3)*vec2(2)+& -! gzz(i,j,k)*vec2(3)*vec2(3) - -! print* -! print*,gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),& -! gyy(i,j,k),gyz(i,j,k),gzz(i,j,k) -! print* -! pause - g2xx = gxx(i,j,k) - norm(1)**2 - g2xy = gxy(i,j,k) - norm(1) * norm(2) - g2xz = gxz(i,j,k) - norm(1) * norm(3) - g2yy = gyy(i,j,k) - norm(2)**2 - g2yz = gyz(i,j,k) - norm(2) * norm(3) - g2zz = gzz(i,j,k) - norm(3)**2 -! print*,g2xx,g2xy,g2xz,g2yy,g2yz,g2zz - -! a = - ( g2xx + g2yy + g2zz ) - b = g2xx * g2yy + g2xx * g2zz + g2yy * g2zz - & - g2xy**2 - g2xz**2 - g2yz**2 -! c = g2xy**2 * g2zz + g2xz**2 * g2yy + g2yz**2 * g2xx - & -! g2xx * g2yy * g2zz - two * g2xy * g2xz * g2yz - -! q = - half * ( a + sign ( one, a ) * sqrt ( a**2 - 4 * b ) ) -! print*,'factors : ',a,b,c,q,b/q,sqrt(b) -! print* -! sqrtdg2 = sqrt ( sqrt (b) ) - sqrtdg2 = sqrt (b) -! call cubic ( a, b, c, r1, r2, r3 ) -! print*,'roots : ',r1,r2,r3 -! print* -! pause -! print*,g2xx,g2xy,g2xz,g2yy,g2yz,g2zz -! print* - -! g22xx = vec1(1)**2*ginv(1,1) + two*vec1(1)*vec1(2)*ginv(1,2) +& -! two*vec1(1)*vec1(3)*ginv(1,3) + vec1(2)**2*ginv(2,2) +& -! two*vec1(2)*vec1(3)*ginv(2,3) + vec1(3)**2 *ginv(3,3) -! g22xy = vec1(1)*vec2(1)*ginv(1,1) + vec1(1)*vec2(2)*ginv(1,2) +& -! vec1(1)*vec2(3)*ginv(1,3) + vec1(2)*vec2(1)*ginv(2,1) +& -! vec1(2)*vec2(2)*ginv(2,2) + vec1(2)*vec2(3)*ginv(2,3) +& -! vec1(3)*vec2(1)*ginv(3,1) + vec1(3)*vec2(2)*ginv(3,2) +& -! vec1(3)*vec2(3)*ginv(3,3) -! g22yy = vec2(1)**2*ginv(1,1) + two*vec2(1)*vec2(2)*ginv(1,2) +& -! two*vec2(1)*vec2(3)*ginv(1,3) + vec2(2)**2*ginv(2,2) +& -! two*vec2(2)*vec2(3)*ginv(2,3) + vec2(3)**2*ginv(3,3) - -! g22xx = vec1(1)**2*gxx(i,j,k) + two*vec1(1)*vec1(2)*gxy(i,j,k) +& -! two*vec1(1)*vec1(3)*gxz(i,j,k) + vec1(2)**2*gyy(i,j,k) +& -! two*vec1(2)*vec1(3)*gyz(i,j,k) + vec1(3)**2 *gzz(i,j,k) -! g22xy = vec1(1)*vec2(1)*gxx(i,j,k) + vec1(1)*vec2(2)*gxy(i,j,k) +& -! vec1(1)*vec2(3)*gxz(i,j,k) + vec1(2)*vec2(1)*gxy(i,j,k) +& -! vec1(2)*vec2(2)*gyy(i,j,k) + vec1(2)*vec2(3)*gyz(i,j,k) +& -! vec1(3)*vec2(1)*gxz(i,j,k) + vec1(3)*vec2(2)*gyz(i,j,k) +& -! vec1(3)*vec2(3)*gzz(i,j,k) -! g22yy = vec2(1)**2*gxx(i,j,k) + two*vec2(1)*vec2(2)*gxy(i,j,k) +& -! two*vec2(1)*vec2(3)*gxz(i,j,k) + vec2(2)**2*gyy(i,j,k) +& -! two*vec2(2)*vec2(3)*gyz(i,j,k) + vec2(3)**2*gzz(i,j,k) - - g22xx = vec1(1)**2*g2xx + two*vec1(1)*vec1(2)*g2xy +& - two*vec1(1)*vec1(3)*g2xz + vec1(2)**2*g2yy +& - two*vec1(2)*vec1(3)*g2yz + vec1(3)**2*g2zz - g22xy = vec1(1)*vec2(1)*g2xx + vec1(1)*vec2(2)*g2xy +& - vec1(1)*vec2(3)*g2xz + vec1(2)*vec2(1)*g2xy +& - vec1(2)*vec2(2)*g2yy + vec1(2)*vec2(3)*g2yz +& - vec1(3)*vec2(1)*g2xz + vec1(3)*vec2(2)*g2yz +& - vec1(3)*vec2(3)*g2zz - g22yy = vec2(1)**2*g2xx + two*vec2(1)*vec2(2)*g2xy +& - two*vec2(1)*vec2(3)*g2xz + vec2(2)**2*g2yy +& - two*vec2(2)*vec2(3)*g2yz + vec2(3)**2*g2zz -! print* -! print*,g22xx,g22xy,g22yy -! print*,g22xx*g22yy-g22xy**2 - - tmp1 = g2yy*g2zz - g2yz**2 - tmp2 = g2xz*g2yz - g2xy*g2zz - tmp3 = g2xy*g2yz - g2xz*g2yy -! print*,tmp1,tmp2,tmp3 - -! print*,tmp1,tmp2,tmp3 -! -! print*, g2xx*tmp1 + g2xy*tmp2 + g2xz*tmp3 -! sqrtdg2 = sqrt ( g2xx*tmp1 + g2xy*tmp2 + g2xz*tmp3 ) -! sqrtdg2 = one / sqrt ( g22xx*g22yy-g22xy**2 ) -! sqrtdg2 = sqrt ( g22xx*g22yy-g22xy**2 ) -! sqrtdg2 = sqrt ( g22xx*g22yy-g22xy**2 ) - -! print*,sqrtdg2,g2xx*tmp1 + g2xy*tmp2 + g2xz*tmp3 -! pause - - local_area = local_area + sqrtdg2 * dirac(f(i,j,k),eps,eps_inv) * & - dfsquare * sqrt(detg) -! sqrt ( dfx(i,j,k)**2 + dfy(i,j,k)**2 + dfz(i,j,k)**2 ) - end if - end do - end do - end do - end if - local_area = dv * local_area -! print*,'local_area = ',local_area,sqrtdg2 - - call CCTK_ReduceLocScalar ( ierr, cctkGH, -1, sum_handle, & - local_area, area, CCTK_VARIABLE_REAL ) - - eh_area = area - print*,'Area = ',eh_area - - return -end subroutine EHFinder_Integrate - - -subroutine find_orthogonal ( n, v1, v2, ginv, idetg ) - - use EHFinder_Constants - - CCTK_REAL, dimension(3) :: n, v1, v2 - CCTK_REAL, dimension(3) :: v1l - CCTK_REAL, dimension(3,3) :: ginv - CCTK_REAL :: idetg - CCTK_REAL, dimension(1) :: maxn - CCTK_REAL :: nv1, nv2 - CCTK_INT, dimension(1) :: maxp - CCTK_INT :: first, second, third - - maxn = maxval (n) - maxp = maxloc (n) - -! print*,'Orth' -! print*,n -! print*,maxn -! print*,maxp - third = maxp(1) - if ( third .eq. 1 ) then - first = 2 - second =3 - else if ( third .eq. 2 ) then - first = 3 - second = 1 - else - first = 1 - second = 2 - end if - -! print*,first,second,third - v1l(first) = -n(third) - v1l(second) = n(first) - v1l(third) = - ( ginv(first,first) * n(first) * v1l(first) + & - ginv(first,second) * n(first) * v1l(second) + & - ginv(second,first) * n(second) * v1l(first) + & - ginv(second,second) * n(second) * v1l(second) + & - ginv(third,first) * n(third) * v1l(first) + & - ginv(third,second) * n(third) * v1l(second) ) / & - ( ginv(first,third) * n(first) + & - ginv(second,third) * n(second) + & - ginv(third,third) * n(third) ) - - nv1 = sqrt ( ginv(1,1) * v1l(1)**2 + & - two * ginv(1,2) * v1l(1) * v1l(2) + & - two * ginv(1,3) * v1l(1) * v1l(3) + & - ginv(2,2) * v1l(2)**2 + & - two * ginv(2,3) * v1l(2) * v1l(3) + & - ginv(3,3) * v1l(3)**2 ) - - v1l = v1l / nv1 - v1(1) = ginv(1,1) * v1l(1) + ginv(1,2) * v1l(2) + ginv(1,3) * v1l(3) - v1(2) = ginv(2,1) * v1l(1) + ginv(2,2) * v1l(2) + ginv(2,3) * v1l(3) - v1(3) = ginv(3,1) * v1l(1) + ginv(3,2) * v1l(2) + ginv(3,3) * v1l(3) - -! print*, v1(1)*v1l(1)+v1(2)*v1l(2)+v1(3)*v1l(3) -! print* - - v2(1) = n(2) * v1l(3) - n(3) * v1l(2) - v2(2) = n(3) * v1l(1) - n(1) * v1l(3) - v2(3) = n(1) * v1l(2) - n(2) * v1l(1) - v2 = v2 * sqrt(idetg) - -! print*, v2 -! print* -! print*,v1(1)*n(1)+v1(2)*n(2)+v1(3)*n(3) -! print*,v2(1)*n(1)+v2(2)*n(2)+v2(3)*n(3) -! print*,v2(1)*v1l(1)+v2(2)*v1l(2)+v2(3)*v1l(3) -! print*,'Orth end' -! pause - return -end subroutine find_orthogonal - -CCTK_REAL function dirac ( x, eps, eps_inv ) - - use EHFinder_Constants - - CCTK_REAL, intent(in) :: x, eps, eps_inv - - if ( abs(x) .lt. eps ) then - dirac = half * ( one + cos(pi*x*eps_inv) ) * eps_inv - else - dirac = zero - end if -end function dirac - -subroutine cubic ( a, b, c, x1, x2, x3 ) - - CCTK_REAL, intent(in) :: a, b, c - CCTK_REAL, intent(out) :: x1, x2, x3 - CCTK_REAL :: q, r, theta - CCTK_REAL, parameter :: third = 1.0d0 / 3.0d0, & - pi2 = 6.2831853071795864769d0 - - q = ( a**2 - 3.0d0 * b ) / 9.0d0 - r = ( 2.0d0 * a**3 - 9.0d0 * a * b + 27.0d0 * c ) / 54.0d0 - theta = acos ( r / sqrt( q**3 ) ) - - x1 = - 2.0d0 * sqrt( q ) * cos( third * theta ) - third * a - x2 = - 2.0d0 * sqrt( q ) * cos( third * ( theta + pi2 ) ) - third * a - x2 = - 2.0d0 * sqrt( q ) * cos( third * ( theta - pi2 ) ) - third * a - - print*,a,b,c - print*,q,r,q**3-r**2 -! pause -end subroutine cubic -- cgit v1.2.3