!------------------------------------------------- ! Parallel (OMP) PM code ! 2015, 1997 Anatoly Klypin (aklypin@nmsu.edu) ! Astronomy Department, NMSU ! ! ! Read PM files and produce: ! density field FI(NGRID.NGRID,NGRID) ! Coodinates of particles XPAR(Nparticles),YPAR,ZPAR in Mpch ! Velocities VX,VY,VZ (km/s) ! Density at particle position pDens(Nparticles) in units of average density ==1 ! ! ! Input: snapshot number (integer). Reads files PMcrd.XXXX.DAT, PMcrs0.XXXX.DAT ... ! Output: Out.XXXX.dat ! Parameters: read from PMcrd file: ! AEXPN - expansion parameter ! iStep - snapshot number ! Om0 - dark matter density ! hubble- Hubble constant ! Box - box size in comoving Mpch ! ! NGRID - mesh size for density field ! NROW - number of particles in 1d ! Nparticles - total number of particles = Nrow**3 ! Memory: requires 7 arrays of sigle precision each of Nparticles length ! requires one single precision array of size Ngrid**3 ! total (7*4*Nparticles+4*Ngrid**3)/1024**3 Gb ! ! Compile: ! ifort -O2 -openmp -i-dynamic -mcmodel=medium -convert big_endian -o PMPread.exe PMPread.f90 ! ! use routine Galaxies to build your analysis tool !------------------------------------------------- ! ! NGRID= number of cells in one dimension ! AEXPN= expansion parameter = 1/(1+z) Module Structures Real*4 :: AEXPN,AEXP0,AMPLT,ASTEP, box, & EKIN,EKIN1,EKIN2,AEU0, & TINTG,Nspecs,Nseed, & Om0,Oml0,hubble, & extras(100), & ENKIN,ENPOT Integer*4 :: iStep,Ngrid,Nrow,moment Integer*8 :: Nparticles Real*4, Allocatable, Dimension(:,:,:) :: FI Real*4, Allocatable, Dimension(:) :: & XPAR,YPAR,ZPAR, & VX,VY,VZ,pDens Real*4, allocatable, Dimension(:) :: Xb,Yb,Zb,VXb,Vyb,Vzb Character*45 :: HEADER end Module Structures !------------------------------------------------- Program PMPread use Structures character*80 :: Name Integer*8 :: count Real*4 :: Memory ! Read data and open files WRITE (*,'(A,$)') 'Enter snapshot number => ' READ (*,*) moment write(Name,'(a,i4.4,a)')'Out.',moment,'.dat' OPEN(17,FILE=TRIM(Name),STATUS='UNKNOWN') Call Timing(0,-1) ! initiate timing counts CALL ReadDataPM count = NGRID Mem_current = Memory(count**3) ALLOCATE(FI(NGRID,NGRID,NGRID)) CALL DENSIT ! Define density CALL PartDENS ! density at paticle positions CALL Rescale ! set coords in Mpch com, V in km/s CALL Galaxies ! make galaxies Call Timing(0,1) ! finish total time Call Timing(0,0) ! print time END Program PMPread !--------------------------------------- ! Read PMfiles: PMcrd/crs SUBROUTINE ReadDataPM() !--------------------------------------- use Structures Character*80 :: Name Logical :: exst Integer*8 :: iCount,ii,ioff,ip Integer*8 :: Ngal,Nrecord,Jpage,NinPage Real*4 :: Memory ! Call Timing(1,-1) ! start reading time write(Name,'(a,i4.4,a)')'PMcrd.',moment,'.DAT' INQUIRE(file=TRIM(Name),EXIST=exst) If(.not.exst)Stop ' File PMcrd.moment.DAT does not exist. Error' Open (4,file =TRIM(Name),form ='UNFORMATTED',status ='UNKNOWN') READ (4) HEADER, & AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW, & TINTG,EKIN,EKIN1,EKIN2,AU0,AEU0, & NROW,NGRID,Nspecs,Nseed,Om0,Oml0,hubble, & Nparticles,extras WRITE (*,'(a,/2x,a,f8.4,4(a,i7))') HEADER, & ' a =',AEXPN, ' step= ',ISTEP, & ' Nrow= ', NROW, ' Ngrid=',NGRID WRITE (17,'(a)') HEADER close(4) Nrecord = 1024**2 Naccess = Nrecord*6 !*4 xR = NGRID +1 boxsize = extras(100) Box = boxsize Npages = (Nparticles-1)/Nrecord+1 ! number of records Nlast = Nparticles - (Npages-1)*Nrecord ! number of particles in last record Jpage = Nrecord write(*,'(a,i10)') ' NROW =',NROW write(*,'(a,i10)') ' Ngal =',Nparticles write(*,'(a,i10)') ' Ngrid =',Ngrid write(*,'(a,f10.1)') ' Box =',Box write(17,'(a,i10)') ' NROW =',NROW write(17,'(a,i10)') ' Ngal =',Nparticles write(17,'(a,i10)') ' Ngrid =',Ngrid write(17,'(a,f10.1)') ' Box =',Box write(17,'(a,f10.4)') ' Hubble =',hubble write(17,'(a,f10.4)') ' Om0 =',Om0 write(17,'(a,i12)') ' Nseed =',Nseed Mem_current = Memory(6_8*Nrecord) Allocate (Xb(Nrecord),Yb(Nrecord),Zb(Nrecord)) Allocate (VXb(Nrecord),VYb(Nrecord),VZb(Nrecord)) Mem_current = Memory(6_8*Nparticles) Allocate (XPAR(Nparticles),YPAR(Nparticles),ZPAR(Nparticles)) Allocate (VX(Nparticles),VY(Nparticles),VZ(Nparticles)) iCount = 0 ifile = 0 jj = 0 write(Name,'(a,i1.1,a,i4.4,a)')'PMcrs',ifile,'.',moment,'.DAT' write(*,'(a)') TRIM(Name) INQUIRE(file=TRIM(Name),EXIST=exst) If(.not.exst)Stop ' File PMcrs0.DAT does not exist. Error' OPEN(UNIT=20,FILE=TRIM(Name),ACCESS='DIRECT', & FORM='unformatted',STATUS='UNKNOWN',RECL=NACCESS) Do ii =1,Npages If(ii==Npages)Then NinPage = Nparticles -(ii-1)*JPAGE ! # particles in the current page Else NinPage = JPAGE EndIf jj = jj +1 If(ii<10.or.ii==Npages)write(*,'(3(a,i9))') ' Reading page= ',ii,' record =',jj,' NinPage= ',NinPage 10 Read(20,REC=jj,iostat=ierr) Xb,Yb,Zb,VXb,VYb,VZb If(ierr /=0)Then close(20) ifile = ifile +1 If(ifile<10)Then write(Name,'(a,i1.1,a,i4.4,a)')'PMcrs',ifile,'.',moment,'.DAT' Else write(Name,'(a,i2.2,a,i4.4,a)')'PMcrs',ifile,'.',moment,'.DAT' EndIf jj = 1 INQUIRE(file=TRIM(Name),EXIST=exst) If(.not.exst)Stop ' Error reading PMcrs files: Did not get all' Open(20,file=TRIM(Name),ACCESS='DIRECT', & FORM='unformatted',STATUS='UNKNOWN',RECL=NACCESS) write(*,'(2i7,2a,3x,i9)') ii,ifile,' Open file = ',TRIM(Name),Ninpage go to 10 end If ioff = (ii-1)*JPAGE !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE (ip) Do ip =1,NinPage If(ip+ioff > Nparticles)STOP 'Attempt to read too many particles ' if(INT(Xb(ip))==Ngrid+1)Xb(ip)=Xb(ip)-1.e-3 if(INT(Yb(ip))==Ngrid+1)Yb(ip)=Yb(ip)-1.e-3 if(INT(Zb(ip))==Ngrid+1)Zb(ip)=Zb(ip)-1.e-3 if(INT(Xb(ip))==Ngrid+1)write(*,*)'Error in boundary: ',INT(Xb(ip)),Xb(ip) if(INT(Yb(ip))==Ngrid+1)write(*,*)'Error in boundary: ',INT(Yb(ip)),Yb(ip) if(INT(Zb(ip))==Ngrid+1)write(*,*)'Error in boundary: ',INT(Zb(ip)),Zb(ip) Xpar(ip+ioff) = Xb(ip) Ypar(ip+ioff) = Yb(ip) Zpar(ip+ioff) = Zb(ip) VX(ip+ioff) = VXb(ip) VY(ip+ioff) = VYb(ip) VZ(ip+ioff) = VZb(ip) end Do end DO close (20) xx = MAXVAL(Xpar) xm = MINVAL(Xpar) write(*,*)' x min/max= ',xm,xx xx = MAXVAL(Ypar) xm = MINVAL(Ypar) write(*,*)' y min/max= ',xm,xx xx = MAXVAL(Zpar) xm = MINVAL(Zpar) write(*,*)' z min/max= ',xm,xx Do ii =1,Nparticles If(Xpar(ii).lt.1.0.or.Ypar(ii).lt.1.0.or.Zpar(ii).lt.1.0)& write(*,'(a,i10,1p,3g14.5)')' Error coord: ',ii,Xpar(ii),Ypar(ii),Zpar(ii) EndDo Mem_current = Memory(6_8*Nrecord) DEALLOCATE (Xb,Yb,Zb,VXb,VYb,VZb) Call Timing(1,1) ! start reading time end SUBROUTINE ReadDataPM !------------------------------------------------------ SUBROUTINE DENSIT !------------------------------------------------------ use Structures real*8 :: XN,YN, & X,Y,Z,D1,D2,D3,T1,T2,T3,T2W,D2W integer*8 :: IN Call Timing(2,-1) ! start timing XN =FLOAT(NGRID)+1.-1.E-7 YN =FLOAT(NGRID) Wpar = YN**3/FLOAT(Nparticles) ! Subtract mean density write(*,*) ' Init density' !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (M1,M2,M3) DO M3=1,NGRID DO M2=1,NGRID DO M1=1,NGRID FI(M1,M2,M3) = 0. END DO END DO END DO !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (IN,X,Y,Z,D1,D2,D3,T1,T2,T3,T2W,D2W) & !$OMP PRIVATE (I,J,K,I1,J1,K1) DO IN=1,Nparticles ! loop over particles X=XPAR(IN) Y=YPAR(IN) Z=ZPAR(IN) I=INT(X) J=INT(Y) K=INT(Z) D1=X-FLOAT(I) D2=Y-FLOAT(J) D3=Z-FLOAT(K) T1=1.-D1 T2=1.-D2 T3=1.-D3 T2W =T2*WPAR D2W =D2*WPAR I1=I+1 IF(I1.GT.NGRID)I1=1 J1=J+1 IF(J1.GT.NGRID)J1=1 K1=K+1 IF(K1.GT.NGRID)K1=1 !$OMP ATOMIC FI(I ,J ,K ) =FI(I ,J ,K ) +T3*T1*T2W !$OMP ATOMIC FI(I1,J ,K ) =FI(I1,J ,K ) +T3*D1*T2W !$OMP ATOMIC FI(I ,J1,K ) =FI(I ,J1,K ) +T3*T1*D2W !$OMP ATOMIC FI(I1,J1,K ) =FI(I1,J1,K ) +T3*D1*D2W !$OMP ATOMIC FI(I ,J ,K1) =FI(I ,J ,K1) +D3*T1*T2W !$OMP ATOMIC FI(I1,J ,K1) =FI(I1,J ,K1) +D3*D1*T2W !$OMP ATOMIC FI(I ,J1,K1) =FI(I ,J1,K1) +D3*T1*D2W !$OMP ATOMIC FI(I1,J1,K1) =FI(I1,J1,K1) +D3*D1*D2W ENDDO D1 = 0. ; D2 =0. !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (M1,M2,M3) REDUCTION(+:D1,D2) DO M3=1,NGRID DO M2=1,NGRID DO M1=1,NGRID D1 = D1 + FI(M1,M2,M3) D2 = D2 + FI(M1,M2,M3)**2 END DO END DO END DO D1 = D1/(float(NGRID))**3 D2 = sqrt(D2/(float(NGRID))**3) Call Timing(2,1) ! start timing write(*,'(a,2es12.4)') ' Finished density: average rms=',D1,D2 write(17,'(10x,a,es12.5,a,es12.5)') ' DM density: average=',D1,' RMS= ',D2 END SUBROUTINE DENSIT !------------------------------------------------------ SUBROUTINE PartDENS !------------------------------------------------------ use Structures real*8 :: XN,YN, & X,Y,Z,D1,D2,D3,T1,T2,T3,T2W,D2W integer*8 :: IN Real*4 :: Memory Call Timing(3,-1) ! start timing Wpar = 1. Mem_current = Memory(1_8*Nparticles) Allocate (pDens(Nparticles)) !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (IN,X,Y,Z,D1,D2,D3,T1,T2,T3,T2W,D2W) & !$OMP PRIVATE (I,J,K,I1,J1,K1) DO IN=1,Nparticles ! loop over particles X=XPAR(IN) Y=YPAR(IN) Z=ZPAR(IN) I=INT(X) J=INT(Y) K=INT(Z) D1=X-FLOAT(I) D2=Y-FLOAT(J) D3=Z-FLOAT(K) T1=1.-D1 T2=1.-D2 T3=1.-D3 T2W =T2*WPAR D2W =D2*WPAR I1=I+1 IF(I1.GT.NGRID)I1=1 J1=J+1 IF(J1.GT.NGRID)J1=1 K1=K+1 IF(K1.GT.NGRID)K1=1 pDens(IN) =FI(I ,J ,K )*T3*T1*T2W + & +FI(I1,J ,K )*T3*D1*T2W & +FI(I ,J1,K )*T3*T1*D2W & +FI(I1,J1,K )*T3*D1*D2W & +FI(I ,J ,K1)*D3*T1*T2W & +FI(I1,J ,K1)*D3*D1*T2W & +FI(I ,J1,K1)*D3*T1*D2W & +FI(I1,J1,K1)*D3*D1*D2W ENDDO D1 = 0. ; D2 =0. !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (M1,M2,M3) REDUCTION(+:D1,D2) DO IN =1,Nparticles D1 = D1 + pDens(IN) D2 = D2 + pDens(IN)**2 END DO D1 = D1/(float(Nparticles)) D2 = sqrt(D2/(float(Nparticles))) Call Timing(3,1) ! start timing write(*,'(a,2es12.4)') ' Finished particles density: average rms=',D1,D2 END SUBROUTINE PartDENS ! !------------------------------------------------------ SUBROUTINE Galaxies ! ! provives access to coordinates, velocities and density ! !------------------------------------------------------ use Structures integer*8 :: IN Call Timing(5,-1) ! start timing write(*,*) ' Start Galaxies',Nparticles !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (IN,X,Y,Z,Ux,Uy,Uz,D) DO IN=1,Nparticles ! loop over particles X = XPAR(IN) Y = YPAR(IN) Z = ZPAR(IN) Ux= VX(IN) Uy= VY(IN) Uz= VZ(IN) D = pDens(IN) ENDDO Call Timing(5,1) ! finish timing write(*,'(a,2es12.4)') ' Finished Galaxies' end SUBROUTINE Galaxies !------------------------------------------------------ SUBROUTINE Rescale !------------------------------------------------------ use Structures integer*8 :: IN Real*8 :: Dx1=0.,Dx2=0.,Dy1=0.,Dy2=0., & Dz1=0.,Dz2=0.,D1=0.,D2=0.,Dpart Call Timing(4,-1) ! start timing Xscale = Box/NGRID Vscale = 100.*Xscale/AEXPN Dpart = float(Nparticles) !write(*,*) ' Scales=',Xscale,Vscale,Dpart !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (IN) DO IN=1,Nparticles ! loop over particles XPAR(IN)=(XPAR(IN)-1.)*Xscale YPAR(IN)=(YPAR(IN)-1.)*Xscale ZPAR(IN)=(ZPAR(IN)-1.)*Xscale VX(IN) = VX(IN)*Vscale VY(IN) = VY(IN)*Vscale VZ(IN) = VZ(IN)*Vscale ENDDO !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (IN) REDUCTION(+:Dx1,Dx2,Dy1,Dy2,Dz1,Dz2,D1,D2) DO IN =1,Nparticles D1 = D1 + pDens(IN) D2 = D2 + pDens(IN)**2 Dx1 = Dx1 + VX(IN) Dx2 = Dx2 + VX(IN)**2 Dy1 = Dy1 + VY(IN) Dy2 = Dy2 + VY(IN)**2 Dz1 = Dz1 + VZ(IN) Dz2 = Dz2 + VZ(IN)**2 END DO D1 = D1/Dpart D2 = sqrt( MAX( D2/(float(Nparticles)- D1**2),1.d-14)) Dx1= Dx1/Dpart Dx2= sqrt(Dx2/Dpart) Dy1= Dy1/Dpart Dy2= sqrt(Dy2/Dpart) Dz1= Dz1/Dpart Dz2= sqrt(Dz2/Dpart) D3d= sqrt(Dx2**2 + Dy2**2 + Dz2**2) Call Timing(4,1) ! finish timing write(17,'(30x,a,5x,a)') 'average','RMS' write(17,'(10x,a,2f9.4)') ' Particle density:',D1,D2 write(17,'(10x,a,2f9.3)') ' Velocity Vx :',Dx1,Dx2 write(17,'(10x,a,2f9.3)') ' (km/s) Vy :',Dy1,Dy2 write(17,'(10x,a,2f9.3)') ' Vz :',Dz1,Dz2 write(17,'(10x,a,f9.3)') ' 3d rms Velocity km/s = ',D3d write(*,*) ' Done rescaling' END SUBROUTINE Rescale ! ---------------------------------------------------- function seconds () ! ---------------------------------------------------- ! ! purpose: returns elapsed time in seconds Integer*8, SAVE :: first=0,rate=0,i0=0 Integer*8 :: i If(first==0)Then CALL SYSTEM_CLOCK(i,rate) first =1 i0 = i seconds = 0. Else CALL SYSTEM_CLOCK(i) seconds = float(i-i0)/float(rate) EndIf end function seconds !-------------------------------------------- subroutine Timing ( ielement , isign ) !-------------------------------------------- ! 0 - total ! 1 - read, 2- density ! 3 - PartDensity, 4- rescaling cooordinates and V ! 5 - bias models, making galaxies use Structures Real*8, SAVE :: CPU(0:5) =0. CPU(ielement) = CPU(ielement) + float(isign) * seconds() If(isign == 0)Then write(*,'(a,6f10.4)') ' Time (min) Total = ',CPU(0)/60. write(*,'(a,6f10.4)') ' Read = ',CPU(1)/60. write(*,'(a,6f10.4)') ' Dens = ',CPU(2)/60. write(*,'(a,6f10.4)') ' pDens = ',CPU(3)/60. write(*,'(a,6f10.4)') ' Scale = ',CPU(4)/60. write(*,'(a,6f10.4)') ' Galax = ',CPU(5)/60. end If end subroutine Timing !-------------------------------------------- Real*4 Function Memory(i_add) !-------------------------------------------- use Structures Real*8, SAVE :: mem =0.001 ! initial memory in Gb Integer*8 :: i_add ! number of 4byte words Real*8 :: add mem = mem + i_add*4./1024.**3 Memory = mem write(*,'(a,f8.3,a)') ' Current Memory = ',Memory,'Gb' end Function Memory