/*@@ @file AHFinder_int.F @date April 1998 @author Miguel Alcubierre @desc Find surface integrals. The integrals are done in parallel, but the number of processors is assumed to be either the square of an integer or a power of two (if this is not the case, I just use less processors). @enddesc @version $Header$ @@*/ #include "cctk.h" #include "cctk_Parameters.h" #include "cctk_Arguments.h" #include "cctk_Functions.h" subroutine AHFinder_int(CCTK_ARGUMENTS) use AHFinder_dat implicit none DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS integer i,j,l,m integer npt,npp integer l_ntheta,l_nphi,theta0,phi0 integer npoints integer auxi,vindex integer ierror integer interp_handle,coord_system_handle integer reduction_handle,param_table_handle character(30) options_string character(128) operator CCTK_INT red_tmp, nchars CCTK_POINTER, dimension(3) :: interp_coords CCTK_INT, dimension(9) :: in_array_indices CCTK_POINTER, dimension(9) :: out_arrays CCTK_INT, dimension(9) :: out_array_type_codes CCTK_REAL LEGEN CCTK_REAL xp,yp,zp,rp CCTK_REAL xw,yw,zw CCTK_REAL cost,sint,cosp,sinp CCTK_REAL trr,ttt,tpp,trt,trp,ttp,ft,fp CCTK_REAL dtheta,dphi,dtp,idtheta,idphi,phistart CCTK_REAL det,idet CCTK_REAL zero,quarter,half,one,two,three,four,pi CCTK_REAL sina,cosa,intw,grad,sigma,h0 CCTK_REAL theta,phi CCTK_REAL aux,aux1,aux2 CCTK_REAL tempv(0:lmax),tempm(lmax,lmax) CCTK_REAL, dimension(3,3) :: gupij CCTK_REAL, allocatable, dimension(:) :: costheta,sintheta CCTK_REAL, allocatable, dimension(:) :: cosphi,sinphi CCTK_REAL, allocatable, dimension(:,:) :: rr,xa,ya,za,da,exp,gradn CCTK_REAL, allocatable, dimension(:,:) :: txx,tyy,tzz,txy,txz,tyz CCTK_REAL, allocatable, dimension(:,:) :: intmask ! Reduction variables CCTK_POINTER_TO_CONST, dimension(1) :: input_array CCTK_POINTER, dimension(1) :: reduction_value CCTK_INT, dimension(1) ::input_array_type ! Variables to be saved for next call. save dtheta,dphi,dtp,idtheta,idphi,phistart save npt,npp,l_ntheta,l_nphi,theta0,phi0 ! Description of variables: ! ! i,j,l,m Counters. ! ! theta Latitude. ! phi Longitude. ! ! npt number of processors in theta direction. ! npp number of processors in phi direction. ! ! l_ntheta Local number of grid points in theta direction. ! l_nphi Local number of grid points in phi direction. ! ! theta0 Local origin in theta direction. ! phi0 Local origin in phi direction. ! ! npoints Number of points to interpolate . ! ! 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. ! ! exp Interpolated expansion. ! ! txx Array with interpolated gxx metric component. ! tyy Array with interpolated gyy metric component. ! tzz Array with interpolated gzz metric component. ! txy Array with interpolated gxy metric component. ! txz Array with interpolated gxz metric component. ! tyz Array with interpolated gyz metric component. ! ! gradn Array with interpolated norm of gradient of horizon function. ! ! intmask Array with interpolated mask. ! ! cost cos(theta) ! sint sin(theta) ! cosp cos(phi) ! sinp sin(phi) ! ! trr Metric component {r,r}. ! ttt Metric component {theta,theta}. ! tpp Metric component {phi,phi}. ! trt Metric component {r,theta}. ! trp Metric component {r,phi}. ! ttp Metric component {theta,phi}. ! ! ft dr/dtheta ! fp dr/dphi ! ! dtheta Grid spacing in theta. ! dphi Grid spacing in phi. ! dtp dtheta*dphi. ! ! idtheta 1/(2 dtheta) ! idphi 1/(2 dphi) ! ! phistart Origin for phi (normally 0, but for some symmetries -pi/2). ! ! da Area element. ! ******************* ! *** NUMBERS *** ! ******************* zero = 0.0D0 quarter = 0.25D0 half = 0.5D0 one = 1.0D0 two = 2.0D0 three = 3.0D0 four = 4.0D0 pi = acos(-one) ! ***************************************************************** ! *** TOTAL NUMBER OF PROCESSORS AND LOCAL PROCESSOR NUMBER *** ! ***************************************************************** myproc = CCTK_MyProc(cctkGH) nprocs = CCTK_nProcs(cctkGH) ! **************************************** ! *** GET REDUCTION HANDLE FOR SUM *** ! **************************************** call CCTK_LocalArrayReductionHandle(reduction_handle,"sum") if (reduction_handle.lt.0) then call CCTK_WARN(1,"Cannot get handle for sum reduction ! Forgot to activate an implementation providing reduction operators ??") end if ! ****************************************************** ! *** THINGS TO DO ON FIRST CALL TO THIS ROUTINE *** ! ****************************************************** if (firstint) then firstint = .false. ! ****************************************** ! *** 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 = quarter*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 ! ********************************************* ! *** DIVIDE 2D DOMAIN AMONG PROCESSORS *** ! ********************************************* ! For cartoon this is trivial. if (cartoon) then npt = nprocs npp = 1 ! Otherwise it is more complicated. At the moment I ! only consider numbers of processors that are either ! the square of an integer, or a power of two. ! If neither of this is the case, then I use fewer ! processors to make sure that I always deal with ! either a square or a power of two. This can ! certainly be improved, but it might not be that ! urgent since in most cases the number of processors ! will be a power of two. else aux = sqrt(dble(nprocs)) ! One processor. if (nprocs.eq.1) then npt = 1 npp = 1 ! Square of an integer. else if (aux.eq.int(aux)) then npt = int(aux) npp = int(aux) ! The number of processors is not the square of an integer. ! We will now check if it is a power of two. else l = 1 m = 1 i = 1 do while (l*m.lt.nprocs) npt = l npp = m if (mod(i,2).ne.0) then m = 2*m else l = 2*l end if i = i + 1 end do ! "nprocs" is a power of two. if (l*m.eq.nprocs) then npt = l npp = m ! "nprocs" is not a power of two. This is where I ! would need to do something clever. At the moment, ! I will just use fewer processors. The number of ! processors I use is the power of two or the square ! of an integer that is closest to "nprocs". else if (npt*npp.lt.int(aux)**2) then npt = int(aux) npp = int(aux) end if end if end if end if ! Now be careful not to use too many procs and end up ! with very few points per proc. if (ntheta/npt.lt.3) npt = ntheta/3 if (.not.cartoon) then if (nphi/npp.lt.3) npp = nphi/3 end if ! Figure out the number of grid points per processor, ! and the "origin" for a given processor. if (myproc.ge.npt*npp) then l_ntheta = 0 l_nphi = 0 theta0 = 0 phi0 = 0 else ! First take the number of grid points per processor ! in a given direction to be equal to the total number ! of grid points divided by the number of processors ! in that direction. l_ntheta = ntheta/npt l_nphi = nphi/npp ! And take the "origin" as the corresponding displacement ! from the real origin. Notice that mod(myproc,npt) ! tells my on which "theta" bin I am, while myproc/npt ! tells me on which "phi" bin I am. i = mod(myproc,npt) j = myproc/npt theta0 = i*l_ntheta phi0 = j*l_nphi ! Now correct all this if the total number of grid points ! in a given direction was not exactly divisible by the ! number of processors on that direction. if (l_ntheta*npt.ne.ntheta) then ! Find residue on theta direction. l = ntheta - l_ntheta*npt ! Distribute residue in first few processors, and ! displace theta0 accordingly. if (i.lt.l) then l_ntheta = l_ntheta + 1 theta0 = theta0 + i else theta0 = theta0 + l end if end if if (l_nphi*npp.ne.nphi) then ! Find residue on phi direction. l = nphi - l_nphi*npp ! Distribute residue in first few processors, and ! displace phi0 accordingly. if (j.lt.l) then l_nphi = l_nphi + 1 phi0 = phi0 + j else phi0 = phi0 + l end if end if end if end if ! ************************************** ! *** ALLOCATE MEMORY FOR ARRAYS *** ! ************************************** allocate(costheta(1:l_ntheta+1)) allocate(sintheta(1:l_ntheta+1)) allocate(cosphi(1:l_nphi+1),sinphi(1:l_nphi+1)) allocate(rr(1:l_ntheta+1,1:l_nphi+1)) allocate(xa(1:l_ntheta+1,1:l_nphi+1)) allocate(ya(1:l_ntheta+1,1:l_nphi+1)) allocate(za(1:l_ntheta+1,1:l_nphi+1)) allocate(da(1:l_ntheta+1,1:l_nphi+1)) allocate(exp(1:l_ntheta+1,1:l_nphi+1)) allocate(gradn(1:l_ntheta+1,1:l_nphi+1)) allocate(txx(1:l_ntheta+1,1:l_nphi+1)) allocate(tyy(1:l_ntheta+1,1:l_nphi+1)) allocate(tzz(1:l_ntheta+1,1:l_nphi+1)) allocate(txy(1:l_ntheta+1,1:l_nphi+1)) allocate(txz(1:l_ntheta+1,1:l_nphi+1)) allocate(tyz(1:l_ntheta+1,1:l_nphi+1)) allocate(intmask(1:l_ntheta+1,1:l_nphi+1)) ! ******************************** ! *** INITIALIZE VARIABLES *** ! ******************************** ! Initialize error flags. interror1 = 0 interror2 = 0 interror3 = 0 ! Initialize surface integrals. intexp = zero intexp2 = zero intarea = zero intexpdel2 = zero ! For flow algorithm, initialize spectral components. if (flow) then hflow0 = zero hflowc = zero hflows = zero cflow0 = zero cflowc = zero cflows = zero nflow0 = zero nflowc = zero nflows = zero end if ! ************************************* ! *** FIND INTERPOLATING POINTS *** ! ************************************* if (myproc.ge.npt*npp) then xa = half*(xmx+xmn) ya = half*(ymx+ymn) za = half*(zmx+zmn) else do j=1,l_nphi+1 ! Find phi. phi = dphi*dble(j-1+phi0) + phistart ! Find sines and cosines of phi. cosp = cos(phi) sinp = sin(phi) cosphi(j) = cosp sinphi(j) = sinp do i=1,l_ntheta+1 ! Find theta. theta = dtheta*dble(i-1+theta0) ! Find sines and cosines of theta. cost = cos(theta) sint = sin(theta) costheta(i) = cost sintheta(i) = sint ! 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 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 interror1 = 1 end if ! Find cartesian coordinates. aux = rp*sint xp = aux*cosp + xc yp = aux*sinp + yc zp = rp*cost + zc ! Save arrays. rr(i,j) = rp xa(i,j) = xp ya(i,j) = yp za(i,j) = zp ! Check if we are within bounds. 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 interror2 = 1 end if end do end do end if ! ******************************************* ! *** REDUCE ERRORS ACROSS PROCESSORS *** ! ******************************************* input_array_type(1) = CCTK_VARIABLE_INT reduction_value(1)= CCTK_PointerTo(red_tmp) input_array(1) = CCTK_PointerTo(interror1) call CCTK_ReduceArraysGlobally(ierror, cctkGH, -1,reduction_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 of norm failed!") end if interror1 = red_tmp input_array_type(1) = CCTK_VARIABLE_INT reduction_value(1)= CCTK_PointerTo(red_tmp) input_array(1) = CCTK_PointerTo(interror2) call CCTK_ReduceArraysGlobally(ierror, cctkGH, -1,reduction_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 of norm failed!") end if interror2 = red_tmp ! If there was an error on any processor then assign ! large values to the integrals and return from the ! subroutine (but remember to deallocate the arrays ! first!). interror = .false. if ((interror1.ne.0).or.(interror2.ne.0)) then intexp = 1.0D10 intexp2 = 1.0D10 intarea = 1.0D10 intexpdel2 = 1.0D10 interror = .true. goto 10 end if ! *********************** ! *** INTERPOLATE *** ! *********************** ! Number of points to interpolate. npoints = (l_ntheta+1)*(l_nphi+1) ! 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::ahf_exp") in_array_indices(7) = vindex call CCTK_VarIndex (vindex, "ahfinder::ahmask") in_array_indices(8) = vindex call CCTK_VarIndex (vindex, "ahfinder::ahfgradn") in_array_indices(9) = 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(exp) out_arrays(8) = CCTK_PointerTo(intmask) out_arrays(9) = CCTK_PointerTo(gradn) out_array_type_codes = CCTK_VARIABLE_REAL ! For minimization we need the interpolated metric, ! the expansion and the mask. if (.not.flow) then call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle, . param_table_handle, coord_system_handle, . npoints, CCTK_VARIABLE_REAL, interp_coords, . 8, in_array_indices, . 8, out_array_type_codes, out_arrays) if (ierror < 0) then call CCTK_WARN (1, "interpolator call returned an error code"); endif ! For N flow, we need the interpolated metric, the expansion, ! the norm of the gradient of the horizon function, and the mask. else if (nw.ne.zero) then call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle, . param_table_handle, coord_system_handle, . npoints, CCTK_VARIABLE_REAL, interp_coords, . 9, in_array_indices, . 9, out_array_type_codes, out_arrays) if (ierror < 0) then call CCTK_WARN (1, "interpolator call returned an error code"); endif ! For H or C flows, we need the interpolated expansion, the ! norm of the gradient of the horizon function, and the mask. else call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle, . param_table_handle, coord_system_handle, . npoints, CCTK_VARIABLE_REAL, interp_coords, . 3, in_array_indices(7), . 3, out_array_type_codes(7), out_arrays(7)) if (ierror < 0) then call CCTK_WARN (1, "interpolator call returned an error code"); endif end if ! release parameter table call Util_TableDestroy (ierror, param_table_handle) ! *************************************** ! *** CHECK IF WE ARE INSIDE MASK *** ! *************************************** ! Check if we are either inside mask, or to close for comfort. ! The way I do this is very simple, I just check the interpolated ! values of the mask. If any of these values is not 1.0, it means ! that a point inside the mask is contaminating the interpolation, ! and this means in turn that we are too close to the mask. ! Notice that in practice I do not ask for the interpolated values ! to be exactly 1.0, since round off error can modify the value ! slightly without really meaning that we are close to the mask. if (.not.CCTK_EQUALS(ahf_mask,'off')) then if (myproc.lt.npt*npp) then do j=1,l_nphi+1 do i=1,l_ntheta+1 if (intmask(i,j).lt.0.99D0) then interror3 = 1 end if end do end do end if input_array_type(1) = CCTK_VARIABLE_INT reduction_value(1)= CCTK_PointerTo(red_tmp) input_array(1) = CCTK_PointerTo(interror3) call CCTK_ReduceArraysGlobally(ierror, cctkGH, -1,reduction_handle, -1, . 1, input_array,0, 0, input_array_type, 1, . input_array_type, reduction_value) interror3 = red_tmp if (interror3.ne.0) then intexp = 1.0D10 intexp2 = 1.0D10 intarea = 1.0D10 intexpdel2 = 1.0D10 interror = .true. goto 10 end if end if ! ************************* ! *** AREA ELEMENTS *** ! ************************* ! Find 2D array for area elements. Notice that ! this is not needed for the flow algorithm! if (.not.flow) then if (myproc.ge.npt*npp) then da = zero else do j=1,l_nphi+1 ! Find sines and cosines of phi. cosp = cosphi(j) sinp = sinphi(j) do i=1,l_ntheta+1 ! Find {rp}. rp = rr(i,j) ! Find sines and cosines of theta. cost = costheta(i) sint = sintheta(i) ! Find spherical metric components. 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. Notice that ! since this is done even in inter-processor boundaries, ! the final result will depend on the number of processors, ! but the differences will converge away at high order. if ((i.ne.1).and.(i.ne.l_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(l_ntheta+1,j) . - four*rr(l_ntheta,j) + rr(l_ntheta-1,j)) end if if (nonaxi) then if ((j.ne.1).and.(j.ne.l_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,l_nphi+1) . - four*rr(i,l_nphi) + rr(i,l_nphi-1)) end if else fp = zero end if ! Find area element. For this I use the fact that if ! we have r = f(theta,phi), then the area element is: ! ! / 2 ! da = | [ gtt + grr (f,t) + 2 grt f,t ] ! \ ! 2 ! [ gpp + grr (f,p) + 2 grp f,p ] - [ gtp ! ! 2 \ 1/2 ! + grr f,t f,p + grt f,p + grp f,t ] | dtheta dphi ! / ! ! where f,t = df/dtheta and f,p = df/dphi. aux = (ttt + trr*ft**2 + two*trt*ft) . * (tpp + trr*fp**2 + two*trp*fp) . - (ttp + ft*(trr*fp + trp) + trt*fp)**2 if (aux.gt.zero) then aux = sqrt(aux) else aux = zero end if ! Multiply with dtheta*dphi. da(i,j) = aux*dtp end do end do end if end if ! ******************************************* ! *** AREA AND INTEGRALS OF EXPANSION *** ! ******************************************* ! To find the integrals I use a second order method: ! I approximate the value of the integrand at each ! cell centre by an average of the values at the ! four cell corners, and then add up all cells. ! Notice that for the flow algorithm I do not need to ! calculate the area nor the integrals of the expansion ! except at the very end, but by then the flag 'flow' has ! been switched off. inside_min_count = 0 inside_min_neg_count = 0 if (.not.flow) then if (myproc.lt.npt*npp) then do j=1,l_nphi do i=1,l_ntheta intarea = intarea + quarter . *(da(i,j ) + da(i+1,j ) . + da(i,j+1) + da(i+1,j+1)) intexp = intexp + quarter . *(da(i ,j )*exp(i ,j ) . + da(i+1,j )*exp(i+1,j ) . + da(i ,j+1)*exp(i ,j+1) . + da(i+1,j+1)*exp(i+1,j+1)) intexp2 = intexp2 + quarter . *(da(i ,j )*exp(i ,j )**2 . + da(i+1,j )*exp(i+1,j )**2 . + da(i ,j+1)*exp(i ,j+1)**2 . + da(i+1,j+1)*exp(i+1,j+1)**2) intexpdel2 = intexpdel2 + quarter . *(da(i ,j )*(exp(i ,j )+trapped_surface_delta)**2 . + da(i+1,j )*(exp(i+1,j )+trapped_surface_delta)**2 . + da(i ,j+1)*(exp(i ,j+1)+trapped_surface_delta)**2 . + da(i+1,j+1)*(exp(i+1,j+1)+trapped_surface_delta)**2) inside_min_count = inside_min_count + 1 if ((exp(i ,j ).lt.0.0D0).and. . (exp(i+1,j ).lt.0.0D0).and. . (exp(i ,j+1).lt.0.0D0).and. . (exp(i+1,j+1).lt.0.0D0)) then inside_min_neg_count = inside_min_neg_count + 1 end if end do end do end if ! Reduce integrals. input_array_type(1) = CCTK_VARIABLE_REAL reduction_value(1)= CCTK_PointerTo(aux) input_array(1) = CCTK_PointerTo(intarea) call CCTK_ReduceArraysGlobally(ierror, cctkGH, -1,reduction_handle, -1, . 1, input_array,0, 0, input_array_type, 1, . input_array_type, reduction_value) intarea=aux input_array_type(1) = CCTK_VARIABLE_REAL reduction_value(1)= CCTK_PointerTo(aux) input_array(1) = CCTK_PointerTo(intexp) call CCTK_ReduceArraysGlobally(ierror, cctkGH, -1,reduction_handle, -1, . 1, input_array,0, 0, input_array_type, 1, . input_array_type, reduction_value) intexp=aux input_array_type(1) = CCTK_VARIABLE_REAL reduction_value(1)= CCTK_PointerTo(aux) input_array(1) = CCTK_PointerTo(intexp2) call CCTK_ReduceArraysGlobally(ierror, cctkGH, -1,reduction_handle, -1, . 1, input_array,0, 0, input_array_type, 1, . input_array_type, reduction_value) intexp2=aux input_array_type(1) = CCTK_VARIABLE_REAL reduction_value(1)= CCTK_PointerTo(aux) input_array(1) = CCTK_PointerTo(intexpdel2) call CCTK_ReduceArraysGlobally(ierror, cctkGH, -1,reduction_handle, -1, . 1, input_array,0, 0, input_array_type, 1, . input_array_type, reduction_value) intexpdel2=aux input_array_type(1) = CCTK_VARIABLE_INT reduction_value(1)= CCTK_PointerTo(red_tmp) input_array(1) = CCTK_PointerTo(inside_min_neg_count) call CCTK_ReduceArraysGlobally(ierror, cctkGH, -1,reduction_handle, -1, . 1, input_array,0, 0, input_array_type, 1, . input_array_type, reduction_value) inside_min_neg_count=red_tmp input_array_type(1) = CCTK_VARIABLE_INT reduction_value(1)= CCTK_PointerTo(red_tmp) input_array(1) = CCTK_PointerTo(inside_min_count) call CCTK_ReduceArraysGlobally(ierror, cctkGH, -1,reduction_handle, -1, . 1, input_array,0, 0, input_array_type, 1, . input_array_type, reduction_value) inside_min_count=red_tmp end if ! ******************************* ! *** SPECTRAL COMPONENTS *** ! ******************************* ! For flow algorithm, find spectral components. ! For this I need to recalculate the Legendre polynomials, ! but I see no way around it. Also, to avoid having to ! define lots of 2D arrays I do the integrals inside ! the loop. This means that I have to take care to see ! if I am on an edge or corner point. Notice also that ! these integrals do not involve the area element. if (flow) then if (myproc.lt.npt*npp) then ! Find h0. if (find_trapped_surface) then h0 = - trapped_surface_delta else h0 = zero end if ! Find integrals. do j=1,l_nphi+1 ! Find phi. phi = dphi*dble(j+phi0) + phistart do i=1,l_ntheta+1 ! Find sines and cosines of theta. cost = costheta(i) sint = sintheta(i) ! Find norm of gradient of horizon function. grad = gradn(i,j) ! Find weight factor. intw = dtp*(exp(i,j)-h0)*sint**2 ! Modify weight factor at corners and edges. if (((j.eq.1).or.(j.eq.l_nphi+1)).and. . ((i.eq.1).or.(i.eq.l_ntheta+1))) then intw = quarter*intw else if (((j.eq.1).or.(j.eq.l_nphi+1)).or. . ((i.eq.1).or.(i.eq.l_ntheta+1))) then intw = half*intw end if ! Find "sigma" for Nakamura flow (see Carstens paper). if (nw.eq.0.0) then sigma = zero else xw = xa(i,j) - xc yw = ya(i,j) - yc zw = za(i,j) - zc ! Calculate metric on surface. det = txx(i,j)*tyy(i,j)*tzz(i,j) . + two*txy(i,j)*txz(i,j)*tyz(i,j) . - txx(i,j)*tyz(i,j)**2 - tyy(i,j)*txz(i,j)**2 . - tzz(i,j)*txy(i,j)**2 idet = one/det gupij(1,1) = idet*(tyy(i,j)*tzz(i,j)-tyz(i,j)**2) gupij(2,2) = idet*(txx(i,j)*tzz(i,j)-txz(i,j)**2) gupij(3,3) = idet*(txx(i,j)*tyy(i,j)-txy(i,j)**2) gupij(1,2) = idet*(txz(i,j)*tyz(i,j)-tzz(i,j)*txy(i,j)) gupij(2,1) = gupij(1,2) gupij(1,3) = idet*(txy(i,j)*tyz(i,j)-tyy(i,j)*txz(i,j)) gupij(3,1) = gupij(1,3) gupij(2,3) = idet*(txy(i,j)*txz(i,j)-txx(i,j)*tyz(i,j)) gupij(3,2) = gupij(2,3) call AHFinder_calcsigma(CCTK_ARGUMENTS, . xw,yw,zw,gupij,sigma) end if ! Monopole components. aux = intw*grad hflow0(0) = hflow0(0) + intw cflow0(0) = cflow0(0) + aux nflow0(0) = nflow0(0) + aux*sigma ! Axisymmetric components. do l=1+stepz,lmax,1+stepz aux = intw*LEGEN(l,0,cost) aux1 = aux*grad hflow0(l) = hflow0(l) + aux cflow0(l) = cflow0(l) + aux1 nflow0(l) = nflow0(l) + aux1*sigma end do ! Non-axisymmetric components. if (nonaxi) then do m=1,lmax aux = dble(m)*phi sina = sin(aux) cosa = cos(aux) do l=m,lmax aux = intw*LEGEN(l,m,cost) aux1 = aux*cosa aux2 = aux1*grad hflowc(l,m) = hflowc(l,m) + aux1 cflowc(l,m) = cflowc(l,m) + aux2 nflowc(l,m) = nflowc(l,m) + aux2*sigma if (.not.refy) then aux1 = aux*sina aux2 = aux1*grad hflows(l,m) = hflows(l,m) + aux1 cflows(l,m) = cflows(l,m) + aux2 nflows(l,m) = nflows(l,m) + aux2*sigma end if end do end do end if end do end do end if ! Reduce integrals. ! Axisymmetric components. auxi = lmax+1 if (hw.ne.zero) then input_array_type(1) = CCTK_VARIABLE_REAL reduction_value(1)= CCTK_PointerTo(tempv) input_array(1) = CCTK_PointerTo(hflow0) call CCTK_ReduceArraysGlobally(ierror, cctkGH, -1,reduction_handle, -1, . 1, input_array,1, auxi, input_array_type, 1, . input_array_type, reduction_value) hflow0 = tempv end if if (cw.ne.zero) then input_array_type(1) = CCTK_VARIABLE_REAL reduction_value(1)= CCTK_PointerTo(tempv) input_array(1) = CCTK_PointerTo(cflow0) call CCTK_ReduceArraysGlobally(ierror, cctkGH, -1,reduction_handle, -1, . 1, input_array,1, auxi, input_array_type, 1, . input_array_type, reduction_value) cflow0 = tempv end if if (nw.ne.zero) then input_array_type(1) = CCTK_VARIABLE_REAL reduction_value(1)= CCTK_PointerTo(tempv) input_array(1) = CCTK_PointerTo(nflow0) call CCTK_ReduceArraysGlobally(ierror, cctkGH, -1,reduction_handle, -1, . 1, input_array,1, auxi, input_array_type, 1, . input_array_type, reduction_value) nflow0 = tempv end if ! Non-axisymmetric components. if (nonaxi) then auxi = lmax*lmax if (hw.ne.zero) then input_array_type(1) = CCTK_VARIABLE_REAL reduction_value(1)= CCTK_PointerTo(tempm) input_array(1) = CCTK_PointerTo(hflowc) call CCTK_ReduceArraysGlobally(ierror, cctkGH, -1,reduction_handle, -1, . 1, input_array,1, auxi, input_array_type, 1, . input_array_type, reduction_value) hflowc = tempm end if if (cw.ne.zero) then input_array_type(1) = CCTK_VARIABLE_REAL reduction_value(1)= CCTK_PointerTo(tempm) input_array(1) = CCTK_PointerTo(cflowc) call CCTK_ReduceArraysGlobally(ierror, cctkGH, -1,reduction_handle, -1, . 1, input_array,1, auxi, input_array_type, 1, . input_array_type, reduction_value) cflowc = tempm end if if (nw.ne.zero) then input_array_type(1) = CCTK_VARIABLE_REAL reduction_value(1)= CCTK_PointerTo(tempm) input_array(1) = CCTK_PointerTo(nflowc) call CCTK_ReduceArraysGlobally(ierror, cctkGH, -1,reduction_handle, -1, . 1, input_array,1, auxi, input_array_type, 1, . input_array_type, reduction_value) nflowc = tempm end if if (.not.refy) then if (hw.ne.zero) then input_array_type(1) = CCTK_VARIABLE_REAL reduction_value(1)= CCTK_PointerTo(tempm) input_array(1) = CCTK_PointerTo(hflows) call CCTK_ReduceArraysGlobally(ierror, cctkGH, -1,reduction_handle, -1, . 1, input_array,1, auxi, input_array_type, 1, . input_array_type, reduction_value) hflows = tempm end if if (cw.ne.zero) then input_array_type(1) = CCTK_VARIABLE_REAL reduction_value(1)= CCTK_PointerTo(tempm) input_array(1) = CCTK_PointerTo(cflows) call CCTK_ReduceArraysGlobally(ierror, cctkGH, -1,reduction_handle, -1, . 1, input_array,1, auxi, input_array_type, 1, . input_array_type, reduction_value) cflows = tempm end if if (nw.ne.zero) then input_array_type(1) = CCTK_VARIABLE_REAL reduction_value(1)= CCTK_PointerTo(tempm) input_array(1) = CCTK_PointerTo(nflows) call CCTK_ReduceArraysGlobally(ierror, cctkGH, -1,reduction_handle, -1, . 1, input_array,1, auxi, input_array_type, 1, . input_array_type, reduction_value) nflows = tempm end if end if end if end if ! ***************************************************** ! *** RESCALE INTEGRALS ACCORDING TO SYMMETRIES *** ! ***************************************************** ! For cartoon multiply the integrals with 2*pi. if (cartoon) then intw = 6.283185307D0 if (.not.flow) then intexp = intw*intexp intexp2 = intw*intexp2 intarea = intw*intarea intexpdel2 = intw*intexpdel2 end if if (flow) then hflow0 = intw*hflow0 cflow0 = intw*cflow0 nflow0 = intw*nflow0 end if ! Other symmetries. else if (refx.or.refy.or.refz) then intw = one if (refx) intw = two*intw if (refy) intw = two*intw if (refz) intw = two*intw if (.not.flow) then intexp = intw*intexp intexp2 = intw*intexp2 intarea = intw*intarea intexpdel2 = intw*intexpdel2 end if ! Spectral components for flow. if (flow) then ! Axisymmetric. hflow0 = intw*hflow0 cflow0 = intw*cflow0 nflow0 = intw*nflow0 ! Non-axisymmetric. if (nonaxi) then hflowc = intw*hflowc cflowc = intw*cflowc nflowc = intw*nflowc if (.not.refy) then hflows = intw*hflows cflows = intw*cflows nflows = intw*nflows end if end if end if end if ! ***************************** ! *** DEALLOCATE MEMORY *** ! ***************************** 10 continue deallocate(costheta,sintheta) deallocate(cosphi,sinphi) deallocate(rr,xa,ya,za) deallocate(da,exp,gradn) deallocate(txx,tyy,tzz,txy,txz,tyz) deallocate(intmask) ! *************** ! *** END *** ! *************** end subroutine AHFinder_int