!--------------------------------------------------------------------------------------------------------- ! Main module for NBODY2: direct summation code with variable time-step ! ! Time-steps can change by factor two ! Maximum time step is dtstep ! Maximum possible level is iMultiple ! Maximum level at given moment is MaxLev ! If iMultiple =0, code uses only single constant time step ! In this case the scheme is Kick-Drift-Kick ! If iMultiple > 0, variable discrete time-stepping. ! scheme of stepping is Drift-Kick-Drift ! level of time-stepping is decide by SelectStep ! Level of each particles can change only by factor 2 ! Conditions for chosing a time-step: ! max(go,g1)dtstep < Alph*Velocity ! max(Vo,V1)dtstep < aStep ! LevPart(i) - level of each particles relative to MaxLev ! dt_particle =dtstep*2**(LevPart-MaxLev) ! LevPart =0 -> smallest step ! LevPart =1 -> twice bigger step and so on !---------------------------------------------------------------------------------------------------- MODULE DataStructures IMPLICIT DOUBLE PRECISION (a-h,o-z) Real(8), Parameter :: aStep = 2.d-3, & ! accuracy of displacements Alph = 5.d-4, & ! accuracy of velocities eps = 1.d-3, & ! force softening eps2 = eps**2 Integer, Parameter :: iMultiple = 4 ! 0 for single time-step ! if /= 0, max alowed level of time stepping Real(8) :: time = 0.d0, & dtstep =2.d-4, & ! Max time-step PrintTime = 0., & dPrintTime = 0.01, & RunTime = 1.57*80. Integer :: MaxLev ! current maximum level REAL*8, ALLOCATABLE, DIMENSION(:) :: xc,yc,zc, & vx,vy,vz, & gx,gy,gz,amass, & timeP, & xa,ya,za, & ! Auxiliary set of {x,v,g} vxa,vya,vza, & gxa,gya,gza Integer, ALLOCATABLE, DIMENSION(:) :: LevPart,iLevPart Logical :: iFlag =.true. ! flag for printing Integer, SAVE :: Np,iStep,Naccel=0 Integer :: iLevel(0:20) ! number of particles on different levels Contains !-------------------------------------------------------------------- ! Find Energy of the system SUBROUTINE Energies(Ekin,Epot) !-------------------------------------------------------------------- REAL*8 acc,rr,Ekin,Epot Ekin =0. Epot =0. !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP REDUCTION(+:Ekin) Do i=1,Np ! kinetic energy Ekin = Ekin +(vx(i)**2+vy(i)**2+vz(i)**2)*amass(i) EndDo Ekin = Ekin/2. !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE(i,j) SCHEDULE(DYNAMIC) & !$OMP REDUCTION(+:Epot) Do i=1,Np-1 ! sum contributions pair-wise Do j=i+1,Np Epot = Epot - amass(i)*amass(j)/ & sqrt( (xc(i)-xc(j))**2+ (yc(i)-yc(j))**2+ (zc(i)-zc(j))**2+ eps2 ) EndDo EndDo Return End SUBROUTINE Energies !-------------------------------------------------------------------- ! write snapshot of the system to file SUBROUTINE WriteSnapshot(center) !-------------------------------------------------------------------- REAL*8 center(3) Character*80 Name ! construct file name using time ! split: time = i1.i2 i1 = INT(time) i2 = MIN(INT((time-i1)*1000.+0.5),999) write(Name,'(a,i3.3,".",i3.3,a)')'DATA/Snap',i1,i2,'.dat' open(10,file=Name) write(10,'(i6,2g11.4,i8,3g12.4)') & Np,dtstep,time,iStep,eps,RunTime, & PrintTime Do i=1,Np write(10,20) xc(i)-center(1), yc(i)-center(2), zc(i)-center(3), & vx(i),vy(i),vz(i),amass(i),LevPart(i) EndDo close(10) 20 format(7g13.5,2x,i5) Return End SUBROUTINE WriteSnapshot !-------------------------------------------------------------------- ! save current state of the system ! use unformatted file 'FileSave' SUBROUTINE SaveMoment(FileSave) !-------------------------------------------------------------------- Character(*) :: FileSave open(30,file=TRIM(FileSave),form='UNFORMATTED') write(30) Np,dtstep,time,iStep,eps,RunTime,PrintTime Do i=1,Np write(30) xc(i), yc(i), zc(i), & vx(i),vy(i),vz(i), & gx(i),gy(i),gz(i), & amass(i),LevPart(i) EndDo close (30) Return End SUBROUTINE SaveMoment !-------------------------------------------------------------------- ! Read current state of the system ! Allocate arrays ! dump data into allocated arrays SUBROUTINE ReadMoment(FileSave) !-------------------------------------------------------------------- Character(*) :: FileSave open(30,file=TRIM(FileSave),form='UNFORMATTED') read(30) Np,dtstep,time,iStep,eps0,RunTime,PrintTime ALLOCATE (xc(Np),yc(Np), zc(Np)) ! allocate arrays ALLOCATE (xa(Np),ya(Np), za(Np)) ALLOCATE (vx(Np),vy(Np), vz(Np)) ALLOCATE (vxa(Np),vya(Np), vza(Np)) ALLOCATE (gx(Np),gy(Np), gz(Np)) ALLOCATE (gxa(Np),gya(Np), gza(Np)) ALLOCATE (timeP(Np),aMass(NP)) ALLOCATE (LevPart(Np),iLevPart(Np)) !$OMP PARALLEL DO DEFAULT(SHARED) Do i=1,Np ! use first touch rule to distribute the data xc(i) =0.; yc(i)=0.; zc(i) =0. xa(i) =0.; ya(i)=0.; za(i) =0. vx(i) =0.; vy(i)=0.; vz(i) =0. vxa(i) =0.; vya(i)=0.; vza(i) =0. gx(i) =0.; gy(i) =0.; gz(i) =0. gxa(i) =0.; gya(i)=0.; gza(i) =0. LevPart(i) =MIN(MaxLev,2) iLevPart(i) =0 timeP(i)= 0. EndDo Do i=1,Np read(30) xc(i), yc(i), zc(i), & vx(i),vy(i),vz(i), & gx(i),gy(i),gz(i), & amass(i),LevPart(i) EndDo close (30) Return End SUBROUTINE ReadMoment !-------------------------------------------------------------------- ! SUBROUTINE Analyze(indx,center) !-------------------------------------------------------------------- PARAMETER (Nbin =18) ! number of bins for mass profile REAL*8, DIMENSION(3) :: center, vaver, wMass, & center2, vaver2, wMass2 REAL*8 Rmax,Raverage, tMass,tMass2,pi Character*80 Name DIMENSION Radii(Nbin),aMrad(Nbin) DIMENSION dDen(-100:Nbin),rDen(-100:Nbin),nDen(-100:Nbin), & sDen(-100:Nbin),uDen(-100:Nbin) DATA Radii/0.05,0.10,0.16,0.2,0.25,0.3,0.35,0.4,0.5, & 0.6,0.80,1.0,2.0,4.0, & 5.0,10.0,15.0,20.0/ pi =4.*atan(1.d0) i1 = INT(time) i2 = MIN(INT((time-i1)*1000.+0.5),999) write(Name,'(a,i3.3,".",i3.3,a)')'DATA/Prof',i1,i2,'.dat' open(9,file=Name) Rtoolarge =Radii(Nbin)*10. drlog = 0.12 ! step in log for dens profile center =0.; center2 =0. vaver =0.; vaver2 =0.; wMass=0.; wMass2=0. If(indx.eq.0)Then write(8,120)(Radii(i),i=1,Nbin) indx = 1 120 format(T48,20f8.2) EndIf ! find center of mass and average velocity Raverage =0. Rmax =0. tMass2 =0. tMass =0. Do i=1,Np aM = aMass(i) tMass = tMass + aM rr = sqrt(xc(i)**2+yc(i)**2+zc(i)**2) center(1) =center(1) +xc(i)*aM center(2) =center(2) +yc(i)*aM center(3) =center(3) +zc(i)*aM vaver(1) = vaver(1) +vx(i)*aM vaver(2) = vaver(2) +vy(i)*aM vaver(3) = vaver(3) +vz(i)*aM If(rr.lt.Rtoolarge)Then tMass2 = tMass2 + aM center2(1) =center2(1) +xc(i)*aM center2(2) =center2(2) +yc(i)*aM center2(3) =center2(3) +zc(i)*aM vaver2(1) = vaver2(1) +vx(i)*aM vaver2(2) = vaver2(2) +vy(i)*aM vaver2(3) = vaver2(3) +vz(i)*aM endif EndDo ! normalize results If(tMass>1.d-10)Then center = center/tMass vaver = vaver/tMass EndIf If(tMass2>1.d-10)Then center2 = center2/tMass2 vaver2 = vaver2/tMass2 EndIf ! find center of the densiest part of the system Dens0 = tMass2/Rtoolarge**3 Do iter=1,150 Rtoolarge =Rtoolarge/1.1 Raverage =0. tMass2 =0. inside = 0 center = 0. vaver = 0. Do i=1,Np aM = aMass(i) dx = xc(i) -center2(1) dy = yc(i) -center2(2) dz = zc(i) -center2(3) rr = sqrt(dx**2+dy**2+dz**2) If(rr.lt.Rtoolarge)Then tMass2= tMass2 + aM inside= inside +1 center(1) =center(1) +xc(i)*aM center(2) =center(2) +yc(i)*aM center(3) =center(3) +zc(i)*aM vaver(1) = vaver(1) +vx(i)*aM vaver(2) = vaver(2) +vy(i)*aM vaver(3) = vaver(3) +vz(i)*aM EndIf EndDo If(tMass2>1.d-10)Then center2 = center/tMass2 vaver2 = vaver/tMass2 EndIf Dens1 = tMass2/Rtoolarge**3 !write(*,'(a,3f8.4,a,2f11.3,i6)')' center=',center2 , ' Rmax=',Rtoolarge,Dens1,inside If(Dens1.lt.Dens0.or.float(inside).lt.float(Np)/4.)exit Dens0 = Dens1 EndDo center =center2 aMrad =0. Rtoolarge =Radii(Nbin) tMass2 =0. Raverage=0. Rmax =0. dDen =0. rDen =0. nDen =0 sDen =0. uDen =0. Do i=1,Np ! get statistics relative to the center aM = aMass(i) dx = xc(i) -center2(1) dy = yc(i) -center2(2) dz = zc(i) -center2(3) rr = sqrt(dx**2+dy**2+dz**2) vv = vx(i)**2+vy(i)**2+vz(i)**2 vr = (vx(i)*dx+vy(i)*dy+vz(i)*dz) /max(rr,1.e-10) ind = INT(log10(rr)/drlog+1.e5)-1e5 ind = max(min(ind,Nbin),-100) dDen(ind) = dDen(ind) + aM rDen(ind) = rDen(ind) + rr*aM nDen(ind) = nDen(ind) + 1 sDen(ind) = sDen(ind) + vr*aM uDen(ind) = uDen(ind) + vv*aM If(rr.lt.Rtoolarge)Then tMass2 =tMass2 + aM Raverage = Raverage + rr**2*aM do j=1,Nbin if(rr.lt.Radii(j))Then aMrad(j) =aMrad(j) +aM goto 10 EndIf EndDo 10 Rmax =max(Rmax,rr) EndIf EndDo Do j=2,Nbin ! make mass cumulative aMrad(j) =aMrad(j) + aMrad(j-1) EndDo Raverage = sqrt(Raverage/tMass2) ! normalize results write(8,130)time,iStep,tMass,tMass2,Rmax,Raverage, & (aMrad(i),i=1,Nbin) 100 format(a,f8.4,a,i8,a,2f8.3,a,2f8.3,20f8.4) 130 format(f8.4,i8,2f8.3,2f8.3,20f8.4) 110 format(i4,2g12.4,3x,2g12.4,1x,g12.3,i6,4x,4f8.3) ! write density profile write(9,100)' t=',time,' Step=',iStep, & ' Mass=',tMass,tMass2, & ' Rmax/Raverage=',Rmax,Raverage write(9,140) 140 format(' bin',' ',7x,'rad_in',8x,' dens',7x,' NumDen', & 7x,' ',6x,' Npart',5x,'Vradial',3x,'Vrms',3x,'M(