! ______________________________Find particles at density maxima ! April 2008 ! A. Klypin (NMSU) ! Module Structures Integer, PARAMETER :: & NROW =1024, & ! maximum number of particles in 1D NGRID =256, & ! zero-level mesh in 1D ! must be equal to ng in ART Nmaxpart = 1024**3, & ! max number of particles for this run ! Nmaxpart must be less than NROW**3 nbyteword = 4, & ! defines length of direct-access:1or4 NPAGE = NROW**2, & ! number of particles in a record NMAX=NGRID/2, & NRECL= NPAGE*6, & ! number of words in a record Nmesh = 512, & ! size of mesh for finding maxima Nproc = 4, & ! number of processors Length = 10000 ! number of maxima per processor REAL, PARAMETER :: pi=3.14159265 REAL, ALLOCATABLE, DIMENSION(:,:,:) :: Fi,Dummy,Filter COMMON / CONTROL/ AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW, & TINTG,EKIN,EKIN1,EKIN2,AU0,AEU0, & NROWC,NGRIDC,Nspecies,Nseed,Om0,Oml0,hubble,Wp5 & ,Ocurv,extras(100) Integer, dimension(Length,Nproc) :: iLx,iLy,iLz Integer, dimension(Nproc) :: nMaxima CHARACTER*45 HEADER COMMON / ROW / XPAR(NPAGE),YPAR(NPAGE),ZPAR(NPAGE), & VX(NPAGE),VY(NPAGE),VZ(NPAGE) DIMENSION RECDAT(NRECL),wspecies(10),lspecies(10) EQUIVALENCE (RECDAT(1),XPAR(1)) & ,(wspecies(1),extras(1)), & (lspecies(1),extras(11)) Contains !--------------------------------------- ! Density Field SUBROUTINE Densit !--------------------------------------- XN =FLOAT(Nmesh)+1.-1.E-7 YN =FLOAT(Nmesh) !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (M1) DO M1=1,Nmesh FI(:,:,M1) = 0. END DO N =0 xmin =1.e+5 xmax =-1.e+5 ymin =1.e+5 ymax =-1.e+5 zmin =1.e+5 zmax =-1.e+5 ! Loop over particles If(Nspecies.eq.0)Then ! old case: 1 complete set Npages =NROW ! Number of pages N_in_last=NPAGE ! Number of particles in last page Else ! multiple masses N_particles =lspecies(Nspecies) ! Total number of particles Npages = (N_particles -1)/NPAGE +1 N_in_last=N_particles -NPAGE*(Npages-1) EndIf write (*,'(2(a,i9),a,g12.5)') ' Nspecies=',Nspecies, & ' Nparticles=',lspecies(Nspecies), & ' Last weight=',wspecies(Nspecies) Nsp_current =1 WPAR = wspecies(1) ! *8. Lspec_next=lspecies(1)+1 ff = float(NGRIDC) DO IROW=1,Npages In_page =NPAGE If(IROW.eq.Npages)In_page =N_in_last If(IROW/50*50.eq.IROW) & write (*,*) ' Read page=',IROW,' N=',N, & ' Npages=',Npages,In_page iL = NPAGE*(IROW-1) ! Loop over particles READ (21,REC=IROW) RECDAT DO IN=1,In_page N = N +1 If(N.ge.Lspec_next)Then ! next specie write (*,*) ' Another specie. Current=',Nsp_current write (*,*) ' Particle=',N,WPAR, & ' lspecies=',lspecies(Nsp_current), & ' next=',lspecies(Nsp_current+1) Nsp_current=Nsp_current+1 Lspec_next =lspecies(Nsp_current)+1 WPAR = wspecies(Nsp_current) EndIf X=(XPAR(IN)-1.)*Nmesh/ff +1. Y=(YPAR(IN)-1.)*Nmesh/ff +1. Z=(ZPAR(IN)-1.)*Nmesh/ff +1. IF(X.ge.Nmesh+1.)X=X-Nmesh IF(Y.ge.Nmesh+1.)Y=Y-Nmesh IF(Z.ge.Nmesh+1.)Z=Z-Nmesh xmin =MIN(xmin,X) ymin =MIN(ymin,Y) zmin =MIN(zmin,Z) xmax =MAX(xmax,X) ymax =MAX(ymax,Y) zmax =MAX(zmax,Z) I=INT(X) J=INT(Y) K=INT(Z) If(I.le.0)write (*,*) ' X:',X,Y,Z,I,' Irow=',IROW,IN,ifile If(J.le.0)write (*,*) ' Y:',X,Y,Z,I,' Irow=',IROW,IN,ifile If(K.le.0)write (*,*) ' Z:',X,Y,Z,I,' Irow=',IROW,IN,ifile 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.Nmesh)I1=1 J1=J+1 IF(J1.GT.Nmesh)J1=1 K1=K+1 IF(K1.GT.Nmesh)K1=1 FI(I ,J ,K ) =FI(I ,J ,K ) +T3*T1*T2W FI(I1,J ,K ) =FI(I1,J ,K ) +T3*D1*T2W FI(I ,J1,K ) =FI(I ,J1,K ) +T3*T1*D2W FI(I1,J1,K ) =FI(I1,J1,K ) +T3*D1*D2W FI(I ,J ,K1) =FI(I ,J ,K1) +D3*T1*T2W FI(I1,J ,K1) =FI(I1,J ,K1) +D3*D1*T2W FI(I ,J1,K1) =FI(I ,J1,K1) +D3*T1*D2W FI(I1,J1,K1) =FI(I1,J1,K1) +D3*D1*D2W ENDDO ENDDO write (*,*) ' Number of particles read=', N write (*,*) ' X:',xmin,xmax write (*,*) ' Y:',ymin,ymax write (*,*) ' Z:',zmin,zmax RETURN END Subroutine Densit !------------------------------------------------------------------------------------- ! Compute mean density and rms SUBROUTINE Dentes !------------------------------------------------------------------------------------ real*8 :: SUMM,SUMM2,DENM, DELR SUMM = 0. SUMM2= 0. !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (i) REDUCTION(+:SUMM,SUMM2) DO K=1,Nmesh DO J=1,Nmesh DO I=1,Nmesh SUMM = SUMM + FI(I,J,K) SUMM2= SUMM2+ FI(I,J,K)**2 ENDDO ENDDO ENDDO Total =Nmesh**3 DENM = SUMM/Total DELR = DSQRT(SUMM2/Total-DENM**2) WRITE (*,150) DELR,DENM WRITE (16,150) DELR,DENM 150 format(20('-'),' Density Distribution',20('-'),/ & ' RMS Delta Rho =',G11.4,/20x, & ' Mean Rho =',G11.4) END Subroutine Dentes !------------------------------------------------------------------------------------- ! Make map of maxima SUBROUTINE MapMaxim !------------------------------------------------------------------------------------ !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (i,j,k) DO K=1,Nmesh DO J=1,Nmesh DO I=1,Nmesh FI(I,J,K) = 0. ENDDO ENDDO ENDDO !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (i,j,k) Do iProc=1,Nproc nn = nMaxima(iProc) if(nn /=0)Then Do i=1,nn Fi(iLx(i,iProc),iLy(i,iProc),iLz(i,iProc)) =1.0 EndDo EndIf EndDo write(*,*) ' --- Done Map --' END Subroutine MapMaxim !------------------------------------------------------------------------------------- ! Make particle list SUBROUTINE ListParticles !------------------------------------------------------------------------------------ iList =0 Scale = Box/NGRID ! Loop over particles If(Nspecies.eq.0)Then ! old case: 1 complete set Npages =NROW ! Number of pages N_in_last=NPAGE ! Number of particles in last page Else ! multiple masses N_particles =lspecies(Nspecies) ! Total number of particles Npages = (N_particles -1)/NPAGE +1 N_in_last=N_particles -NPAGE*(Npages-1) EndIf write (*,'(2(a,i9),a,g12.5)') ' Nspecies=',Nspecies, & ' Nparticles=',lspecies(Nspecies), & ' Last weight=',wspecies(Nspecies) Nsp_current =1 WPAR = wspecies(1) ! *8. Lspec_next=lspecies(1)+1 ff = float(NGRIDC) N = 0 DO IROW=1,Npages In_page =NPAGE If(IROW.eq.Npages)In_page =N_in_last If(IROW/50*50.eq.IROW) & write (*,*) ' Read page=',IROW,' N=',N, & ' Npages=',Npages,iList iL = NPAGE*(IROW-1) ! Loop over particles READ (21,REC=IROW) RECDAT DO IN=1,In_page N = N +1 If(N.ge.Lspec_next)Then ! next specie write (*,*) ' Another specie. Current=',Nsp_current write (*,*) ' Particle=',N,WPAR, & ' lspecies=',lspecies(Nsp_current), & ' next=',lspecies(Nsp_current+1) Nsp_current=Nsp_current+1 Lspec_next =lspecies(Nsp_current)+1 EndIf X=(XPAR(IN)-1.)*Nmesh/ff +1. Y=(YPAR(IN)-1.)*Nmesh/ff +1. Z=(ZPAR(IN)-1.)*Nmesh/ff +1. IF(X.ge.Nmesh+1.)X=X-Nmesh IF(Y.ge.Nmesh+1.)Y=Y-Nmesh IF(Z.ge.Nmesh+1.)Z=Z-Nmesh I=INT(X) J=INT(Y) K=INT(Z) If(Fi(I,J,K)>0.5.and.Fi(I,J,K)<3.5)Then iList =iList +1 Fi(I,J,K) = Fi(I,J,K) +1. X=(XPAR(IN)-1.)*Scale Y=(YPAR(IN)-1.)*Scale Z=(ZPAR(IN)-1.)*Scale write(16,'(i11,3f10.4,3x,3i4)') N,X, Y, Z,I,J,K EndIf EndDo EndDo write(*,*) ' --- Done List Particles --:',iList END Subroutine ListParticles !------------------------------------------------------------------------------------- ! Make a test SUBROUTINE Mock !------------------------------------------------------------------------------------ write(*,*) ' pi=',pi !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (i,j,k) DO K=1,Nmesh DO J=1,Nmesh DO I=1,Nmesh ! FI(I,J,K) =1./ ((i-Nmesh/2.)**2+(j-Nmesh/2.)**2+(k-Nmesh/2.)**2+0.01) FI(I,J,K) =sin(2.*2.*pi*i/float(Nmesh))*sin(4.*2.*pi*j/float(Nmesh))* & sin(8.*2.*pi*k/float(Nmesh)) ENDDO ENDDO ENDDO !Fi =0. Fi(64,64,64) =1.0 Fi(44,64,64) =10.0 Fi(24,64,64) =100.0 write(*,*) ' --- Done Mock --' END Subroutine Mock !!------------------------------------------------------------------------------------- ! Add buffer around the density cube SUBROUTINE Expand !------------------------------------------------------------------------------------ ! expand the box periodically to make buffer layers ! of width 1 cell arounf the box Do j=1,Nmesh Do k=1,Nmesh FI(0,j,k) =FI(Nmesh,j,k) FI(j,0,k) =FI(j,Nmesh,k) FI(j,k,0) =FI(j,k,Nmesh) FI(Nmesh+1,j,k) =FI(1,j,k) FI(j,Nmesh+1,k) =FI(j,1,k) FI(j,k,Nmesh+1) =FI(j,k,1) EndDo EndDo Do j=1,Nmesh FI(0,0,j) =FI(Nmesh,Nmesh,j) FI(0,j,0) =FI(Nmesh,0,Nmesh) FI(j,0,0) =FI(0,Nmesh,Nmesh) FI(Nmesh+1,0,j) =FI(1,Nmesh,j) FI(0,j,0) =FI(Nmesh,0,Nmesh) FI(j,0,0) =FI(j,0,Nmesh) ! front panel: Y=0 FI(j,0,Nmesh+1) =FI(j,0,1) FI(0,0,j) =FI(Nmesh,0,j) FI(Nmesh+1,0,j) =FI(1,0,j) FI(j,Nmesh+1,0) =FI(j,Nmesh+1,Nmesh) ! back panel: Y=N+1 FI(j,Nmesh+1,Nmesh+1) =FI(j,Nmesh+1,1) FI(0,Nmesh+1,j) =FI(Nmesh,Nmesh+1,j) FI(Nmesh+1,Nmesh+1,j) =FI(1,Nmesh+1,j) FI(0,j,0) =FI(Nmesh,j,0) ! left low FI(Nmesh+1,j,0) =FI(1,j,0) ! right low FI(0,j,Nmesh+1) =FI(Nmesh,j,Nmesh+1) ! left up FI(Nmesh+1,j,Nmesh+1) =FI(1,j,Nmesh+1) ! right up EndDo FI(0,0,0) = FI(Nmesh,0,0) ! fill 8 corners FI(Nmesh+1,0,0) = FI(1,0,0) FI(0,0,Nmesh+1) = FI(Nmesh,0,Nmesh+1) FI(Nmesh+1,0,Nmesh+1) = FI(1,0,Nmesh+1) FI(0,Nmesh+1,0) = FI(Nmesh,Nmesh+1,0) FI(Nmesh+1,Nmesh+1,0) = FI(1,Nmesh+1,0) FI(0,Nmesh+1,Nmesh+1) = FI(Nmesh,Nmesh+1,Nmesh+1) FI(Nmesh+1,Nmesh+1,Nmesh+1) = FI(1,Nmesh+1,Nmesh+1) END Subroutine Expand !------------------------------------------------------------------------------------- ! Smooth density field with gaussian filter ! Rf is in cell units SUBROUTINE Smooth(Rf) !------------------------------------------------------------------------------------ real*8 :: SUMM,SUMM2,DENM, DELR sigmaLimit =2. ! truncation radius for filter ALLOCATE (Dummy(0:Nmesh+1,0:Nmesh+1,0:Nmesh+1),Stat=iStatus) If(iStatus.ne.0)Stop ' not enough memory for Density Field' Nlimit = INT(Rf*sigmaLimit+0.5) iOffset =128*Nmesh-1 ! offset for periodical boundary conditions ALLOCATE (Filter(-Nlimit:Nlimit,-Nlimit:Nlimit,-Nlimit:Nlimit),Stat=iStatus) If(iStatus.ne.0)Stop ' not enough memory for Filter' write(*,*) ' Nlimit=',Nlimit,Rf write(*,'(17f8.4)') Fi(50:66,64,64) write(*,'(17f8.4)') Fi(50:66,48,48) SUMM = 0. ! Set the filter !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (i,j,k) REDUCTION(+:SUMM) DO k=-Nlimit,Nlimit DO j=-Nlimit,Nlimit DO i=-Nlimit,Nlimit Filter(i,j,k) = exp(-0.5*(i**2+j**2+k**2)/Rf**2) SUMM = SUMM + Filter(i,j,k) ENDDO ENDDO ENDDO Filter =Filter/SUMM ! normalize it write(*,'(16f8.4)')(Filter(i,0-1,0),i=-Nlimit,Nlimit) write(*,'(16f8.4)')(Filter(i,0,0),i=-Nlimit,Nlimit) write(*,'(16f8.4)')(Filter(i,1,0),i=-Nlimit,Nlimit) SUMM = SUM(Filter) write(*,*) ' Sum of filter weights=',SUMM SUMM = 0. ! Smooth the density field !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (i,j,k,i1,j1,k1,i2,j2,k2,SUMM) DO k=1,Nmesh if(mod(k,50)==0)write(*,*) ' k=',k,' out of ',Nmesh DO j=1,Nmesh DO i=1,Nmesh SUMM =0. DO k1=-Nlimit,Nlimit k2 =mod(k+k1+iOffset,Nmesh)+1 DO j1=-Nlimit,Nlimit j2 =mod(j+j1+iOffset,Nmesh)+1 DO i1=-Nlimit,Nlimit i2 =mod(i+i1+iOffset,Nmesh)+1 SUMM = SUMM + Filter(i1,j1,k1)*FI(i2,j2,k2) !If(k==64.or.k==48) write(*,'(9i5,2x,g12.4)')i,j,k,i1,j1,k1,i2,j2,k2,SUMM !If(SUMM >1.e-3) write(*,'(9i5,2x,g12.4)')i,j,k,i1,j1,k1,i2,j2,k2,SUMM ENDDO ENDDO ENDDO Dummy(i,j,k) =SUMM !write(*,'(2f8.4,3x,3i5)') Fi(i,j,k),Dummy(i,j,k),i,j,k ENDDO ENDDO ENDDO !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (i,j,k) DO k=1,Nmesh DO j=1,Nmesh DO i=1,Nmesh Fi(i,j,k) = Dummy(i,j,k) ENDDO ENDDO ENDDO write(*,'(17f8.4)') Fi(50:66,64,64) write(*,'(17f8.4)') Fi(50:66,48,48) DEALLOCATE(Dummy,Filter) write(16,'(5x,a,f8.2,a,f8.2)')'Gauss filter=',Rf,' trancation=',Rf*SigmaLimit END Subroutine Smooth !------------------------------------------------------------------------------------- ! Make a list of density maxima SUBROUTINE ListMaxima !------------------------------------------------------------------------------------ nMaxima =0 iLx =0 iLy =0 iLz =0 nn =0 Chunk = (Nmesh+1.)/Nproc !!$OMP PARALLEL DO DEFAULT(SHARED) & !!$OMP PRIVATE (iProc,i,j,k,i1,i2,j1,j2,k1,k2,fmax,ff,iLast,iFirst) REDUCTION(+:nn) DO iProc =1,Nproc iFirst = (iProc-1)*Chunk+1 iLast = iProc*Chunk if(iProc==Nproc)iLast=Nmesh write(*,'(a,i4,a,6i4)') ' Proc=',iProc,' chunk=',iLast-iFirst+1,iFirst,iLast DO k=iFirst,iLast !write(*,*) ' k=',k,iProc nn = nn + 1 DO j=1,Nmesh DO i=1,Nmesh ff =FI(i,j,k) fmax =0. do k1=-1,1 do j1=-1,1 do i1=-1,1 fmax = max(FI(i+i1,j+j1,k+k1),fmax) enddo enddo enddo If(fmax .le. ff.and.ff>1.e-9)Then ! local maximum nMaxima(iProc) =nMaxima(iProc)+1 If(nMaxima(iProc)> Length)Stop ' Too many maxima. Increase Length' iLx(nMaxima(iProc),iProc) = i iLy(nMaxima(iProc),iProc) = j iLz(nMaxima(iProc),iProc) = k !write (50,'(/5x,2i7,3x,3i6,3f9.4)') iProc,nMaxima(iProc),i,j,k,ff,fmax ! If(i>10.and.i10.and.j',A45,/ & 1X,' A=',F8.3,' A0=',F8.3,' Ampl=',F8.3,' Step=',F8.3,/ & 1X,' I =',I4,' WEIGHT=',F8.3,' Ekin=',3E12.3,/ & 1X,' Nrow=',I4,' Ngrid=',I4,' Nrecl=',I9,/ & 1x,' Omega_0=',F7.3,' OmLam_0=',F7.4,' Hubble=',f7.3,/ & 1x,' Omega_curvature=',F7.3) IF(NROW.NE.NROWC) THEN WRITE (*,*) ' NROW in PARAMETER and in TAPE-FILE are different' write (*,*) ' Nrow=',NROW,' NROWC=',NROWC ENDIF IF(NGRID.NE.NGRIDC) THEN WRITE (*,*) ' NGRID in PARAMETER and in TAPE-FILE are different:' write (*,*) ' Ngrid=',NGRID,' NgridC=',NGRIDC ENDIF ! Open work file on disk write (*,*) ' start openning PMcrs0' 10 NBYTE = NRECL*4 NACCES= NBYTE / nbyteword OPEN(UNIT=21,FILE='PMcrs0.DAT',ACCESS='DIRECT', & STATUS='UNKNOWN',RECL=NACCES) Close (9) write (*,*) ' done openning PMcrd ----' RETURN 20 write (*,*) ' Error reading the header file: Abort' stop END Subroutine RDTAPE Real FUNCTION Input(Line) character(*) :: Line write(*,'(a,$)')Line read(*,*) Input End Function Input End Module Structures !-------------------------------------------------------------------------------------------------- ! Program MaximaDensity use Structures ! ! Read data and open files Open(16,file='ListMaxima.dat') CALL RDTAPE write (*,*) ' RDTAPE is done' Box = extras(100) If(Box.le.1.e-3)Then Box =INPUT(' Enter box size in Mpc/h =') EndIf ! ALLOCATE (FI(0:Nmesh+1,0:Nmesh+1,0:Nmesh+1),Stat=iStatus) If(iStatus.ne.0)Stop ' not enough memory for Density Field' !Call Mock CALL Densit ! Find density CALL Dentes ! Find rms density CALL Smooth(2.0) CALL Dentes ! Find rms density CALL Expand ! add buffer write (*,*) ' Density is ok' Call ListMaxima Call MapMaxim Call ListParticles STOP END Program MaximaDensity