!----------------------------------------------------------------- !----------------------------------------------------------------- Module Structures Real*4 :: boxsize,AEXPN,AEXP0,AMPLT,ASTEP, & Om0,Oml0,hubble,MassOne, & Box Integer*4 :: iCIC, & ! CIC (1)/ NGP(0) jstep, & ! time-step Np, & ! number of particles per record iThreads, & ! number of OMP threads Ngrid, & ! grid size NgridP ! grid for projection character*45 :: HEADER,Letter Real*4, Allocatable,Dimension(:) :: Xp,Yp,Zp, & ! read buffers VXp,VYp,Vzp, & Wp Real*4, Allocatable,Dimension(:,:) :: SurfDenX,SurfDenY,SurfDenZ ! density/potential TYPE :: RecEntry REAL*4 :: x,y,z,vxx,vyy,vzz Integer*8 :: Id0 end type RecEntry Type(RecEntry), Allocatable, Dimension(:) :: RRec end Module Structures ! !----------------------------------------------------------------- Module PMssBoundary Integer*4, parameter :: MaxDomains =2048 Integer*4, Dimension(MaxDomains) :: iLeft,iRight, jLeft,jRight,kLeft,kRight Real*4, Dimension(MaxDomains) :: xLeft,xRight, yLeft,yRight,zLeft,zRight Integer*8, Dimension(MaxDomains) :: NinNode Integer*4 :: Nodes,Nx,Ny,Nz,nRec end Module PMssBoundary !----------------------------------------------------------------- Program Project ! ! !----------------------------------------------------------------- Use Structures Real*4 :: ww Character*120 :: Line dd =seconds() ! initialize timing NgridP = 30000 jstep = 170 Allocate(SurfDenX(NgridP,NgridP)) Allocate(SurfDenY(NgridP,NgridP)) Allocate(SurfDenZ(NgridP,NgridP)) Call DENSITdm ! Read and find density for all particles in PMss CALL DENTES(DELR) ! Statistics !---------- write(17,'(f8.4,f8.2,i5,g12.5)')AEXPN,Box,NgridP,MassOne write(19,'(f8.4,f8.2,i5,g12.5)')AEXPN,Box,NgridP,MassOne write(21,'(f8.4,f8.2,i5,g12.5)')AEXPN,Box,NgridP,MassOne write(18)AEXPN,Box,NgridP,MassOne write(20)AEXPN,Box,NgridP,MassOne write(22)AEXPN,Box,NgridP,MassOne Do j=1,5 write(17,'(20g12.4)') SurfDenZ(1:10,j) write(19,'(20g12.4)') SurfDenY(1:10,j) write(21,'(20g12.4)') SurfDenX(1:10,j) EndDo Do j=1,NgridP write(18) SurfDenZ(:,j) write(20) SurfDenY(:,j) write(22) SurfDenX(:,j) EndDo close(17) close(18) STOP END Program Project !--------------------------------------- ! Density Field: ! Read files in parallel ! Each thread updates its own domain ! SUBROUTINE DENSITdm !--------------------------------------- use Structures use PMssBoundary Real*8 :: ssx,ssy,ssz Integer*8 :: iCount Character*120 :: Xi,Line write(XI,'(2a,i3.3,a,i6.6,a)')'ProjectXY','.',jstep,'.',NGRIDP,'.dat' open(17,file=TRIM(XI)) write(XI,'(2a,i3.3,a,i6.6,a)')'ProjectXYbinary','.',jstep,'.',NGRIDP,'.dat' open(18,file=TRIM(XI),form='binary') write(XI,'(2a,i3.3,a,i6.6,a)')'ProjectXZ','.',jstep,'.',NGRIDP,'.dat' open(19,file=TRIM(XI)) write(XI,'(2a,i3.3,a,i6.6,a)')'ProjectXZbinary','.',jstep,'.',NGRIDP,'.dat' open(20,file=TRIM(XI),form='binary') write(XI,'(2a,i3.3,a,i6.6,a)')'ProjectYZ','.',jstep,'.',NGRIDP,'.dat' open(21,file=TRIM(XI)) write(XI,'(2a,i3.3,a,i6.6,a)')'ProjectYZbinary','.',jstep,'.',NGRIDP,'.dat' open(22,file=TRIM(XI),form='binary') Call InitReadPMss ! read header. Sets boundaries NinNode(:) = 0 Allocate(RRec(Np)) write(*,*) ' Allocated Record: ',Np ! Initialize mesh !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (M1,M2) DO M2=1,NGRID DO M1=1,NGRID SurfDenX(M1,M2) = 0. SurfDenY(M1,M2) = 0. SurfDenZ(M1,M2) = 0. END DO END DO iCount = 0 write(*,*) ' Start reading data' Do iFile =1, Nodes Call ReadHeader(iFile) Call DENSfile(iFile,iCount) write(*,*) ' File finished =',iFile,iCount NinNode(iFile) = iCount End Do ssx = 0. ; ssy = 0. ; ssz = 0. !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (M1,M2,M3) REDUCTION(+:ssx,ssy,ssz) DO M2=1,NGRIDP DO M1=1,NGRIDP ssx = ssx + SurfDenX(M1,M2) ssy = ssy + SurfDenY(M1,M2) ssz = ssz + SurfDenZ(M1,M2) END DO END DO iCount = 0 Do iFile =1,Nodes iCount = iCount + NinNode(iFile) End Do write (17,*) ' Aexpn = ', AEXPN write (17,*) ' Box = ', Box write (17,*) ' Number of particles read = ', iCount write (17,*) ' Density/cell = ', ssx/float(NGRIDP)**2 write (17,*) ' Density/cell = ', ssy/float(NGRIDP)**2 write (17,*) ' Density/cell = ', ssz/float(NGRIDP)**2 End SUBROUTINE DENSITdm !--------------------------------------- ! Density Field: only one file SUBROUTINE DENSfile(iFile,iCount) !--------------------------------------- use Structures use PMssBoundary Logical :: inside Character*80 :: Name Logical :: exst Integer*8 :: iCount,iDummy Real*8 :: X,Y,Z,ss,W,D1,D2,D3,T1,T2,T3, & xL,xR,yL,yR,zL,zR, & T2W,D2W k = (ifile-1)/(Nx*Ny)+1 j = (ifile- (k-1)*Nx*Ny-1)/Nx +1 i = ifile- (k-1)*Nx*Ny-(j-1)*Nx ii = i +(j-1)*Nx +(k-1)*Nx*Ny if(ii/=iFile)stop ' Error in iFile -> i,j,k mapping' xL = xLeft(i) ! boundaries for this box xR = xRight(i) yL = yLeft(j) yR = yRight(j) zL = zLeft(k) zR = zRight(k) iL = iLeft(i) iR = iRight(i) jL = jLeft(j) jR = jRight(j) kL = kLeft(k) kR = kRight(k) ii = 0 jcnt = 0 W = 1. ss = NgridP/boxsize ! scale factor from Mpc to NGRID !write(*,'(a,i3,3g12.4,2i6)')' Inside DENSfile ',iFile,ss,xL,xR,iL,iR Do !---------------------- Loop over all records read(30,iostat=ierr)Nrecord !if(mod(jcnt,100)==0) & !write(*,'(a,3i8,2i10,2i10)') ' ifile=',iFile,ierr,Nrecord,ii,iCount if(ierr/= 0)exit Read(30)(RRec(m)%x,m=1,Nrecord), & (RRec(m)%y,m=1,Nrecord), & (RRec(m)%z,m=1,Nrecord) Read(30)X ! dummy reads Read(30)X ! Read(20+iOMP)(Xb(m) ,Yb(m) ,Zb(m), & ! VXb(m),VYb(m),VZb(m), & ! Ib(m),m=1,Nrecord) ii = ii + Nrecord jcnt = 0 !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (jj,x,y,z,i,j,k,D1,D2,D3,T1,T2,T3,T2W,D2W,I1,J1,K1) & !$OMP REDUCTION(+:jcnt) DO jj=1, Nrecord x = RRec(jj)%x !Xb(jj) y = RRec(jj)%y !Yb(jj) z = RRec(jj)%z !Zb(jj) If(x.ge.xL.and.x.lt.xR)Then If(y.ge.yL.and.y.lt.yR)Then If(z.ge.zL.and.z.lt.zR)Then X= ss*X +1.00001d0 Y= ss*Y +1.00001d0 Z= ss*Z +1.00001d0 I=INT(X) J=INT(Y) K=INT(Z) !---------------------------------------- CIC D1=X-FLOAT(I) D2=Y-FLOAT(J) D3=Z-FLOAT(K) T1=1.-D1 T2=1.-D2 T3=1.-D3 I1=I+1 J1=J+1 K1=K+1 If(K.ge.kL.and.K.le.kR)Then If(J.ge.jL.and.J.le.jR)Then If(I.ge.iL.and.I.le.iR)Then jcnt = jcnt +1 If(I==0.or.I>NGridP)write(*,'(a,3g12.4,i6)') 'Error X: ',x,y,z,i If(J==0.or.J>NGridP)write(*,'(a,3g12.4,i6)') 'Error Y: ',x,y,z,J If(K==0.or.K>NGridP)write(*,'(a,3g12.4,i6)') 'Error Z: ',x,y,z,K !$OMP Atomic SurfDenZ(I ,J ) =SurfDenZ(I ,J) +T1*T2 !$OMP Atomic SurfDenZ(I1,J ) =SurfDenZ(I1,J) +D1*T2 !$OMP Atomic SurfDenZ(I ,J1 ) =SurfDenZ(I ,J1) +T1*D2 !$OMP Atomic SurfDenZ(I1,J1 ) =SurfDenZ(I1,J1) +D1*D2 !$OMP Atomic SurfDenY(I ,K ) =SurfDenY(I ,K) +T1*T3 !$OMP Atomic SurfDenY(I1,K ) =SurfDenY(I1,K) +D1*T3 !$OMP Atomic SurfDenY(I ,K1 ) =SurfDenY(I ,K1) +T1*D3 !$OMP Atomic SurfDenY(I1,K1 ) =SurfDenY(I1,K1) +D1*D3 !$OMP Atomic SurfDenX(J ,K ) =SurfDenX(J ,K) +T2*T3 !$OMP Atomic SurfDenX(J1,K ) =SurfDenX(J1,K) +D2*T3 !$OMP Atomic SurfDenX(J ,K1 ) =SurfDenX(J ,K1) +T2*D3 !$OMP Atomic SurfDenX(J1,K1 ) =SurfDenX(J1,K1) +D2*D3 end If end If end If end If end If end If ENDDO iCount = iCount + jcnt END Do close(30) end SUBROUTINE DENSfile !----------------------------------------------------- ! Compute mean density and rms SUBROUTINE DENTES(DELR) !----------------------------------------------------- Use Structures real*8 :: SUM1,SUM2 SUM1 = 0. SUM2= 0. Nn = 0 Am = 0. !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (j,i) & !$OMP REDUCTION(+:SUM1,SUM2) DO J=1,NGRIDP DO I=1,NGRIDP SUM1 = SUM1 + SurfDenX(I,J) SUM2= SUM2 + SurfDenX(I,J)**2 ENDDO ENDDO Total =(FLOAT(NGRIDP))**2 DENM = SUM1/Total DELR = DSQRT(SUM2/Total-DENM**2) !DensAv=PARTW*float(NROW)**3/(float(NGRID))**3 WRITE (*,150) DELR,DENM,SUM1!,DensAv WRITE (17,150) DELR,DENM,SUM1!,DensAv 150 format(20x,'Density is in units partle number per cell', & /20x, & ' RMS =',G11.4,/20x, & ' Mean =',G13.5,' Sum =',G13.5) END SUBROUTINE DENTES ! ! ! ---------------------------------------------------- real function seconds () ! ---------------------------------------------------- ! ! purpose: returns elapsed time in seconds Integer, SAVE :: first=0,rate=0 Real , SAVE :: t0=0. !real*8 dummy !real tarray(2) If(first==0)Then CALL SYSTEM_CLOCK(i,rate) first =1 t0 =float(i)/float(rate) seconds = 0. Else ! seconds = etime(tarray) CALL SYSTEM_CLOCK(i) seconds = float(i)/float(rate)-t0 EndIf end function seconds !------------------------------------------- ! Routine called by DENSITdm ! - Initializes gadget PMss reading ! - Sets boundaries for each file/domain SUBROUTINE InitReadPMss Use Structures Use PMssBoundary Character*80 :: Name Logical :: exst !--- read info from first file to check if everything is ok node = 1 write(Name,'(a,i4.4,a,i4.4,a)') 'PMss/PMss.',jstep,'.',node,'.DAT' Inquire(file=TRIM(Name),exist=exst) If(.not.exst)Then write(*,'(2a)') ' Attempt to read: ',TRIM(Name) Stop ' File not found' End If open(1,file=Trim(Name),form='unformatted') Read(1)HEADER Read(1)AEXPN,ASTEP,ISTEP,NROWC,NGRIDC,Nspecies, & Nseed,Om0,Oml0,hubble,Box,MassOne write(*,*) ' A =',AEXPN write(*,'(a,1p,5g12.4)') ' Omegas =',Om0,Oml0,hubble Read(1) i_node,Nx,Ny,Nz,dBuffer ,nBuffer Read(1) xL,xR,yL,yR,zL,zR Read(1) Np Read(1) Nrecord backspace (1) Ngrid = NGRIDC Nodes = Nx*Ny*Nz nRec = Nrecord Np = Nrecord write (*,'(a,f8.4)') '# Aexp = ',AEXPN write (*,'(a,f8.4)') '# Om0 = ',Om0 write (*,'(a,f9.4)') '# Box (Mpch) = ',Box write (*,'(a,i5)') '# Number of files = ',Nodes write (*,'(a,f8.4)') '# Buffer (Mpch) = ',dBuffer write (*,'(a,1p,g12.4)') '# Mass One (Msunh) = ',MassOne write (*,'(a,f12.6)')'# Npartcle/file 1e6 = ',float(Np)/1.e6 write (*,'(a,i12)')'# Npartcle/record = ',nRec boxsize = Box close(1) dx = (NgridP+1.)/Nx ! Set boundaries of domains dy = (NgridP+1.)/Ny dz = (NgridP+1.)/Nz qx = Box/(NgridP) write(*,*) ' Step of domains qx =',qx write(*,*) ' Left/Right X 1st =',xL,xR Do i=1,Nx iLeft(i) =(i-1)*dx +1 iRight(i) =MIN(INT(i*dx),NGRIDP) xLeft(i) = (iLeft(i) -2)*qx xRight(i) = (iRight(i))*qx end Do Do i=1,Ny jLeft(i) =(i-1)*dy +1 jRight(i) =MIN(INT(i*dy),NGRIDP) yLeft(i) = (jLeft(i) -2)*qx yRight(i) = (jRight(i))*qx print *,i,yLeft(i),yRight(i) end Do Do i=1,Nz kLeft(i) =(i-1)*dz +1 kRight(i) =MIN(INT(i*dz),NGRIDP) zLeft(i) = (kLeft(i) -2)*qx zRight(i) = (kRight(i))*qx end Do write(*,'(a,25i9)')'X Left = ',(iLeft(i), i=1,Nx) write(*,'(a,25i9)')' Right= ',(iRight(i),i=1,Nx) write(*,'(a,25f9.3)')' Left = ',(xLeft(i), i=1,Nx) write(*,'(a,25f9.3)')' Right= ',(xRight(i),i=1,Nx) End SUBROUTINE InitReadPMss !------------------------------------------- ! Routine called by DENSITdm ! ! SUBROUTINE ReadHeader(iFile) Use Structures Use PMssBoundary Character*120 :: Name Logical :: exst !--- read info from first file to check if everything is ok write(Name,'(a,i4.4,a,i4.4,a)') 'PMss/PMss.',jstep,'.',iFile,'.DAT' Inquire(file=TRIM(Name),exist=exst) If(.not.exst)Then write(*,'(2a)') ' Attempt to read: ',TRIM(Name) Stop ' File not found' End If Nn = 0 open(30,file=Trim(Name),form='unformatted') Read(30)HEADER Read(30) x Read(30) i Read(30) x Read(30) Nn write(*,'(a,i6,a,i10,3x,3i6)') ' Init finished for iFile=',iFile, & ' #particles=',Nn,iFile end SUBROUTINE ReadHeader