/*@@ @file num_integrate.F @date October 1999 @author Miguel Alcubierre @desc Integrate geodesics in numerical spacetime. @enddesc @version $Id$ @@*/ #include "cctk.h" #include "cctk_Parameters.h" #include "cctk_Arguments.h" #include "cctk_Functions.h" /*@@ @routine num_integrate @date August 1997 @author Miguel Alcubierre @desc This routine integrates geodesics in the numerical spacetime using a second order Runge-Kutta method (well ... sort of). The geodesics are integrated using coordinate time as a parameter. It should be called after the main code has updated the geometry (including lapse and shift). @enddesc @calls timegeodesic_num_dt @@*/ subroutine num_integrate(CCTK_ARGUMENTS) use arrays implicit none DECLARE_CCTK_PARAMETERS DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_FUNCTIONS integer i,ierr integer nx,ny,nz integer npoints integer ierror,alp_index integer gxx_index,gyy_index,gzz_index integer gxy_index,gxz_index,gyz_index integer vindex,interp_handle,coord_system_handle,param_table_handle character(30) options_string 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 xmn,ymn,zmn CCTK_REAL xmx,ymx,zmx CCTK_REAL zero,half,one CCTK_REAL dt,dth CCTK_REAL rm2 CCTK_REAL rr,sxy,cost,cosp,sint,sinp,gqq CCTK_REAL, dimension (lgeos) :: XS,YS,ZS,PXS,PYS,PZS,DTAU CCTK_REAL, dimension (lgeos) :: A,G1,G2,G3,G4,G5,G6 ! Description of variables: ! ! i,j,k Counters ! ! xmn Minimum value of x on the grid. ! ymn Minimum value of y on the grid. ! zmn Minimum value of z on the grid. ! ! xmx Maximum value of x on the grid. ! ymx Maximum value of y on the grid. ! zmx Maximum value of z on the grid. ! ! XS Array that holds x position of the geodesics ! on input to the subroutine timegeodesic_dt, ! and its time derivative on output. ! ! YS Array that holds y position of the geodesics ! on input to the subroutine timegeodesic_dt, ! and its time derivative on output. ! ! ZS Array that holds z position of the geodesics ! on input to the subroutine timegeodesic_dt, ! and its time derivative on output. ! ! PXS Array that holds x mementum of the geodesics ! on input to the subroutine timegeodesic_dt, ! and its time derivative on output. ! ! PYS Array that holds y momentum of the geodesics ! on input to the subroutine timegeodesic_dt, ! and its time derivative on output. ! ! PZS Array that holds z momentum of the geodesics ! on input to the subroutine timegeodesic_dt, ! and its time derivative on output. ! ! DTAU Array that contains time derivatives of the ! proper time on output from timegeodesic_dt. ! ! dth dt/2 ! ! rm2 Square of the Minkowski radius of geodesics ! (only used in test4). ! ************************************* ! *** IF iteration = 0 RETURN *** ! ************************************* ! Because of the calling structure of CACTUS, this routine ! is called even after time step 0. This should not happen ! because the geodesics will not be sincronized properly ! any more. To avoid it, I make the routine return here ! if iteration = 0. if (cctk_iteration.eq.0) return ! ************************** ! *** DEFINE NUMBERS *** ! ************************** zero = 0.0d0 half = 0.5d0 one = 1.0d0 ! Grid parameters. nx = cctk_lsh(1) ny = cctk_lsh(2) nz = cctk_lsh(3) dt = cctk_delta_time dth = half*dt call UpdateGeodesics(cctkGH, lgeos, geotau, XS,YS,ZS, PXS, PYS, PZS, geoflag) ! ******************************************* ! *** FIND {xmn,xmx,ymn,ymx,zmn,zmx} *** ! ******************************************* 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") ! ******************************** ! *** DO A FULL EULER STEP *** ! ******************************** ! Do a full Euler step for the position and store the ! predicted position in {XS,YS,ZS}. XS = geoxn + dt*geodtxn YS = geoyn + dt*geodtyn ZS = geozn + dt*geodtzn ! Do a full Euler step for the momentum and store the ! predicted momentum in {PXS,PYS,PZS}. PXS = geopxn + dt*geodtpxn PYS = geopyn + dt*geodtpyn PZS = geopzn + dt*geodtpzn ! Check for geodesics that have left the grid. do i=1,lgeos if ((geoflag(i).lt.half).or. & (XS(i).le.xmn).or.(XS(i).ge.xmx).or. & (YS(i).le.ymn).or.(YS(i).ge.ymx).or. & (ZS(i).le.zmn).or.(ZS(i).ge.zmx)) then geoflag(i) = zero XS(i) = geoxn(i) YS(i) = geoyn(i) ZS(i) = geozn(i) PXS(i) = zero PYS(i) = zero PZS(i) = zero end if end do ! ******************************************* ! *** FIND PREDICTED 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. call num_dt(CCTK_ARGUMENTS,XS,YS,ZS,PXS,PYS,PZS,DTAU,lgeos) ! ******************************************* ! *** FIND THE AVERAGE OF THE OLD AND *** ! *** PREDICTED TIME DERIVATIVES *** ! ******************************************* geodtxn = half*(geodtxn + XS) geodtyn = half*(geodtyn + YS) geodtzn = half*(geodtzn + ZS) geodtpxn = half*(geodtpxn + PXS) geodtpyn = half*(geodtpyn + PYS) geodtpzn = half*(geodtpzn + PZS) ! ********************************************** ! *** ADVANCE A FULL TIME STEP USING THE *** ! *** AVERAGED TIME DERIVATIVES *** ! ********************************************** ! Advance the position. XS = geoxn + dt*geodtxn YS = geoyn + dt*geodtyn ZS = geozn + dt*geodtzn ! Advance the momentum. PXS = geopxn + dt*geodtpxn PYS = geopyn + dt*geodtpyn PZS = geopzn + dt*geodtpzn ! Check for geodesics that have left the grid. do i=1,lgeos if ((geoflag(i).lt.half).or. & (XS(i).le.xmn).or.(XS(i).ge.xmx).or. & (YS(i).le.ymn).or.(YS(i).ge.ymx).or. & (ZS(i).le.zmn).or.(ZS(i).ge.zmx)) then geoflag(i) = zero XS(i) = geoxn(i) YS(i) = geoyn(i) ZS(i) = geozn(i) PXS(i) = zero PYS(i) = zero PZS(i) = zero end if end do ! *************************************************** ! *** UPDATE THE POSITION AND MOMENTUM ARRAYS *** ! *************************************************** geoxn = XS geoyn = YS geozn = ZS geopxn = PXS geopyn = PYS geopzn = PZS ! ************************************* ! *** FIND NEW 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. 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 ! ********************************** ! *** UPDATE THE PROPER TIME *** ! ********************************** ! Find the increment in proper time. geodtau = dth*(geodttau + DTAU) ! Update the proper time. geotau = geotau + geodtau ! Update the time derivative of the proper time. geodttau = DTAU call OutputGeodesics(cctkGH, lgeos, geotau, geoxn, geoyn, geozn, geopxn, geopyn, geopzn, geoflag) if (out_xgraph.eq.1 .and. CCTK_MyProc(cctkGH).eq.0) then ! ****************** ! *** OUTPUT *** ! ****************** open(1,file='geoxn',form='formatted', . status='old',position='append') open(2,file='geoyn',form='formatted', . status='old',position='append') open(3,file='geozn',form='formatted', . status='old',position='append') 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 ! ***************************** ! *** CONSISTENCY TESTS *** ! ***************************** ! Here I do some consistency tests for the solution. if (num_test.eq.1) then XS = geoxn YS = geoyn ZS = geozn PXS = geopxn PYS = geopyn PZS = geopzn ! Interpolate to find lapse and metric. npoints = lgeos ! 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_InterpHandle (interp_handle, "uniform cartesian") 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(XS) interp_coords(2) = CCTK_PointerTo(YS) interp_coords(3) = CCTK_PointerTo(ZS) call CCTK_VarIndex (vindex, "admbase::alp") in_array_indices(1) = vindex call CCTK_VarIndex (vindex, "admbase::gxx") in_array_indices(2) = vindex call CCTK_VarIndex (vindex, "admbase::gyy") in_array_indices(3) = vindex call CCTK_VarIndex (vindex, "admbase::gzz") in_array_indices(4) = vindex call CCTK_VarIndex (vindex, "admbase::gxy") in_array_indices(5) = vindex call CCTK_VarIndex (vindex, "admbase::gxz") in_array_indices(6) = vindex call CCTK_VarIndex (vindex, "admbase::gyz") in_array_indices(7) = vindex out_arrays(1) = CCTK_PointerTo(A) out_arrays(2) = CCTK_PointerTo(G1) out_arrays(3) = CCTK_PointerTo(G2) out_arrays(4) = CCTK_PointerTo(G3) out_arrays(5) = CCTK_PointerTo(G4) out_arrays(6) = CCTK_PointerTo(G5) out_arrays(7) = CCTK_PointerTo(G6) out_array_type_codes = CCTK_VARIABLE_REAL 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 call Util_TableDestroy (ierror, param_table_handle) ! TEST I: The first test is for consistency between covariant ! and contravariant momentum. ! ! mu 0 / k .k \ ! p = gamma p => p - p g | b + x | = 0 ! i i mu i ik \ / ! ! ! Notice that this should be true for any curve, regardless ! of whether it is a geodesic or not. Notice also that since ! the time derivative of the position was calculated in terms ! of the other quantities, this condition should be satisfied ! up to round of error regardless of resolution (unless there ! is a bug). ! ! I calculate the difference above and store it in the arrays ! {test1_x,test1_y,test1_z}. test1_x = PXS - (G1*geodtxn + G4*geodtyn + G5*geodtzn)/DTAU test1_y = PYS - (G4*geodtxn + G2*geodtyn + G6*geodtzn)/DTAU test1_z = PZS - (G5*geodtxn + G6*geodtyn + G3*geodtzn)/DTAU c if (shift_active.ne.SHIFT_INACTIVE) then c call CCTK_INFO('Shift terms missing in TEST I') c end if ! TEST II: The second test is to check if the curves are ! in fact timelike and are properly normalized. ! ! 0 k ! p p + p p = - 1 ! 0 k ! ! 0 / 2 0 .k \ 0 2 i / j .j \ ! => p | - a p + p x | + (p ) g b | b + x | + 1 = 0 ! \ k / ij \ / ! ! I calculate the left hand side above and store it in the ! array {test2}. Again, since the time derivative of the ! position was calculated in terms of the other quantities, ! this condition should be satisfied up to round of errors ! regardless of the resolution. test2 = (PXS*geodtxn + PYS*geodtyn + PZS*geodtzn & - A**2/DTAU)/DTAU + one c if (shift_active.ne.SHIFT_INACTIVE) then c call CCTK_INFO('Shift terms missing in TEST II') c end if ! TEST III: Only valid in model "flatty". Because the ! initial lapse is spherically symmetric, the coordinate ! speed of the geodesics must point in the radial direction. ! This implies that: ! . . ! x . x - |x| |x| = 0 ! ! where the dot product and the norms are those for flat space. ! I expect these to converge to high order, depending on the ! specific numerical scheme. ! ! I calculate the difference and store it in {test3}. test3 = (XS*geodtxn + YS*geodtyn + ZS*geodtzn) & - dsqrt((XS**2 + YS**2 + ZS**2) & *(geodtxn**2 + geodtyn**2 + geodtzn**2)) ! TEST IV: Only valid in model "flatty". Because the ! initial lapse is spherically symmetric, the angular ! metric component g_{theta,theta} will always be equal ! to the square of the Minkowski radius. Now since the ! geodesics remain in the same place in Minkowski coordinates, ! we must have: ! 2 ! g - rm = 0 ! qq ! ! where g_qq is short for g_{theta,theta} and r is the Minkowski ! radius of the geodesics. This is a non-trivial test that ! really checks that the geodesics do what they should in a case ! where we know their exact location (in Minkowski coordinates ! they have not moved at all). I expect this relation to converge ! to second order. ! ! I calculate the difference and store it in {test4}. do i=1,lgeos sxy = dsqrt(XS(i)**2 + YS(i)**2) rr = dsqrt(sxy**2 + ZS(i)**2) ! Find sines and cosines. if (sxy.eq.zero) then ! Axis. cost = one sint = zero sinp = zero cosp = one else ! Away from axis. cost = ZS(i)/rr sint = sxy/rr sinp = YS(i)/sxy cosp = XS(i)/sxy end if ! Find g_{theta,theta} at the position of ! the geodesics. gqq = rr**2*(sint**2*G3(i) & + cost**2*(cosp**2*G1(i) + sinp**2*G2(i)) & + 2.0D0*cost*(cost*cosp*sinp*G4(i) & - sint*(cosp*G5(i) + sinp*G6(i)))) ! Find the Minkowski radius of the geodesic. ! This is obtained from the initial positions. rm2 = geoxni(i)**2 + geoyni(i)**2 + geozni(i)**2 ! Find difference between g_{theta,theta} and ! r**2. The difference should be zero since ! the Minkowski radius of the geodesics has not ! changed. test4(i) = gqq - rm2 end do ! Write test arrays to file. open(1,file='test1_x',form='formatted', . status='old',position='append') open(2,file='test1_y',form='formatted', . status='old',position='append') open(3,file='test1_z',form='formatted', . status='old',position='append') open(4,file='test2',form='formatted', . status='old',position='append') open(5,file='test3',form='formatted', . status='old',position='append') open(6,file='test4',form='formatted', . status='old',position='append') 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 ! *************** ! *** END *** ! *************** return end