/*@@ @file AHFinder_gau.F @date October 1998 @author Miguel Alcubierre @desc Find gaussian curvature of surface. The gaussian curvature is defined as R/2, where R is the Ricci scalar of the induced 2-geometry of the surface. As an extra "goody", this routine also integrates the equatorial circumference of the horizon and two polar circumferences at phi=0 and phi=pi/2. @enddesc @@*/ #include "cctk.h" #include "cctk_Parameters.h" #include "cctk_Arguments.h" #include "cctk_Functions.h" subroutine AHFinder_gau(CCTK_ARGUMENTS) use AHFinder_dat implicit none DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS logical firstgau logical firstcal(4) integer i,j,k,l,m,n,p integer npoints,vindex integer param_table_handle,interp_handle,coord_system_handle,sum_handle integer ierror character(len=128) :: operator CCTK_INT nchars CCTK_INT rerror,error1,error2 character(30) options_string CCTK_REAL LEGEN CCTK_REAL theta,phi,xp,yp,zp,rp CCTK_REAL cost,sint,cosp,sinp CCTK_REAL dtheta,dphi,dtp,idtheta,idphi,phistart CCTK_REAL trr,ttt,tpp,trt,trp,ttp,ft,fp 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 c CCTK_REAL dt,idx,idy,idz c CCTK_REAL i2dx,i2dy,i2dz,idxx,idyy,idzz,idxy,idxz,idyz 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_POINTER, dimension(3) :: interp_coords CCTK_INT, dimension(7) :: in_array_indices CCTK_POINTER, dimension(7) :: out_arrays CCTK_INT, dimension(7) :: out_array_type_codes 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, dimension(3) :: ahf_centroid, rahf_centroid CCTK_REAL :: m1p, m2p, x1p, x2p, y1p, y2p, z1p, z2p CCTK_REAL :: r12, r22, sigma CCTK_REAL, dimension(nx,ny,nz) :: fac CCTK_INT :: use_att, apower logical :: use_rot_att logical :: str_comp logical :: gaussf_exists external str_comp ! CCTK_INT :: CCTK_IsFunctionAliased character(len=32) :: tmpstr character(len=200) :: gaussf ! Reduction variables CCTK_POINTER_TO_CONST, dimension(3) :: input_array CCTK_POINTER, dimension(1) :: reduction_value CCTK_POINTER, dimension(3) :: reduction_value3 CCTK_INT, dimension(1) ::input_array_type CCTK_INT, dimension(3) ::input_array_type3 ! Declarations for macros. c#include "EinsteinBase/ADMMacros/src/macro/ADM_Spacing_declare.h" CCTK_REAL :: dt CCTK_REAL :: idx, idy, idz CCTK_REAL :: i2dx, i2dy, i2dz CCTK_REAL :: i12dx, i12dy, i12dz CCTK_REAL :: idxx, idxy, idxz, idyy, idyz, idzz CCTK_REAL :: i12dxx, i12dyy, i12dzz CCTK_REAL :: i36dxy, i36dxz, i36dyz #include "EinsteinBase/ADMMacros/src/macro/UPPERMET_declare.h" #include "EinsteinBase/ADMMacros/src/macro/CHR2_declare.h" #include "EinsteinBase/ADMMacros/src/macro/RICCI_declare.h" #include "EinsteinBase/ADMMacros/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. ! ! phistart Origin for phi (normally 0, but for some symmetries -pi/2). ! ! 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. ! ! 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}. ! ! 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 = acos(-one) ! Number of points to interpolate. if (myproc.eq.0) then npoints = (ntheta+1)*(nphi+1) else npoints = 1 end if ! ****************************************************** ! *** Get the reduction handle for sum operation *** ! ****************************************************** call CCTK_LocalArrayReductionHandle(sum_handle,"sum") if (sum_handle.lt.0) then call CCTK_WARN(1,"Cannot get handle for sum reduction ! Forgot to activate an implementation providing reduction operators ??") end if ! ************************************** ! *** ALLOCATE MEMORY FOR ARRAYS *** ! ************************************** if (myproc.eq.0) then allocate(rr(1:ntheta+1,1:nphi+1)) allocate(xa(1:ntheta+1,1:nphi+1)) allocate(ya(1:ntheta+1,1:nphi+1)) allocate(za(1:ntheta+1,1:nphi+1)) allocate(txx(1:ntheta+1,1:nphi+1)) allocate(tyy(1:ntheta+1,1:nphi+1)) allocate(tzz(1:ntheta+1,1:nphi+1)) allocate(txy(1:ntheta+1,1:nphi+1)) allocate(txz(1:ntheta+1,1:nphi+1)) allocate(tyz(1:ntheta+1,1:nphi+1)) allocate(g11(1:ntheta+1,1:nphi+1)) allocate(g22(1:ntheta+1,1:nphi+1)) allocate(g12(1:ntheta+1,1:nphi+1)) else allocate(rr(1,1)) allocate(xa(1,1)) allocate(ya(1,1)) allocate(za(1,1)) allocate(txx(1,1)) allocate(tyy(1,1)) allocate(tzz(1,1)) allocate(txy(1,1)) allocate(txz(1,1)) allocate(tyz(1,1)) allocate(g11(1,1)) allocate(g22(1,1)) allocate(g12(1,1)) 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 {phistart,dtheta,dphi,dtp}. phistart = zero 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 phistart = - half*pi else if (refy) then dphi = half*dphi end if dtp = dtheta*dphi idtheta = half/dtheta idphi = half/dphi end if ! Initialize circumferences. circ_eq = zero meri_p1 = zero meri_p2 = zero #include "EinsteinBase/ADMMacros/src/macro/ADM_Spacing.h" ! ******************************************************* ! *** 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 "EinsteinBase/ADMMacros/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 "EinsteinBase/ADMMacros/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 "EinsteinBase/ADMMacros/src/macro/RICCI_guts.h" #include "EinsteinBase/ADMMacros/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 - 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)) aux = aux + trxi**2 - xi2 ! 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 "EinsteinBase/ADMMacros/src/macro/UPPERMET_undefine.h" #include "EinsteinBase/ADMMacros/src/macro/CHR2_undefine.h" #include "EinsteinBase/ADMMacros/src/macro/RICCI_undefine.h" #include "EinsteinBase/ADMMacros/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(ierror,cctkGH,"ahfinder::ahfinder_gauss") ! ************************************* ! *** FIND INTERPOLATING POINTS *** ! ************************************* ! Initialize {error1,error2}. error1 = 0 error2 = 0 ! Notice that everything here is done on processor zero! if (myproc.eq.0) then avgx = 0.0D0 avgy = 0.0D0 avgz = 0.0D0 do j=1,nphi+1 do i=1,ntheta+1 ! Find {theta,phi}. theta = dtheta*dble(i-1) phi = dphi*dble(j-1) + phistart ! 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 ! Add to sums for average position. avgx = avgx + xp avgy = avgy + yp avgz = avgz + zp end do end do ! Find position of horizon centroid. avgx = avgx/dble(npoints) avgy = avgy/dble(npoints) avgz = avgz/dble(npoints) ahf_centroid(1) = avgx ahf_centroid(2) = avgy ahf_centroid(3) = avgz ! Take symmetries into account. if (refx) ahf_centroid(1) = xc if (refy) ahf_centroid(2) = yc if (refz) ahf_centroid(3) = zc ! Other processors. else xa = half*(xmx+xmn) ya = half*(ymx+ymn) za = half*(zmx+zmn) ahf_centroid(1) = zero ahf_centroid(2) = zero ahf_centroid(3) = zero end if ! Reduce the errors across processors (all processors must ! know about this since all will participate on the interpolation ! below). input_array_type(1) = CCTK_VARIABLE_INT reduction_value(1)= CCTK_PointerTo(rerror) input_array(1) = CCTK_PointerTo(error1) call CCTK_ReduceArraysGlobally(ierror, cctkGH, -1,sum_handle, -1, . 1, input_array,0, 0, input_array_type, 1, . input_array_type, reduction_value) if (ierror.ne.0) then call CCTK_WARN(1,"Reduction failed!") end if error1 = rerror input_array(1) = CCTK_PointerTo(error2) call CCTK_ReduceArraysGlobally(ierror, cctkGH, -1,sum_handle, -1, . 1, input_array,0, 0, input_array_type, 1, . input_array_type, reduction_value) if (ierror.ne.0) then call CCTK_WARN(1,"Reduction failed!") end if error2 = rerror input_array(1) = CCTK_PointerTo(ahf_centroid(1)) input_array(2) = CCTK_PointerTo(ahf_centroid(2)) input_array(3) = CCTK_PointerTo(ahf_centroid(3)) reduction_value3(1) = CCTK_PointerTo(rahf_centroid(1)) reduction_value3(2) = CCTK_PointerTo(rahf_centroid(2)) reduction_value3(3) = CCTK_PointerTo(rahf_centroid(3)) input_array_type3(1) = CCTK_VARIABLE_REAL input_array_type3(2) = CCTK_VARIABLE_REAL input_array_type3(3) = CCTK_VARIABLE_REAL call CCTK_ReduceArraysGlobally(ierror, cctkGH, -1,sum_handle, -1, . 3, input_array,0, 0, input_array_type3, 3, . input_array_type3, reduction_value3) if (ierror.ne.0) then call CCTK_WARN(1,"Reduction failed!") end if if ( ( horizon_to_announce_centroid .gt. 0 ) .and. . ( horizon_to_announce_centroid .eq. mfind ) ) then call CCTK_IsFunctionAliased(ierror, . "SetDriftCorrectPosition") if ( ierror .eq. 1 ) then call CCTK_INFO("Announcing centroid to DriftCorrect") call SetDriftCorrectPosition ( cctkGH, ahf_centroid(1), . ahf_centroid(2), . ahf_centroid(3) ) end if end if if ( ( horizon_to_output_centroid .gt. 0 ) .and. . ( horizon_to_output_centroid .eq. mfind ) ) then ahf_centroid_x = rahf_centroid(1) ahf_centroid_y = rahf_centroid(2) ahf_centroid_z = rahf_centroid(3) if (verbose) then write(*,*) write(6,*) 'AHFinder_gau: ahf_centroid_x = ', . ahf_centroid_x write(6,*) 'AHFinder_gau: ahf_centroid_y = ', . ahf_centroid_y write(6,*) 'AHFinder_gau: ahf_centroid_z = ', . ahf_centroid_z end if end if ! 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 deallocate(rr) deallocate(xa,ya,za) deallocate(txx,tyy,tzz,txy,txz,tyz) deallocate(g11,g22,g12) return end if ! *********************** ! *** INTERPOLATE *** ! *********************** ! Here I should really only interpolate the gaussian curvature. ! The interpolation of the metric is only needed for the old ! gaussian curvature calculation (which I dont use any more), ! and for the calculation of the circumferences (which really ! should be done in the routine AHFinder_int.F). I need to fix ! this sometime soon! (October 2000). ! parameter, local interpolator, and coordinate system handle. param_table_handle = -1 interp_handle = -1 coord_system_handle = -1 options_string = "order = " // char(ichar('0') + interpolation_order) call Util_TableCreateFromString (param_table_handle, options_string) if (param_table_handle .lt. 0) then call CCTK_WARN(0,"Cannot create parameter table for interpolator") endif call CCTK_FortranString (nchars, interpolation_operator, operator) call CCTK_InterpHandle (interp_handle, operator) if (interp_handle .lt. 0) then call CCTK_WARN(0,"Cannot get handle for interpolation ! Forgot to activate an implementation providing interpolation operators ??") endif call CCTK_CoordSystemHandle (coord_system_handle, "cart3d") if (coord_system_handle .lt. 0) then call CCTK_WARN(0,"Cannot get handle for cart3d coordinate system ! Forgot to activate an implementation providing coordinates ??") endif ! fill in the input/output arrays for the interpolator interp_coords(1) = CCTK_PointerTo(xa) interp_coords(2) = CCTK_PointerTo(ya) interp_coords(3) = CCTK_PointerTo(za) call CCTK_VarIndex (vindex, "admbase::gxx") in_array_indices(1) = vindex call CCTK_VarIndex (vindex, "admbase::gyy") in_array_indices(2) = vindex call CCTK_VarIndex (vindex, "admbase::gzz") in_array_indices(3) = vindex call CCTK_VarIndex (vindex, "admbase::gxy") in_array_indices(4) = vindex call CCTK_VarIndex (vindex, "admbase::gxz") in_array_indices(5) = vindex call CCTK_VarIndex (vindex, "admbase::gyz") in_array_indices(6) = vindex call CCTK_VarIndex (vindex, "ahfinder::ahfgauss") in_array_indices(7) = vindex out_arrays(1) = CCTK_PointerTo(txx) out_arrays(2) = CCTK_PointerTo(tyy) out_arrays(3) = CCTK_PointerTo(tzz) out_arrays(4) = CCTK_PointerTo(txy) out_arrays(5) = CCTK_PointerTo(txz) out_arrays(6) = CCTK_PointerTo(tyz) out_arrays(7) = CCTK_PointerTo(gaussian) out_array_type_codes = CCTK_VARIABLE_REAL ! Interpolation. call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle, . param_table_handle, coord_system_handle, . npoints, CCTK_VARIABLE_REAL, interp_coords, . 7, in_array_indices, . 7, out_array_type_codes, out_arrays) if (ierror < 0) then call CCTK_WARN (1, "interpolator call returned an error code"); endif ! release parameter table call Util_TableDestroy (ierror, param_table_handle) ! ********************************* ! *** FIND SPHERICAL METRIC *** ! ********************************* if (myproc.eq.0) then do j=1,nphi+1 do i=1,ntheta+1 ! Find {theta,phi}. theta = dtheta*dble(i-1) phi = dphi*dble(j-1) + phistart ! 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.1).and.(i.ne.ntheta+1)) then ft = idtheta*(rr(i+1,j) - rr(i-1,j)) else if (i.eq.1) then ft = - idtheta*(three*rr(1,j) . - four*rr(2,j) + rr(3,j)) else ft = + idtheta*(three*rr(ntheta,j) . - four*rr(ntheta,j) + rr(ntheta-1,j)) end if if (nonaxi) then if ((j.ne.1).and.(j.ne.nphi+1)) then fp = idphi*(rr(i,j+1) - rr(i,j-1)) else if (j.eq.1) then fp = - idphi*(three*rr(i,1) . - four*rr(i,2) + rr(i,3)) else fp = + idphi*(three*rr(i,nphi) . - four*rr(i,nphi) + rr(i,nphi-1)) 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 end if ! ************************************ ! *** CALCULATE CIRCUMFERENCES *** ! ************************************ if (myproc.eq.0) then ! Equatorial circumference. if (refz) then do j=1,nphi circ_eq = circ_eq + dphi . *sqrt(abs(half*(g22(ntheta+1,j) + g22(ntheta+1,j+1)))) end do else aux = half*pi/dtheta i = int(aux) aux = aux - int(aux) if (cartoon) then circ_eq = 6.283185307D0*sqrt(g22(i,0)) else do j=1,nphi circ_eq = circ_eq + dphi . *((one-aux)*sqrt(abs(half*(g22(i,j) + g22(i,j+1)))) . + aux*sqrt(abs(half*(g22(i+1,j) + g22(i+1,j+1))))) end do end if end if ! Meridians. do i=1,ntheta meri_p1 = meri_p1 + dtheta . *sqrt(abs(half*(g11(i,1) + g11(i+1,1)))) end do if (cartoon) then meri_p2 = meri_p1 else if (refx.and.refy) then do i=1,ntheta meri_p2 = meri_p2 + dtheta . *sqrt(abs(half*(g11(i,nphi+1) + g11(i+1,nphi+1)))) end do else aux = half*pi/dphi j = int(aux) aux = aux - int(aux) do i=1,ntheta meri_p2 = meri_p2 + dtheta . *((one-aux)*sqrt(abs(half*(g11(i,j) + g11(i,j)))) . + aux*sqrt(abs(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 is 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 (we would need higher ! order interpolation). I did not want to delete this code ! since it might be useful in the future, so here 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.0) then ! Interior points. do j=2,nphi do i=2,ntheta ! 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-1) 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(1,:) = two*gaussian(2,:) - gaussian(3,:) gaussian(ntheta+1,:) = two*gaussian(ntheta,:) . - gaussian(ntheta-1,:) gaussian(:,1) = two*gaussian(:,2) - gaussian(:,3) gaussian(:,nphi+1) = two*gaussian(:,nphi) . - gaussian(:,nphi-1) end if end if ! ************************************ ! *** INITIALIZE FIRSTGAU FLAG *** ! ************************************ if (find3) then firstgau = firstcal(mfind) firstcal(mfind) = .false. else firstgau = firstcal(4) firstcal(4) = .false. end if ! *********************************** ! *** SAVE GAUSSIAN CURVATURE *** ! *********************************** if (myproc.eq.0) then if (find3) then if (mfind.eq.1) then gaussf = filestr(1:nfile)//"/ahf_0.gauss" else if (mfind.eq.2) then gaussf = filestr(1:nfile)//"/ahf_1.gauss" else gaussf = filestr(1:nfile)//"/ahf_2.gauss" end if else gaussf = filestr(1:nfile)//"/ahf.gauss" end if inquire(file=gaussf, exist=gaussf_exists) if (firstgau) then if (find3) then if (mfind.eq.3) then firstgau = .false. end if else firstgau = .false. end if end if ! On first call replace file. if (.not.gaussf_exists) then call AHFinder_WriteGaussHeader(CCTK_PASS_FTOF, gaussf) ! On subsequent calls, append to file. else open(11,file=gaussf,form='formatted', status='old', . position='append') write(11,*) write(11,"(A11,I4)") '# Time step',cctk_iteration write(11,"(A6,ES11.3)") '# Time',cctk_time write(11,"(A6,I4)") '# Call',ahf_ncall end if ! Save data points. do i=1,ntheta+1 do j=1,nphi+1 write(11,"(ES11.3)") gaussian(i,j) end do end do ! Close file. close(11) end if ! ***************************** ! *** DEALLOCATE MEMORY *** ! ***************************** deallocate(rr) deallocate(xa,ya,za) deallocate(txx,tyy,tzz,txy,txz,tyz) deallocate(g11,g22,g12) ! *************** ! *** END *** ! *************** end subroutine AHFinder_gau ! ****************************************** ! *** String comparison utility function *** ! ****************************************** function str_comp(a1,a2) logical :: str_comp character(len=*),intent(in) :: a1, a2 if ( ( verify(trim(a1),trim(a2)).eq.0 ) .and. . ( len_trim(a1) .eq. len_trim(a2) ) ) then str_comp = .true. else str_comp = .false. end if end function str_comp