implicit none integer i,j,k,n,nmax,kdum,kp1,kp2,nsteps parameter (nmax=9000) real*8 a(8),TWOPI,ri,ran2 parameter (TWOPI=2.0*3.1415927) real*8 tcomp,tcomp0 integer npipe,g6_npipes,g6_set_j_particle,g6calc_lasthalf real*8 x,xdot,f,fdot,phi,body,eps2 parameter (eps2=1.E-6) common /pardata/ x(3,nmax),xdot(3,nmax),body(nmax),f(3,nmax), * fdot(3,nmax),phi(nmax),n real*8 time,told,dt,result,etot0 common /edata/ etot0,time,dt real*8 xsg,ysg,zsg,dis c GRAPE6 constants integer clusid,g6tunit,g6xunit,ipmax parameter (clusid = 0) parameter (g6tunit = 51) parameter (g6xunit = 51) parameter (ipmax = 100) c variables used in connection with GRAPE real*8 xgrape(3,ipmax),vgrape(3,ipmax),fold(3,ipmax), * fdold(3,ipmax),phiold(ipmax),acc(3,ipmax),jerk(3,ipmax), * potg(ipmax),h2 parameter (h2=0.d0) ! Neighbour radius (not needed here) integer index(ipmax),nsend c choose number of stars n = 8192 ! Number of stars kdum = -1 ! Random number seed c setup of a Plummer model (from NBODY4) do i = 1,n body(i)=1.d0/(1.d0*n) 30 A(1) = RAN2(KDUM) IF (A(1).LT.1.0D-10) GO TO 30 RI = (A(1)**(-0.6666667) - 1.0)**(-0.5) * A(2) = RAN2(KDUM) A(3) = RAN2(KDUM) X(3,I) = (1.0 - 2.0*A(2))*RI X(1,I) = SQRT(RI**2 - X(3,I)**2)*COS(TWOPI*A(3)) X(2,I) = SQRT(RI**2 - X(3,I)**2)*SIN(TWOPI*A(3)) 32 A(4) = RAN2(KDUM) A(5) = RAN2(KDUM) A(6) = A(4)**2*(1.0 - A(4)**2)**3.5 IF (0.1*A(5).GT.A(6)) GO TO 32 * A(8) = A(4)*SQRT(2.0)/(1.0 + RI**2)**0.25 A(6) = RAN2(KDUM) A(7) = RAN2(KDUM) XDOT(3,I) = (1.0 - 2.0*A(6))*A(8) XDOT(1,I) = SQRT(A(8)**2 - XDOT(3,I)**2)*COS(TWOPI*A(7)) XDOT(2,I) = SQRT(A(8)**2 - XDOT(3,I)**2)*SIN(TWOPI*A(7)) end do c initializing GRAPE c open GRAPE call g6_open(clusid) c get number of pipes npipe = g6_npipes() write (*,*) 'Got Grape6 - npipes = ',npipe c Set time and space units on GRAPE6 call g6_set_tunit(g6tunit) call g6_set_xunit(g6xunit) c estimate values for F,FDOT and PHI on host call estimate(1,n) c setting integration constants time = 0.d0 dt = 1.d0/2.d0**10 etot0 = 0.d0 c get cpu time elapsed so far call cputim(tcomp0) c start of integrator c set new time 10 time = time + dt c set time on grape call g6_set_ti(clusid, time) c submitting background particles to GRAPE, note the i-1 do i=1,n result = g6_set_j_particle(clusid, i-1, i, time, * dt, body(i), 0.d0, fdot(1,i), f(1,i), * xdot(1,i), x(1,i)) end do c calculate f and fdot of all particles on GRAPE kp1 = 1 20 kp2 = min(kp1+npipe-1,n) do i=kp1,kp2 j = i-kp1+1 do k=1,3 xgrape(k,j) = x(k,i) vgrape(k,j) = xdot(k,i) index(j) = i fold(k,j) = f(k,i) fdold(k,j) = fdot(k,i)/5.d0 phiold(j) = phi(i) end do end do nsend = kp2-kp1+1 call g6calc_firsthalf(clusid, n, nsend, index, xgrape, vgrape, * fold, fdold, phiold, eps2, h2) result=g6calc_lasthalf(clusid,n, nsend, index, xgrape, vgrape, * eps2, h2, acc, jerk, potg) do i=kp1,kp2 j = i-kp1+1 do k=1,3 f(k,i) = acc(k,j) fdot(k,i) = jerk(k,j) phi(i) = potg(j) end do end do kp1=kp2+1 if (kp2.lt.n) goto 20 c advance velocities do i=1,n do k=1,3 xdot(k,i)=xdot(k,i)+dt*f(k,i) end do end do c check for output if (dmod(time*128.d0,1.d0).eq.0.d0) then call auswert ! Check total energy call cputim(tcomp) if (time.gt.0.0) * write (*,'(a,f10.5,a,f12.5,a)') "Time/NBODY: ", * (tcomp-tcomp0)/(time-told)," Min Speed: ", * 57.d0*nsteps*(1.d0*n)**2/(60.d0*(tcomp-tcomp0))/1.E9, * " GFLOPS" tcomp0 = tcomp told = time nsteps = 0 end if c advance positions do i=1,n do k=1,3 x(k,i)=x(k,i)+dt*xdot(k,i) end do end do nsteps = nsteps+1 goto 10 stop end subroutine auswert implicit none integer i,j,k,n,nmax parameter (nmax=9000) real*8 x,xdot,f,fdot,phi,body,time,xdot2(3,nmax) common /pardata/ x(3,nmax),xdot(3,nmax),body(nmax),f(3,nmax), * fdot(3,nmax),phi(nmax),n real*8 epot,ekin,etot0,dt common /edata/ etot0,time,dt c move velocities back by half a step do i=1,n do k=1,3 xdot2(k,i)=xdot(k,i)-0.5*dt*f(k,i) end do end do epot = 0.d0 ekin = 0.d0 do i=1,n epot = epot + 0.5d0*body(i)*phi(i) ekin = ekin + 0.5*body(i)*(xdot2(1,i)**2+xdot2(2,i)**2+ * xdot2(3,i)**2) end do if (etot0.eq.0.d0) etot0 = epot+ekin write(*,'(/4(a,f12.7))') "T = ",time," etot: ",epot+ekin, * " de: ",(epot+ekin-etot0)/etot0," vir: ",-epot/(2.d0*ekin) etot0 = epot+ekin return end FUNCTION RAN2(IDUM) * * * Random number generator (Press p. 195). * --------------------------------------- * PARAMETER (M=714025,IA=1366,IC=150889,RM=1./M) COMMON/RAND2/ IY,IFF,IR(97) * DATA IFF /0/ * * IF (IDUM.LT.0.OR.IFF.EQ.0) THEN IFF = 1 IDUM = MOD(IC-IDUM,M) DO 11 J = 1,97 IDUM = MOD(IA*IDUM+IC,M) IR(J) = IDUM 11 CONTINUE IDUM = MOD(IA*IDUM+IC,M) IY = IDUM END IF J = 1 + (97*IY)/M IF (J.GT.97.OR.J.LT.1) WRITE (6,12) J, IDUM 12 FORMAT (/,' TROUBLES IN RAN2 J IDUM ',2I12) IY = IR(J) RAN2 = IY*RM IDUM = MOD(IA*IDUM+IC,M) IR(J) = IDUM RETURN END c determine cpu-time SUBROUTINE CPUTIM(TCOMP) real*8 TCOMP REAL*4 TARRAY(2) integer ierr,group,isize,rank TCOMP = ETIME(TARRAY)/60.d0 RETURN END c estimate f, fdot and phi for first use on GRAPE subroutine estimate (i1,i2) implicit none integer i,j,k,i1,i2,n,nmax,nmx parameter (nmax=9000) real*8 x,xdot,f,fdot,phi,body,eps2,vxsg,vysg,vzsg parameter (eps2=1.E-6) common /pardata/ x(3,nmax),xdot(3,nmax),body(nmax),f(3,nmax), * fdot(3,nmax),phi(nmax),n real*8 dis,xsg,ysg,zsg,rsg if (i.lt.10000) then nmx = n else nmx = 1000 end if do i=i1,i2 do k=1,3 f(k,i) = 0.d0 fdot(k,i) = 0.d0 end do phi(i) = 0.d0 do j=1,nmx if (j.ne.i) then dis = sqrt((x(1,i)-x(1,j))**2+(x(2,i)-x(2,j))**2+ * (x(3,i)-x(3,j))**2+eps2) xsg = x(1,i)-x(1,j) ysg = x(2,i)-x(2,j) zsg = x(3,i)-x(3,j) vxsg = xdot(1,i)-xdot(1,j) vysg = xdot(2,i)-xdot(2,j) vzsg = xdot(3,i)-xdot(3,j) rsg = vxsg*xsg+vysg*ysg+vzsg*zsg f(1,i) = f(1,i)-body(j)*xsg/dis**3 f(2,i) = f(2,i)-body(j)*ysg/dis**3 f(3,i) = f(3,i)-body(j)*zsg/dis**3 c Note: Usage of FDOT not necessary in leapfrog integrator. Only done c to demonstrate its use for higher order integrators fdot(1,i)= fdot(1,i) - body(j)*vxsg/dis**3+3.d0* * body(j)*rsg*xsg/dis**5 fdot(2,i)= fdot(2,i) - body(j)*vysg/dis**3+3.d0* * body(j)*rsg*ysg/dis**5 fdot(3,i)= fdot(3,i) - body(j)*vzsg/dis**3+3.d0* * body(j)*rsg*zsg/dis**5 phi(i)=phi(i)-body(j)/dis end if end do end do do i=1,n do k=1,3 f(k,i) = f(1,i)*(1.d0*n)/(1.d0*nmx) fdot(k,i) = fdot(1,i)*(1.d0*n)/(1.d0*nmx) end do phi(i) = phi(i)*(1.d0*n)/(1.d0*nmx) end do return end