/*@@ @file num_init.F @date October 1999 @author Miguel Alcubierre @desc Initial position and covariant momentum of the geodesics. @enddesc @@*/ #include "cctk.h" #include "cctk_Parameters.h" #include "cctk_Arguments.h" #include "cctk_Functions.h" /*@@ @routine num_init @date August 1997 @author Miguel Alcubierre @desc Gives initial position and covariant momentum of the geodesics, and initializes the proper time array. The initial positions of the geodesics are determined by the values of the parameter "geopos" (see file timegeodesic_param.h). I also choose the initial geodesics to point in the normal direction to the hypersurfaces. @enddesc @calls @calledby @history @endhistory @@*/ subroutine num_init(CCTK_ARGUMENTS) use arrays implicit none DECLARE_CCTK_PARAMETERS DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_FUNCTIONS integer i,j,k,l,m integer nx0,ny0,nz0,sample integer ngeos,res,nprev integer nprocs,myproc integer ierr CCTK_REAL dx,dy,dz CCTK_REAL zero,half,one CCTK_REAL xmn,ymn,zmn,xmx,ymx,zmx CCTK_REAL dnx,dny,dnz CCTK_REAL, DIMENSION (:), ALLOCATABLE :: XS,YS,ZS,PXS,PYS,PZS,DTAU ! Description of variables: ! ! i,j,k,l,m Counters ! ! ngeos Total number of geodesics ! ! xmn Minimum value of x on the grid ! xmx Maximum value of x on the grid ! ! ymn Minimum value of y on the grid ! ymx Maximum value of y on the grid ! ! zmn Minimum value of z on the grid ! zmx Maximum value of z on the grid ! ! sample Down sample size for: geopos = grid ! ! dnx Down-sampled number of points in x direction ! dny Down-sampled number of points in y direction ! dnz Down-sampled number of points in z direction ! ************************* ! *** PRINT MESSAGE *** ! ************************* call CCTK_INFO('Initializing timelike geodesics') ! *************************** ! *** DEFINE NUMBERS *** ! *************************** zero = 0.0d0 half = 0.5d0 one = 1.0d0 ! Grid parameters. dx = cctk_delta_space(1) dy = cctk_delta_space(2) dz = cctk_delta_space(3) ! ******************************************* ! *** FIND {xmn,xmx,ymn,ymx,zmn,zmx} *** ! ******************************************* ! Get range of coordinates. call CCTK_CoordRange(ierr,cctkGH,xmn,xmx,-1,"x","cart3d") call CCTK_CoordRange(ierr,cctkGH,ymn,ymx,-1,"y","cart3d") call CCTK_CoordRange(ierr,cctkGH,zmn,zmx,-1,"z","cart3d") ! Move safely away from boundaries. xmn = xmn + 2.0d0*dx ymn = ymn + 2.0d0*dy zmn = zmn + 2.0d0*dz xmx = xmx - 2.0d0*dx ymx = ymx - 2.0d0*dy zmx = zmx - 2.0d0*dz ! ************************************************* ! *** DISTRIBUTE GEODESICS AMONG PROCESSORS *** ! ************************************************* ! Find out total number of processors and identity of ! local processor. nprocs = CCTK_nProcs(cctkGH) myproc = CCTK_MyProc(cctkGH) ! Get the total number of geodesics. if (.not.CCTK_EQUALS(geopos,'grid')) then ! If the distribution of geodesics is not a down-sampled ! 3D array, the total number of geodesics will be given ! directly by the parameter "geodesics" ngeos = geodesics else ! If the distribution of geodesics is a down-sampled 3D array, ! the value of "geodesics" is ignored, and the total number of ! geodesics has to be computed from the total number of grid ! points. ! Find global number of grid points. nx0 = cctk_gsh(1) ny0 = cctk_gsh(2) nz0 = cctk_gsh(3) ! Find down-sampling ratio. sample = downsample ! Find total number of geodesics. dnx = (nx0-1)/sample + 1 dny = (ny0-1)/sample + 1 dnz = (nz0-1)/sample + 1 ngeos = dnx*dny*dnz end if ! Calculate the local number of geodesics per processor "lgeos". ! First, just do an integer division. lgeos = ngeos/nprocs ! Now, figure out how many geodesics are left over. res = ngeos - lgeos*nprocs ! And assign those geodesics one to each processor until we ! run out of geodesics. if (myproc+1.le.res) then lgeos = lgeos + 1 end if ! In order to assign the initial positions of the geodesics, ! we need to know how many were given to the previous processors. if (myproc+1.le.res) then nprev = (ngeos/nprocs + 1)*myproc else nprev = (ngeos/nprocs + 1)*res . + (ngeos/nprocs)*(myproc-res) end if ! ************************************** ! *** ALLOCATE MEMORY FOR ARRAYS *** ! ************************************** allocate(geoflag(lgeos)) allocate(geotau(lgeos),geodtau(lgeos),geodttau(lgeos)) allocate(geoxn(lgeos),geoyn(lgeos),geozn(lgeos)) allocate(geodtxn(lgeos),geodtyn(lgeos),geodtzn(lgeos)) allocate(geopxn(lgeos),geopyn(lgeos),geopzn(lgeos)) allocate(geodtpxn(lgeos),geodtpyn(lgeos),geodtpzn(lgeos)) allocate(XS(lgeos),YS(lgeos),ZS(lgeos),DTAU(lgeos)) allocate(PXS(lgeos),PYS(lgeos),PZS(lgeos)) ! ********************************** ! *** INITIALIZE PROPER TIME *** ! ********************************** geotau = zero ! ****************************************** ! *** INITIAL POSITIONS OF GEODESICS *** ! ****************************************** ! The initial position of the geodesics depend ! on the values of the parameter "geopos". ! x axis. if (CCTK_EQUALS(geopos,'x00')) then do i=1,lgeos j = i + nprev geoxn(i) = xmn + dble(j)*(xmx - xmn)/dble(ngeos) end do geoyn = zero geozn = zero call CCTK_INFO('Initial position of geodesics along x axis') ! y axis. else if (CCTK_EQUALS(geopos,'0y0')) then do i=1,lgeos j = i + nprev geoyn(i) = ymn + dble(j)*(ymx - ymn)/dble(ngeos) end do geoxn = zero geozn = zero call CCTK_INFO('Initial position of geodesics along y axis') ! z axis. else if (CCTK_EQUALS(geopos,'00z')) then do i=1,lgeos j = i + nprev geozn(i) = zmn + dble(j)*(zmx - zmn)/dble(ngeos) end do geoxn = zero geoyn = zero call CCTK_INFO('Initial position of geodesics along z axis') ! xy diagonal. else if (CCTK_EQUALS(geopos,'xy0')) then geozn = zero do i=1,lgeos j = i + nprev geoxn(i) = xmn + dble(j)*(xmx - xmn)/dble(ngeos) geoyn(i) = ymn + dble(j)*(ymx - ymn)/dble(ngeos) end do call CCTK_INFO('Initial position of geodesics along xy diagonal') ! xz diagonal. else if (CCTK_EQUALS(geopos,'x0z')) then geoyn = zero do i=1,lgeos j = i + nprev geoxn(i) = xmn + dble(j)*(xmx - xmn)/dble(ngeos) geozn(i) = zmn + dble(j)*(zmx - zmn)/dble(ngeos) end do call CCTK_INFO('Initial position of geodesics along xz diagonal') ! yz diagonal. else if (CCTK_EQUALS(geopos,'0yz')) then geoxn = zero do i=1,lgeos j = i + nprev geoyn(i) = ymn + dble(j)*(ymx - ymn)/dble(ngeos) geozn(i) = zmn + dble(j)*(zmx - zmn)/dble(ngeos) end do call CCTK_INFO('Initial position of geodesics along yz diagonal') ! xyz diagonal. else if (CCTK_EQUALS(geopos,'xyz')) then do i=1,lgeos j = i + nprev geoxn(i) = xmn + dble(j)*(xmx - xmn)/dble(ngeos) geoyn(i) = ymn + dble(j)*(ymx - ymn)/dble(ngeos) geozn(i) = zmn + dble(j)*(zmx - zmn)/dble(ngeos) end do call CCTK_INFO('Initial position of geodesics along xyz diagonal') ! Geodesics in 3D with downsample parameter. else if (CCTK_EQUALS(geopos,'grid')) then call CCTK_INFO('Initial position of geodesics as down-sampled 3D array') do l=1,lgeos m = l + nprev - 1 k = m/(dnx*dny) j = (m - k*dnx*dny)/dnx i = m - dnx*(k*dny + j) geoxn(l) = xmn + dble(i*sample)*(xmx - xmn)/dble(nx0-1) geoyn(l) = ymn + dble(j*sample)*(ymx - ymn)/dble(ny0-1) geozn(l) = zmn + dble(k*sample)*(zmx - zmn)/dble(nz0-1) end do ! Read from file. else if (CCTK_EQUALS(geopos,'file')) then call CCTK_INFO('Reading positions of geodesics from file') call CCTK_WARN(0,'Option not yet implemented!') else call CCTK_WARN(0,'Invalid value for geopos') end if ! ***************************************** ! *** INITIAL MOMENTUM OF GEODESICS *** ! ***************************************** ! Notice that the variables I use here correspond to the ! spatial components of the covariant momentum. Since the ! geodesics are moving initially in the normal direction to ! the hypersurfaces, those components are zero. geopxn = zero geopyn = zero geopzn = zero ! ************************************ ! *** INITIAL TIME DERIVATIVES *** ! ************************************ ! On input, {XS,YS,ZS} and {PXS,PYS,PZS} are ! the position and momentum respectively. ! On output, they are the time derivatives of ! the position and momentum respectively, and ! DTAU is the time derivative of the proper time. XS = geoxn YS = geoyn ZS = geozn PXS = geopxn PYS = geopyn PZS = geopzn call num_dt(CCTK_ARGUMENTS,XS,YS,ZS,PXS,PYS,PZS,DTAU,lgeos) geodtxn = XS geodtyn = YS geodtzn = ZS geodtpxn = PXS geodtpyn = PYS geodtpzn = PZS geodttau = DTAU ! *************************************** ! *** FLAG ALL GEODESICS AS VALID *** ! *************************************** ! All geodesics are initially inside the grid, so all flags ! are set to 1. Notice that using a real here instead of a ! logical is a gross waste of memory, but those are the only ! global arrays supported at the moment by cactus. geoflag = one if (out_xgraph.eq.1 .and. CCTK_MyProc(cctkGH).eq.0) then ! ****************** ! *** OUTPUT *** ! ****************** open(1,file='geoxn',form='formatted',status='replace') open(2,file='geoyn',form='formatted',status='replace') open(3,file='geozn',form='formatted',status='replace') write(1,"(A7,ES14.6)") """Time =",cctk_time write(2,"(A7,ES14.6)") """Time =",cctk_time write(3,"(A7,ES14.6)") """Time =",cctk_time do i=1,lgeos write(1,"(I10,ES14.6)") i,geoxn(i) write(2,"(I10,ES14.6)") i,geoyn(i) write(3,"(I10,ES14.6)") i,geozn(i) end do write(1,*) write(2,*) write(3,*) close(1) close(2) close(3) end if ! ****************************************************** ! *** IF TESTS ARE ACTIVE SAVE INITIAL POSITIONS *** ! ****************************************************** if (num_test.eq.1) then call CCTK_INFO('Test of numerical geodesics activated') ! Allocate memory for test arrays. allocate(geoxni(lgeos),geoyni(lgeos),geozni(lgeos)) allocate(test1_x(lgeos),test1_y(lgeos),test1_z(lgeos)) allocate(test2(lgeos),test3(lgeos),test4(lgeos)) ! Save initial position of geodesics. geoxni = geoxn geoyni = geoyn geozni = geozn ! Initialize test arrays. test1_x = zero test1_y = zero test1_z = zero test2 = zero test3 = zero test4 = zero ! Write test arrays to file. open(1,file='test1_x',form='formatted',status='replace') open(2,file='test1_y',form='formatted',status='replace') open(3,file='test1_z',form='formatted',status='replace') open(4,file='test2',form='formatted',status='replace') open(5,file='test3',form='formatted',status='replace') open(6,file='test4',form='formatted',status='replace') write(1,"(A7,ES14.6)") """Time =",cctk_time write(2,"(A7,ES14.6)") """Time =",cctk_time write(3,"(A7,ES14.6)") """Time =",cctk_time write(4,"(A7,ES14.6)") """Time =",cctk_time write(5,"(A7,ES14.6)") """Time =",cctk_time write(6,"(A7,ES14.6)") """Time =",cctk_time do i=1,lgeos write(1,"(I10,ES14.6)") i,test1_x(i) write(2,"(I10,ES14.6)") i,test1_y(i) write(3,"(I10,ES14.6)") i,test1_z(i) write(4,"(I10,ES14.6)") i,test2(i) write(5,"(I10,ES14.6)") i,test3(i) write(6,"(I10,ES14.6)") i,test4(i) end do write(1,*) write(2,*) write(3,*) write(4,*) write(5,*) write(6,*) close(1) close(2) close(3) close(4) close(5) close(6) end if ! ************************************************** ! *** IF WE WANT INVARIANTS, ALLOCATE MEMORY *** ! ************************************************** if (findnum_IJ.eq.1) then allocate(geonumI_r(lgeos),geonumI_i(lgeos)) allocate(geonumJ_r(lgeos),geonumJ_i(lgeos)) end if ! *************** ! *** END *** ! *************** return end