diff options
Diffstat (limited to 'src/AHFinder_gau.F')
| -rw-r--r-- | src/AHFinder_gau.F | 1168 |
1 files changed, 1168 insertions, 0 deletions
diff --git a/src/AHFinder_gau.F b/src/AHFinder_gau.F new file mode 100644 index 0000000..ef3e6cb --- /dev/null +++ b/src/AHFinder_gau.F @@ -0,0 +1,1168 @@ +c/*@@ +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 + +#ifdef MPI + include 'mpif.h' +#endif + + logical firstgau + logical firstcal(4) + + integer handle,x_index,y_index,z_index + CCTK_INT ierror + CCTK_REAL maxval(3),minval(3) + + + 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 + +! Declarations for macros. + +!#include "./thorn_StandardEinstein/src/macro/UPPERMET_declare.h" +!#include "./thorn_StandardEinstein/src/macro/CHR2_declare.h" +!#include "./thorn_StandardEinstein/src/macro/RICCI_declare.h" +!#include "./thorn_StandardEinstein/src/macro/TRRICCI_declare.h" + +#include "../../Einstein/src/macro/UPPERMET_declare.h" +#include "../../Einstein/src/macro/CHR2_declare.h" +#include "../../Einstein/src/macro/RICCI_declare.h" +#include "../../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 + + +! ************************************** +! *** 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 + + +! ********************************* +! *** INITIALIZE PARAMETERS *** +! ********************************* + +! Find grid boundaries. +! call CCTK_ReductionHandle(handle,"minimum") +! call CCTK_GetCoordIndex(x_index,"x") +! call CCTK_GetCoordIndex(y_index,"y") +! call CCTK_GetCoordIndex(z_index,"z") + +! call CCTK_Reduce(cctkGH,ierror,-1,handle,3, +! . CCTK_VARIABLE_REAL,minval,3, +! . x_index,y_index,z_index) + +! xmn = gf_min(x) +! ymn = gf_min(y) +! zmn = gf_min(z) + minval = -9.45D0 + + xmn = minval(1) + ymn = minval(2) + zmn = minval(3) + +! call CCTK_ReductionHandle(handle,"maximum") +! call CCTK_GetCoordIndex(x_index,"x") +! call CCTK_GetCoordIndex(y_index,"y") +! call CCTK_GetCoordIndex(z_index,"z") + +! call CCTK_Reduce(cctkGH,ierror,-1,handle,3, +! . CCTK_VARIABLE_REAL,maxval,3, +! . x_index,y_index,z_index) + +! xmx = gf_max(x) +! ymx = gf_max(y) +! zmx = gf_max(z) + maxval =9.45D0 + + 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 + print *, 'This combination of symmetries has not' + print *, '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 "../../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 "../../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 "../../Einstein/src/macro/RICCI_guts.h" +#include "../../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 "../../Einstein/src/macro/UPPERMET_undefine.h" +#include "../../Einstein/src/macro/CHR2_undefine.h" +#include "../../Einstein/src/macro/RICCI_undefine.h" +#include "../../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. + +#ifdef MPI + call synconefunc(ahfgauss) +#endif + +! ************************************* +! *** 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 don't 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). + +#ifdef MPI + + call mpi_allreduce(error1,i,1,MPI_INTEGER, + . MPI_SUM,MPI_COMM_WORLD,ierror) + error1 = i + call mpi_allreduce(error2,i,1,MPI_INTEGER, + . MPI_SUM,MPI_COMM_WORLD,ierror) + error2 = i + +#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 + + call start_getpointsgroup() + call getpoints(gxx,xa,ya,za,txx,npoints) + call getpoints(gyy,xa,ya,za,tyy,npoints) + call getpoints(gzz,xa,ya,za,tzz,npoints) + call getpoints(gxy,xa,ya,za,txy,npoints) + call getpoints(gxz,xa,ya,za,txz,npoints) + call getpoints(gyz,xa,ya,za,tyz,npoints) + call getpoints(ahfgauss,xa,ya,za,gaussian,npoints) + call cleanup_getpointsgroup() + + +! ********************************* +! *** 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 didn't 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 (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 + + +! ***************************** +! *** 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 + + |