/*@@ c @file AHFinder_gau.F c @date October 1998 c @author Miguel Alcubierre c @desc c Find gaussian curvature of surface. The gaussian c curvature is defined as -R/2, where R is the Ricci c scalar of the induced 2-geometry of the surface. c c As an extra "goody", this routine also integrates c the equatorial circumference of the horizon and c two polar circumferences at phi=0 and phi=pi/2. c @enddesc c@@*/ c Note that including cactus.h will also include AHFinder.h !#include "cactus.h" #include "cctk.h" #include "cctk_parameters.h" #include "cctk_arguments.h" subroutine AHFinder_gau(CCTK_FARGUMENTS,circ_eq,meri_p1,meri_p2) use AHFinder_dat implicit none DECLARE_CCTK_FARGUMENTS DECLARE_CCTK_PARAMETERS c#ifdef MPI c include 'mpif.h' c#endif logical firstgau logical firstcal(4) integer gxx_index,gyy_index,gzz_index,gxy_index,gxz_index,gyz_index integer ahfgauss_index integer handle,x_index,y_index,z_index CCTK_REAL maxval(3),minval(3) INTEGER istat CCTK_INT i,j,k,l,m,n,p CCTK_INT npoints CCTK_INT error1,error2,ierror CCTK_INT getoutpfx CCTK_REAL LEGEN CCTK_REAL theta,phi,xp,yp,zp,rp CCTK_REAL cost,sint,cosp,sinp CCTK_REAL dtheta,dphi,dtp,idtheta,idphi CCTK_REAL xmn,ymn,zmn,xmx,ymx,zmx CCTK_REAL trr,ttt,tpp,trt,trp,ttp,ft,fp CCTK_REAL circ_eq,meri_p1,meri_p2 CCTK_REAL deta,ideta CCTK_REAL nlx,nly,nlz,nux,nuy,nuz CCTK_REAL trxi,xi2 CCTK_REAL zero,half,one,two,three,four,pi CCTK_REAL aux,sina,cosa CCTK_REAL, dimension(3,3) :: ug,xi CCTK_REAL, dimension(2,2) :: ga,ua CCTK_REAL, dimension(2,2,2) :: d1a,gammad,gamma CCTK_REAL, dimension(2,2,2,2) :: d2a,gamma2 CCTK_REAL, allocatable, dimension(:,:) :: rr,xa,ya,za CCTK_REAL, allocatable, dimension(:,:) :: txx,tyy,tzz,txy,txz,tyz CCTK_REAL, allocatable, dimension(:,:) :: g11,g22,g12 CCTK_REAL, allocatable, dimension(:,:) :: gaussian character*200 gaussf c Reduction related things INTEGER reduction_handle ! Declarations for macros. #include "CactusEinstein/Einstein/src/Einstein.h" #include "CactusEinstein/Einstein/src/macro/UPPERMET_declare.h" #include "CactusEinstein/Einstein/src/macro/CHR2_declare.h" #include "CactusEinstein/Einstein/src/macro/RICCI_declare.h" #include "CactusEinstein/Einstein/src/macro/TRRICCI_declare.h" ! Data. data firstgau / .true. / data firstcal(1) / .true. / data firstcal(2) / .true. / data firstcal(3) / .true. / data firstcal(4) / .true. / save firstcal save firstgau ! Description of variables: ! ! i,j,k, ! l,m,n,p Counters. ! ! npoints Number of points to interpolate. ! ! error1 Different from zero if radius is negative. ! error2 Different from zero if outside computational domain. ! ! theta Latitude. ! phi Longitude. ! ! cost cos(theta) ! sint sin(theta) ! cosp cos(phi) ! sinp sin(phi) ! ! dtheta Grid spacing in theta. ! dphi Grid spacing in phi. ! dtp dtheta*dphi. ! ! idtheta 1/(2 dtheta) ! idphi 1/(2 dphi) ! ! xp x coordinate from surface centre. ! yp y coordinate from surface centre. ! zp z coordinate from surface centre. ! ! rp Radius. ! rr Radius array. ! ! xa Array with x coordinates from grid centre. ! ya Array with y coordinates from grid centre. ! za Array with z coordinates from grid centre. ! ! xmn Minimum value of x over the grid. ! xmx Maximum value of x over the grid. ! ! ymn Minimum value of y over the grid. ! ymx Maximum value of y over the grid. ! ! zmn Minimum value of z over the grid. ! zmx Maximum value of z over the grid. ! ! txx Interpolated gxx metric component. ! tyy Interpolated gyy metric component. ! tzz Interpolated gzz metric component. ! txy Interpolated gxy metric component. ! txz Interpolated gxz metric component. ! tyz Interpolated gyz metric component. ! ! trr 3-metric component {r,r}. ! ttt 3-metric component {theta,theta}. ! tpp 3-metric component {phi,phi}. ! trt 3-metric component {r,theta}. ! trp 3-metric component {r,phi}. ! ttp 3-metric component {theta,phi}. ! ! ft dr/dtheta ! fp dr/dphi ! ! g11 2-metric component {theta,theta}. ! g22 2-metric component {phi,phi}. ! g12 2-metric component {theta,phi}. ! ! circ_eq Equatorial circumference. ! meri_p1 Length of meridian at phi=0. ! meri_p2 Length of meridian at phi=pi/2. ! ! gaussian Gaussian curvature. ! ! ug 3x3 array with inverse 3-metric. ! ! xi_ij 3x3 array for extrinsic curvature of level sets ! of horizon function. ! trxi Trace of xi. ! xi2 Square of xi. ! ! ga_ij 2x2 array with angular metric, ! ua_ij 2x2 array with inverse angular metric. ! ! deta Determinant of angular metric. ! ideta 1/deta ! ! d1a_ijk d ga ! i kl ! ! d2a_ijkl d d ga ! i j kl ! ! gammad Christoffel symbol with 3 indices down. ! ! gamma Christoffel symbol with first index up. ! ! gamma2 Product of Christoffel symbols. ! ************************** ! *** DEFINE NUMBERS *** ! ************************** zero = 0.0D0 half = 0.5D0 one = 1.0D0 two = 2.0D0 three = 3.0D0 four = 4.0D0 pi = 3.141592654D0 ! ********************************************** ! * Get the reduction handle for sum operation * ! ********************************************** call CCTK_ReductionArrayHandle(reduction_handle,"sum"); if (reduction_handle.lt.0) then call CCTK_WARN(1,"Cannot get reduction handle.") endif ! ************************************** ! *** ALLOCATE MEMORY FOR ARRAYS *** ! ************************************** if (myproc.eq.0) then allocate(rr(0:ntheta,0:nphi)) allocate(xa(0:ntheta,0:nphi)) allocate(ya(0:ntheta,0:nphi)) allocate(za(0:ntheta,0:nphi)) allocate(txx(0:ntheta,0:nphi)) allocate(tyy(0:ntheta,0:nphi)) allocate(tzz(0:ntheta,0:nphi)) allocate(txy(0:ntheta,0:nphi)) allocate(txz(0:ntheta,0:nphi)) allocate(tyz(0:ntheta,0:nphi)) allocate(g11(0:ntheta,0:nphi)) allocate(g22(0:ntheta,0:nphi)) allocate(g12(0:ntheta,0:nphi)) allocate(gaussian(0:ntheta,0:nphi)) else allocate(rr(0:0,0:0)) allocate(xa(0:0,0:0)) allocate(ya(0:0,0:0)) allocate(za(0:0,0:0)) allocate(txx(0:0,0:0)) allocate(tyy(0:0,0:0)) allocate(tzz(0:0,0:0)) allocate(txy(0:0,0:0)) allocate(txz(0:0,0:0)) allocate(tyz(0:0,0:0)) allocate(g11(0:0,0:0)) allocate(g22(0:0,0:0)) allocate(g12(0:0,0:0)) allocate(gaussian(0:0,0:0)) end if rr = zero xa = zero ya = zero za = zero txx = zero tyy = zero tzz = zero txy = zero txz = zero tyz = zero g11 = zero g22 = zero g12 = zero gaussian = zero ! ********************************* ! *** INITIALIZE PARAMETERS *** ! ********************************* ! Find grid boundaries. call CCTK_ReductionHandle(handle,"minimum") call CCTK_CoordIndex(x_index,"x") call CCTK_CoordIndex(y_index,"y") call CCTK_CoordIndex(z_index,"z") call CCTK_Reduce(ierror,cctkGH,-1,handle,1, . CCTK_VARIABLE_REAL,minval,3, . x_index,y_index,z_index) xmn = minval(1) ymn = minval(2) zmn = minval(3) call CCTK_ReductionHandle(handle,"maximum") call CCTK_CoordIndex(x_index,"x") call CCTK_CoordIndex(y_index,"y") call CCTK_CoordIndex(z_index,"z") call CCTK_Reduce(ierror,cctkGH,-1,handle,1, . CCTK_VARIABLE_REAL,maxval,3, . x_index,y_index,z_index) xmx = maxval(1) ymx = maxval(2) zmx = maxval(3) ! Find {dtheta,dphi,dtp}. dtheta = pi/dble(ntheta) if (cartoon) then dphi = zero dtp = dtheta idphi = one idtheta = half/dtheta else dphi = two*pi/dble(nphi) if (refz) dtheta = half*dtheta if (refx.and.refy) then dphi = 0.25D0*dphi else if (refx) then dphi = half*dphi call CCTK_INFO('This combination of symmetries has not been properly implemented yet!') else if (refy) then dphi = half*dphi end if dtp = dtheta*dphi idtheta = half/dtheta idphi = half/dphi end if ! ******************************************************* ! *** FIND GAUSSIAN CURVATURE AS 3D GRID FUNCTION *** ! ******************************************************* ! Here I calculate the gaussian curvature of the level sets ! of the horizon function as a 3D grid function. Later I will ! interpolate it onto the horizon. This is cleaner than my old ! method (see below). It depends on less interpolations, and ! the interpolations are done only after all derivatives have ! been calculated. do k=2,nz-1 do j=2,ny-1 do i=2,nx-1 ! Find covariant normal unit vector. aux = one/ahfgradn(i,j,k) nlx = ahfgradx(i,j,k)*aux nly = ahfgrady(i,j,k)*aux nlz = ahfgradz(i,j,k)*aux ! Find upper spatial metric using standard macro. #include "CactusEinstein/Einstein/src/macro/UPPERMET_guts.h" ug(1,1) = UPPERMET_UXX ug(2,2) = UPPERMET_UYY ug(3,3) = UPPERMET_UZZ ug(1,2) = UPPERMET_UXY ug(1,3) = UPPERMET_UXZ ug(2,3) = UPPERMET_UYZ ug(2,1) = ug(1,2) ug(3,1) = ug(1,3) ug(3,2) = ug(2,3) ! Find contravariant normal unit vector. nux = ug(1,1)*nlx + ug(1,2)*nly + ug(1,3)*nlz nuy = ug(2,1)*nlx + ug(2,2)*nly + ug(2,3)*nlz nuz = ug(3,1)*nlx + ug(3,2)*nly + ug(3,3)*nlz ! Find Christoffel symbols using standard macros. #include "CactusEinstein/Einstein/src/macro/CHR2_guts.h" ! Find extrinsic curvature of the 2-surfaces defined ! as the level sets of the horizon function. The ! extrinsic curvature is defined as: ! __ ! xi = - \/ n ! ij (i j) ! __ ! where \/ is the 3-dimensional covariant derivative and ! n_i is the unit normal vector to the level sets of the ! horizon function. xi(1,1) = - (ahfgradx(i+1,j,k)/ahfgradn(i+1,j,k) . - ahfgradx(i-1,j,k)/ahfgradn(i-1,j,k))/(two*dx) . + (CHR2_XXX*nlx + CHR2_YXX*nly + CHR2_ZXX*nlz) xi(2,2) = - (ahfgrady(i,j+1,k)/ahfgradn(i,j+1,k) . - ahfgrady(i,j-1,k)/ahfgradn(i,j-1,k))/(two*dy) . + (CHR2_XYY*nlx + CHR2_YYY*nly + CHR2_ZYY*nlz) xi(3,3) = - (ahfgradz(i,j,k+1)/ahfgradn(i,j,k+1) . - ahfgradz(i,j,k-1)/ahfgradn(i,j,k-1))/(two*dz) . + (CHR2_XZZ*nlx + CHR2_YZZ*nly + CHR2_ZZZ*nlz) xi(1,2) = - (ahfgradx(i,j+1,k)/ahfgradn(i,j+1,k) . - ahfgradx(i,j-1,k)/ahfgradn(i,j-1,k))/(four*dy) . - (ahfgrady(i+1,j,k)/ahfgradn(i+1,j,k) . - ahfgrady(i-1,j,k)/ahfgradn(i-1,j,k))/(four*dx) . + (CHR2_XXY*nlx + CHR2_YXY*nly + CHR2_ZXY*nlz) xi(1,3) = - (ahfgradx(i,j,k+1)/ahfgradn(i,j,k+1) . - ahfgradx(i,j,k-1)/ahfgradn(i,j,k-1))/(four*dz) . - (ahfgradz(i+1,j,k)/ahfgradn(i+1,j,k) . - ahfgradz(i-1,j,k)/ahfgradn(i-1,j,k))/(four*dx) . + (CHR2_XXZ*nlx + CHR2_YXZ*nly + CHR2_ZXZ*nlz) xi(2,3) = - (ahfgrady(i,j,k+1)/ahfgradn(i,j,k+1) . - ahfgrady(i,j,k-1)/ahfgradn(i,j,k-1))/(four*dz) . - (ahfgradz(i,j+1,k)/ahfgradn(i,j+1,k) . - ahfgradz(i,j-1,k)/ahfgradn(i,j-1,k))/(four*dy) . + (CHR2_XYZ*nlx + CHR2_YYZ*nly + CHR2_ZYZ*nlz) xi(2,1) = xi(1,2) xi(3,1) = xi(1,3) xi(3,2) = xi(2,3) ! Find trace of xi. aux = zero do l=1,3 do m=1,3 aux = aux + ug(l,m)*xi(l,m) end do end do trxi = aux ! Find square of xi. aux = zero do l=1,3 do m=1,3 do n=1,3 do p=1,3 aux = aux + ug(l,m)*ug(n,p)*xi(l,n)*xi(m,p) end do end do end do end do xi2 = aux ! Find 3-Ricci and 3-Ricci scalar using standard macros. #include "CactusEinstein/Einstein/src/macro/RICCI_guts.h" #include "CactusEinstein/Einstein/src/macro/TRRICCI_guts.h" ! Find 2-Ricci scalar using the contracted Gauss-Codazzi ! relations: ! ! (2) (3) a b (3) 2 ab ! R = R - 2 n n R + (tr xi) - xi xi ! ab ab aux = TRRICCI_TRRICCI + trxi**2 - xi2 aux = aux - two*(nux**2*RICCI_RXX + nuy**2*RICCI_RYY . + nuz**2*RICCI_RZZ + two*(nux*nuy*RICCI_RXY . + nux*nuz*RICCI_RXZ + nuy*nuz*RICCI_RYZ)) ! Find gaussian curvature. The gaussian curvature is ! defined as R/2, with R the Ricci scalar of the surfaces. ahfgauss(i,j,k) = half*aux ! Undefine macros! #include "CactusEinstein/Einstein/src/macro/UPPERMET_undefine.h" #include "CactusEinstein/Einstein/src/macro/CHR2_undefine.h" #include "CactusEinstein/Einstein/src/macro/RICCI_undefine.h" #include "CactusEinstein/Einstein/src/macro/TRRICCI_undefine.h" end do end do end do ! Boundaries. ahfgauss(1,:,:) = ahfgauss(2,:,:) ahfgauss(:,1,:) = ahfgauss(:,2,:) ahfgauss(:,:,1) = ahfgauss(:,:,2) ahfgauss(nx,:,:) = ahfgauss(nx-1,:,:) ahfgauss(:,ny,:) = ahfgauss(:,ny-1,:) ahfgauss(:,:,nz) = ahfgauss(:,:,nz-1) ! Synchronize. call CCTK_SyncGroup(cctkGH,"ahfinder::ahfinder_gauss") ! call synconefunc(ahfgauss) ! ************************************* ! *** FIND INTERPOLATING POINTS *** ! ************************************* ! Initialize {error1,error2}. error1 = 0 error2 = 0 ! Notice that everything here is done on processor zero! if (myproc.eq.0) then do j=0,nphi do i=0,ntheta ! Find {theta,phi}. theta = dtheta*dble(i) phi = dphi*dble(j) ! Find sines and cosines. cost = cos(theta) sint = sin(theta) cosp = cos(phi) sinp = sin(phi) ! Find radius rp. rp = c0(0) do l=1+stepz,lmax,1+stepz rp = rp + c0(l)*LEGEN(l,0,cost) end do ! Notice how the sum over m is first. This will allow ! me to use the recursion relations to avoid having to ! start from scratch every time. Also, I sum over all ! l s even if I do not want some terms. This is ! because in order to use the recursion relations I ! need all polynomials. if (nonaxi) then do m=1,lmax aux = dble(m)*phi sina = sin(aux) cosa = cos(aux) do l=m,lmax aux = LEGEN(l,m,cost) rp = rp + aux*cc(l,m)*cosa if (.not.refy) then rp = rp + aux*cs(l,m)*sina end if end do end do end if ! Check for negative radius. if (rp.le.zero) then error1 = 1 end if ! Save radius array. rr(i,j) = rp ! Find {xa,ya,za} xa(i,j) = rp*sint*cosp + xc ya(i,j) = rp*sint*sinp + yc za(i,j) = rp*cost + zc ! Check if we are within bounds. xp = xa(i,j) yp = ya(i,j) zp = za(i,j) if ((xp.gt.xmx).or.(xp.lt.xmn).or. . (yp.gt.ymx).or.(yp.lt.ymn).or. . (zp.gt.zmx).or.(zp.lt.zmn)) then error2 = 1 end if end do end do else xa = half*(xmx+xmn) ya = half*(ymx+ymn) za = half*(zmx+zmn) end if ! Reduce the errors across processors (all processors must ! know about this since all will participate on the interpolation ! below). call CCTK_ReduceLocalScalar(ierror, cctkGH, -1, reduction_handle, $ error1, i, CCTK_VARIABLE_REAL) if (ierror.ne.0) then call CCTK_WARN(1,"Reduction failed!"); endif error1 = i call CCTK_ReduceLocalScalar(ierror, cctkGH, -1, reduction_handle, $ error2, i, CCTK_VARIABLE_REAL) if (ierror.ne.0) then call CCTK_WARN(1,"Reduction failed!"); endif error2 = i cc#ifdef MPI cc call mpi_allreduce(error1,i,1,MPI_INTEGER, cc . MPI_SUM,MPI_COMM_WORLD,ierror) cc error1 = i cc call mpi_allreduce(error2,i,1,MPI_INTEGER, cc . MPI_SUM,MPI_COMM_WORLD,ierror) cc error2 = i cc#endif ! If there was an error then return from the subroutine ! (but remember to deallocate the arrays first!). if ((error1.ne.0).or.(error2.ne.0)) then goto 10 end if ! *********************** ! *** INTERPOLATE *** ! *********************** if (myproc.eq.0) then npoints = (ntheta+1)*(nphi+1) else npoints = 1 end if ! new interp call CCTK_GetInterpHandle (handle, "simple") call CCTK_VarIndex (gxx_index, "einstein::gxx") call CCTK_VarIndex (gyy_index, "einstein::gyy") call CCTK_VarIndex (gzz_index, "einstein::gzz") call CCTK_VarIndex (gxy_index, "einstein::gxy") call CCTK_VarIndex (gxz_index, "einstein::gxz") call CCTK_VarIndex (gyz_index, "einstein::gyz") call CCTK_VarIndex (ahfgauss_index, "ahfinder::ahfgauss") call CCTK_InterpGF (ierror, cctkGH, handle, npoints, 3, 7, 7, $ xa, ya, za, $ CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, $ CCTK_VARIABLE_REAL, $ gxx_index,gyy_index,gzz_index, $ gxy_index,gxz_index,gyz_index,ahfgauss_index, $ txx,tyy,tzz,txy,txz,tyz,gaussian, $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL ) ! ********************************* ! *** FIND SPHERICAL METRIC *** ! ********************************* if (myproc.eq.zero) then do j=0,nphi do i=0,ntheta ! Find {theta,phi}. theta = dtheta*dble(i) phi = dphi*dble(j) ! Find {rp}. rp = rr(i,j) ! Find sines and cosines. cost = cos(theta) sint = sin(theta) cosp = cos(phi) sinp = sin(phi) ! Find metric in spherical coordinates. trr = sint**2*(txx(i,j)*cosp**2 . + tyy(i,j)*sinp**2) + tzz(i,j)*cost**2 . + two*sint*(txy(i,j)*sint*cosp*sinp . + cost*(txz(i,j)*cosp + tyz(i,j)*sinp)) ttt = rp**2*(cost**2*(txx(i,j)*cosp**2 . + tyy(i,j)*sinp**2) + tzz(i,j)*sint**2 . + two*cost*(txy(i,j)*cost*cosp*sinp . - sint*(txz(i,j)*cosp + tyz(i,j)*sinp))) tpp = (rp*sint)**2*(txx(i,j)*sinp**2 + tyy(i,j)*cosp**2 . - two*txy(i,j)*cosp*sinp) trt = rp*(sint*cost*(txx(i,j)*cosp**2 + tyy(i,j)*sinp**2 . - tzz(i,j) + two*txy(i,j)*sinp*cosp) . + (cost**2-sint**2)*(txz(i,j)*cosp + tyz(i,j)*sinp)) trp = rp*sint*(sint*(cosp*sinp*(tyy(i,j) - txx(i,j)) . + txy(i,j)*(cosp**2 - sinp**2)) . + cost*(cosp*tyz(i,j) - sinp*txz(i,j))) ttp = rp**2*sint*(cost*(sinp*cosp*(tyy(i,j) - txx(i,j)) . + (cosp**2 - sinp**2)*txy(i,j)) . + sint*(sinp*txz(i,j) - cosp*tyz(i,j))) ! Find derivatives {ft,fp}. For interior points I ! use centered differences, and for the boundary points ! I use second order one-sided differences. if ((i.ne.0).and.(i.ne.ntheta)) then ft = idtheta*(rr(i+1,j) - rr(i-1,j)) else if (i.eq.0) then ft = - idtheta*(three*rr(0,j) . - four*rr(1,j) + rr(2,j)) else ft = + idtheta*(three*rr(ntheta,j) . - four*rr(ntheta-1,j) + rr(ntheta-2,j)) end if if (nonaxi) then if ((j.ne.0).and.(j.ne.nphi)) then fp = idphi*(rr(i,j+1) - rr(i,j-1)) else if (j.eq.0) then fp = - idphi*(three*rr(i,0) . - four*rr(i,1) + rr(i,2)) else fp = + idphi*(three*rr(i,nphi) . - four*rr(i,nphi-1) + rr(i,nphi-2)) end if else fp = zero end if ! Find induced metric on the surface. g11(i,j) = ttt + trr*ft**2 + two*trt*ft g22(i,j) = tpp + trr*fp**2 + two*trp*fp g12(i,j) = ttp + trr*ft*fp + trp*ft + trt*fp end do end do ! ************************************ ! *** CALCULATE CIRCUMFERENCES *** ! ************************************ ! Find length of equator and meridians circ_eq = zero meri_p1 = zero meri_p2 = zero ! Equatorial circumference. if (refz) then do j=0,nphi-1 circ_eq = circ_eq + dphi . *dsqrt(half*(g22(ntheta,j) + g22(ntheta,j+1))) end do else aux = half*pi/dtheta i = int(aux) aux = aux - int(aux) if (cartoon) then circ_eq = 6.283185307D0*dsqrt(g22(i,0)) else do j=0,nphi-1 circ_eq = circ_eq + dphi . *((one-aux)*dsqrt(half*(g22(i,j) + g22(i,j+1))) . + aux*dsqrt(half*(g22(i+1,j) + g22(i+1,j+1)))) end do end if end if ! Meridians do i=0,ntheta-1 meri_p1 = meri_p1 + dtheta . *dsqrt(half*(g11(i,0) + g11(i+1,0))) end do if (cartoon) then meri_p2 = meri_p1 else if (refx.and.refy) then do i=0,ntheta-1 meri_p2 = meri_p2 + dtheta . *dsqrt(half*(g11(i,nphi) + g11(i+1,nphi))) end do else aux = half*pi/dphi j = int(aux) aux = aux - int(aux) do i=0,ntheta-1 meri_p2 = meri_p2 + dtheta . *((one-aux)*dsqrt(half*(g11(i,j) + g11(i,j))) . + aux*dsqrt(half*(g11(i,j+1) + g11(i,j+1)))) end do end if ! Rescale the integrals according to symmetries. if (.not.cartoon) then if (refx) circ_eq = two*circ_eq if (refy) circ_eq = two*circ_eq end if if (refz) then meri_p1 = two*meri_p1 meri_p2 = two*meri_p2 end if end if ! ************************************************* ! *** OLD CALCULATION OF GAUSSIAN CURVATURE *** ! ************************************************* ! This was the old calculation of the gaussian curvature. ! It is correct as far as I know, but it depends on taking ! second derivatives of the interpolated metric, and this ! seems not to behave well numerically. I did not want to ! delete the code since it might be useful in the future, ! so I just jump over it. if (.false.) then ! The Gaussian curvature of the surface is defined as R/2 ! where R is the Ricci scalar of the two-geometry. I can ! now calculate the Ricci scalar using the spherical metric ! components {ttt,tpp,ttp}. ! Initialize arrays. gaussian = zero ga(1,1) = one ga(2,2) = one ga(1,2) = zero ga(2,1) = zero ua(1,1) = one ua(2,2) = one ua(1,2) = zero ua(2,1) = zero d1a = zero d2a = zero gamma = zero gamma2 = zero gammad = zero if (myproc.eq.zero) then ! Interior points. do j=1,nphi-1 do i=1,ntheta-1 ! Find angular metric. ga(1,1) = g11(i,j) ga(2,2) = g22(i,j) ga(1,2) = g12(i,j) ga(2,1) = g12(i,j) ! Find determinant of angular metric. deta = ga(1,1)*ga(2,2) - ga(1,2)**2 ideta = one/deta ! Find inverse angular metric. ua(1,1) = + ga(2,2)*ideta ua(2,2) = + ga(1,1)*ideta ua(1,2) = - ga(1,2)*ideta ua(2,1) = - ga(2,1)*ideta ! First derivatives of angular metric. d1a(1,1,1) = idtheta*(g11(i+1,j) - g11(i-1,j)) d1a(1,2,2) = idtheta*(g22(i+1,j) - g22(i-1,j)) d1a(1,1,2) = idtheta*(g12(i+1,j) - g12(i-1,j)) d1a(1,2,1) = d1a(1,1,2) if (nonaxi) then d1a(2,1,1) = idphi*(g11(i,j+1) - g11(i,j-1)) d1a(2,2,2) = idphi*(g22(i,j+1) - g22(i,j-1)) d1a(2,1,2) = idphi*(g12(i,j+1) - g12(i,j-1)) d1a(2,2,1) = d1a(2,1,2) else d1a(2,1,1) = zero d1a(2,2,2) = zero d1a(2,1,2) = zero d1a(2,2,1) = zero end if ! Second derivatives of angular metric. d2a(1,1,1,1) = four*idtheta**2*(g11(i+1,j) . - two*g11(i,j) + g11(i-1,j)) d2a(1,1,2,2) = four*idtheta**2*(g22(i+1,j) . - two*g22(i,j) + g22(i-1,j)) d2a(1,1,1,2) = four*idtheta**2*(g12(i+1,j) . - two*g12(i,j) + g12(i-1,j)) if (nonaxi) then d2a(2,2,1,1) = four*idphi**2*(g11(i,j+1) . - two*g11(i,j) + g11(i,j-1)) d2a(2,2,2,2) = four*idphi**2*(g22(i,j+1) . - two*g22(i,j) + g22(i,j-1)) d2a(2,2,1,2) = four*idphi**2*(g12(i,j+1) . - two*g12(i,j) + g12(i,j-1)) d2a(1,2,1,1) = (g11(i+1,j+1) - g11(i+1,j-1) . - g11(i-1,j+1) + g11(i-1,j-1)) . *idtheta*idphi d2a(1,2,2,2) = (g22(i+1,j+1) - g22(i+1,j-1) . - g22(i-1,j+1) + g22(i-1,j-1)) . *idtheta*idphi d2a(1,2,1,2) = (g12(i+1,j+1) - g12(i+1,j-1) . - g12(i-1,j+1) + g12(i-1,j-1)) . *idtheta*idphi else d2a(2,2,1,1) = zero d2a(2,2,2,2) = zero d2a(2,2,1,2) = zero d2a(1,2,1,1) = zero d2a(1,2,2,2) = zero d2a(1,2,1,2) = zero end if d2a(1,1,2,1) = d2a(1,1,1,2) d2a(2,2,2,1) = d2a(2,2,1,2) d2a(1,2,2,1) = d2a(1,2,1,2) d2a(2,1,1,1) = d2a(1,2,1,1) d2a(2,1,2,2) = d2a(1,2,2,2) d2a(2,1,1,2) = d2a(1,2,1,2) d2a(2,1,2,1) = d2a(1,2,2,1) ! Add correction terms to derivatives of g_{phi,phi}. ! Near the axis, it turns out that the lower order ! contributions to the derivatives of g_{phi,phi} ! cancel out, and the gaussian curvature is made only ! of higher order terms that the centered finite difference ! approximations get wrong! So here I add higher order ! corrections to make sure that I get correctly the ! first two terms of the Taylor series (assuming that ! for small theta g_{phi,phi} = f(r,phi) sin(theta)^2). theta = dtheta*dble(i) d1a(1,2,2) = d1a(1,2,2) + (two/three)*sin(theta) . *(g22(i+1,j) - two*g22(i,j) + g22(i-1,j)) d2a(1,1,2,2) = d2a(1,1,2,2) + one/three . *(g22(i+1,j) - two*g22(i,j) + g22(i-1,j)) if (nonaxi) then d2a(1,2,2,2) = d2a(1,2,2,2) + (two/three)*sin(theta) . *((g22(i+1,j+1) - two*g22(i,j+1) + g22(i-1,j+1)) . - (g22(i+1,j-1) - two*g22(i,j-1) + g22(i-1,j-1))) . *idphi d2a(2,1,2,2) = d2a(1,2,2,2) end if ! Find Christoffel symbols. do l=1,2 do m=1,2 do n=1,2 gammad(l,m,n) = half*(d1a(m,n,l) + d1a(n,m,l) . - d1a(l,m,n)) end do end do end do do l=1,2 do m=1,2 do n=1,2 aux = zero do k=1,2 aux = aux + ua(l,k)*gammad(k,m,n) end do gamma(l,m,n) = aux end do end do end do do k=1,2 do l=1,2 do m=1,2 do n=1,2 aux = zero do p=1,2 aux = aux + gamma(p,k,l)*gammad(p,m,n) end do gamma2(k,l,m,n) = aux end do end do end do end do ! Find gaussian curvature. aux = zero do k=1,2 do l=1,2 do m=1,2 do n=1,2 aux = aux + ua(k,l)*ua(m,n) . *(d2a(k,n,l,m) - d2a(k,l,m,n) . + gamma2(k,m,l,n) - gamma2(k,l,m,n)) end do end do end do end do gaussian(i,j) = half*aux end do end do ! Boundaries. gaussian(0,:) = two*gaussian(1,:) - gaussian(2,:) gaussian(ntheta,:) = two*gaussian(ntheta-1,:) . - gaussian(ntheta-2,:) gaussian(:,0) = two*gaussian(:,1) - gaussian(:,2) gaussian(:,nphi) = two*gaussian(:,nphi-1) . - gaussian(:,nphi-2) end if end if ! ************************************ ! *** INITIALIZE FIRSTGAU FLAG *** ! ************************************ if (find3) then firstgau = firstcal(mfind+1) firstcal(mfind+1) = .false. else firstgau = firstcal(4) firstcal(4) = .false. end if ! *********************************** ! *** SAVE GAUSSIAN CURVATURE *** ! *********************************** if (.false.) then if (myproc.eq.0) then if (find3) then if (mfind.eq.0) then gaussf = filestr(1:nfile)//"/mahf_0.gauss" else if (mfind.eq.1) then gaussf = filestr(1:nfile)//"/mahf_1.gauss" else if (mfind.eq.2) then gaussf = filestr(1:nfile)//"/mahf_2.gauss" end if else gaussf = filestr(1:nfile)//"/mahf.gauss" end if ! On first call replace file. if (firstgau) then if (find3) then if (mfind.eq.2) then firstgau = .false. end if else firstgau = .false. end if open(1,file=gaussf,form='formatted', . status='replace') write(1,"(A20)") '# GAUSSIAN CURVATURE' write(1,"(A1)") '#' write(1,"(A35)") '# The data is written in a loop as:' write(1,"(A1)") '#' write(1,"(A17)") '# do i=0,ntheta-1' write(1,"(A18)") '# do j=0,nphi-1' write(1,"(A27)") '# write gaussian(i,j)' write(1,"(A11)") '# end do' write(1,"(A8)") '# end do' write(1,"(A1)") '#' write(1,"(A40)") '# theta and phi are subdivided uniformly' write(1,"(A26)") '# according to symmetries:' write(1,"(A1)") '#' write(1,"(A36)") '# phi=[0,2 pi] (refx=refy=.false.)' write(1,"(A44)") '# phi=[0,pi] (refx=.false., refy=.true.)' write(1,"(A35)") '# phi=[0,pi/2] (refx=refy=.true.)' write(1,"(A1)") '#' write(1,"(A31)") '# theta=[0,pi] (refz=.false.)' write(1,"(A30)") '# theta=[0,pi/2] (refz=.true.)' write(1,"(A1)") '#' write(1,"(A9,L1)") '# refx = ',refx write(1,"(A9,L1)") '# refy = ',refy write(1,"(A9,L1)") '# refz = ',refz write(1,"(A1)") '#' write(1,"(A11,I7)") '# ntheta = ',(ntheta+1) write(1,"(A11,I7)") '# nphi = ',(nphi+1) write(1,*) write(1,"(A11,I4)") '# Time step',cctk_iteration write(1,"(A6,ES11.3)") '# Time',cctk_time write(1,"(A6,I4)") '# Call',ahf_ncall ! On subsequent calls, append to file. else open(1,file=gaussf,form='formatted', . status='old',position='append') write(1,*) write(1,"(A11,I4)") '# Time step',cctk_iteration write(1,"(A6,ES11.3)") '# Time',cctk_time write(1,"(A6,I4)") '# Call',ahf_ncall end if ! Save data points. do i=0,ntheta do j=0,nphi write(1,"(ES11.3)") gaussian(i,j) end do end do ! Close file. close(1) end if endif ! ***************************** ! *** DEALLOCATE MEMORY *** ! ***************************** 10 continue deallocate(rr) deallocate(xa,ya,za) deallocate(txx,tyy,tzz,txy,txz,tyz) deallocate(g11,g22,g12) deallocate(gaussian) ! *************** ! *** END *** ! *************** end subroutine AHFinder_gau