!-------------------------------------------------------------------------------- ! Split simulation box in smaller chunks ! for BDM ! (nx,ny,nz) - number of boxes in each direction ! Module SetArrs Integer, PARAMETER :: & NROW = 4096, & NGRID = 256, & nx = 4, & !number of regions in X ny = 4, & nz = 4, & Nfiles = nx*ny*nz, & ! number of files Nrecord = 500000, & ! max number of particles in buffer nbyteword = 4, & NPAGE = NROW**2, & ! # particles in a record NRECL = NPAGE*6 Integer*8,PARAMETER ::Np=2048_8**3 ! max number of particles COMMON / ROW / XPAR(NPAGE),YPAR(NPAGE),ZPAR(NPAGE), & VX(NPAGE),VY(NPAGE),VZ(NPAGE) COMMON / CONTROL/ AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW, & TINTG,EK,EK1,EK2,AU0,AEU0, & NROWC,NGRIDC,Nspecies,Nseed,Om0,Oml0,hubble,Wp5, & Ocurv,extras(100) COMMON / HEADDR/ HEADER CHARACTER*45 HEADER Real*4, PARAMETER :: Rhalo=2.2 Real, DIMENSION(Np) :: Xp,Yp,Zp,VXp,VYp,Vzp Real*4, DIMENSION(Nrecord,Nfiles) :: xb,yb,zb,vxb,vyb,vzb Integer*8, DIMENSION(Nrecord,Nfiles) :: id_part Character*80 :: fNames Real :: Box,RECDAT(NRECL),wspecies(10) Integer*8 :: lspecies(10) Integer*8 :: jCount(Nfiles),iCount(Nfiles) EQUIVALENCE (wspecies(1),extras(1)), & (lspecies(1),extras(11)) !$OMP THREADPRIVATE(/ROW/) Contains !-------------------------------------------------------------------------------------------- ! ! Open Files, read control information ! Subroutine OpenPM Character :: fname*120 Open(4,file ='PMcrd.DAT',form ='UNFORMATTED',status ='UNKNOWN') READ (4,err=20,end=20) HEADER, & AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW, & TINTG,EKIN,EKIN1,EKIN2,AU0,AEU0, & NROWC,NGRIDC,Nspecies,Nseed,Om0,Oml0,hubble,Wp5 & ,Ocurv,extras WRITE (*,100) HEADER, & AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW, & EKIN,EKIN1,EKIN2, & NROWC,NGRID,NRECL,Om0,Oml0,hubble, & Ocurv 100 FORMAT(1X,'Header=>',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,NGRID (PMparamters.h)=',NROW,NGRID write (*,*) ' NROW,NGRID (PMcrd.DAT) =',NROWC,NGRIDC ENDIF IF(NGRID.NE.NGRIDC) THEN WRITE (*,*) & ' NGRID in PARAMETER and in TAPE-FILE are different:' write (*,*) ' Ngrid=',NGRID,' NgridC=',NGRIDC ENDIF write(*,*) ' Number of particles =',lspecies(Nspecies) Box = extras(100) CLOSE (4) NACCES= NRECL*4 / nbyteword Jfiles =1 If(lspecies(Nspecies)> 1024**3)Jfiles=8 Do i =1,Jfiles write(fname,'(a,i1,a)')'PMcrs',i-1,'.DAT' OPEN(UNIT=20+i,FILE=TRIM(fname),ACCESS='DIRECT', & FORM='unformatted',STATUS='UNKNOWN',RECL=NACCES) EndDo RETURN 20 write (*,*) ' Error reading the header file: Abort' stop end Subroutine OpenPM end Module SetArrs !-------------------------------------------------------------------------------- ! Program SplitBox use SetArrs Character :: Line*80 integer*8 :: Lsmall,N_particles,JPAGE,iL,IN,N dR =Rhalo ! boundary width in Mpch units Call OpenPM Xscale= Box/Ngrid ! Scale for comoving coordinates dR = dR/Xscale Call ReadPntP(N) Call Split(N,dR) Stop end Program SplitBox !--------------------------------------------------------- ! Read particles SUBROUTINE Split(N,dR) !-------------------------------------------------------------- use SetArrs Integer*8 :: N,i,N_particles,Lsmall,Jpage,iL,IN Real*8 :: xmin,xmax,ymin,ymax,zmin,zmax Real*4 :: bndry(2,3,Nfiles) dx = float(NGRID)/float(nx) dy = float(NGRID)/float(ny) dz = float(NGRID)/float(nz) ! define boundaries of regions kk =0 Do k=1,nz bz1 = 1.-dR+ (k-1.)*dz bz2 = 1.+dR+ k *dz Do j=1,ny by1 = 1.-dR+ (j-1.)*dy by2 = 1.+dR+ j *dy Do i=1,nx bx1 = 1.-dR+ (i-1.)*dx bx2 = 1.+dR+ i *dx kk =kk +1 bndry(1,1,kk) =bx1 bndry(2,1,kk) =bx2 bndry(1,2,kk) =by1 bndry(2,2,kk) =by2 bndry(1,3,kk) =bz1 bndry(2,3,kk) =bz2 EndDo EndDo EndDo If(kk/= Nfiles)Stop ' Error in Nfiles' write(*,*) ' File x y z' Do k=1,Nfiles write(*,'(i3,5x,3(2x,2f8.2))')k,((bndry(i,j,k),i=1,2),j=1,3) EndDo ! Open files and write control information Do k=1,Nfiles write(fNames,'(a,i4.4,a,i4.4,a)') 'PMs.',k,'.',ISTEP,'.DAT' open(100+k,file=fNames,form='unformatted') write(100+k)HEADER write(100+k)AEXPN,ASTEP,ISTEP,NROWC,NGRIDC,Nspecies, & Nseed,Om0,Oml0,hubble,Box write(100+k) k,Nx,Ny,Nz write(100+k) ((bndry(i,j,k),i=1,2),j=1,3) ! write(100+k,*)HEADER ! write(100+k,100)AEXPN,ASTEP,ISTEP,NROWC,NGRIDC,Nspecies, & ! Nseed,Om0,Oml0,hubble,Box ! write(100+k,'(4i5)') k,Nx,Ny,Nz ! write(100+k,'(6f8.3)') ((bndry(i,j,k),i=1,2),j=1,3) EndDo 100 format(f9.4,f9.5,i6,2i6,i4,i10,4f8.3) jCount = 0 ! initialize counters: vector assignement iCount = 0 ! Find particles for each region. ! Once buffer is full, dump to a file !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (i,xx,yy,zz,i2,j2,k2,ii,x2,y2,z2,ind_x,ind_y,ind_z) Do i=1,N xx = Xp(i) yy = Yp(i) zz = Zp(i) if(mod(i,1000000)==0)write(*,*) 'Npart=',i ii = 0 i2 = 0 j2 = 0 k2 = 0 x2 =0.; y2 =0.; z2 =0. ind_x = INT((xx-1.)/ dx)+1 ind_y = INT((yy-1.)/ dy)+1 ind_z = INT((zz-1.)/ dz)+1 If((ind_x-1)*dx+1.+dR > xx)Then ! X left boundary ii = ii+1 ! number of other boxes to put the particle in i2 = ind_x -1 ; x2 = xx If(i2<1)Then ! periodical boundary conditions i2 = nx ; x2 =x2 + NGRID endIf Else If(ind_x*dx+1.-dR < xx)Then ! right boundary ii = ii+1 ! number of other boxes to put the particle in i2 = ind_x +1 ; x2 = xx If(i2>nx)Then ! periodical boundary conditions i2 = 1 ; x2 =x2 - NGRID endIf EndIf If((ind_y-1)*dy+1.+dR > yy)Then ! Y left boundary ii = ii+1 ! number of other boxes to put the particle in j2 = ind_y -1 ; y2 = y If(j2<1)Then ! periodical boundary conditions j2 = ny ; y2 =y2 + NGRID endIf Else If(ind_y*dy+1.-dR < yy)Then ! right boundary ii = ii+1 ! number of other boxes to put the particle in j2 = ind_y +1 ; y2 = yy If(j2> ny)Then ! periodical boundary conditions j2 = 1 ; y2 =y2 - NGRID endIf EndIf If((ind_z-1)*dz+1.+dR > zz)Then ! Z left boundary ii = ii+1 ! number of other boxes to put the particle in k2 = ind_z -1 ; z2 = zz If(k2<1)Then ! periodical boundary conditions k2 = nz ; z2 =z2 + NGRID endIf Else If(ind_z*dz+1.-dR < zz)Then ! right boundary ii = ii+1 ! number of other boxes to put the particle in k2 = ind_z +1 ; z2 = zz If(k2> nz)Then ! periodical boundary conditions k2 = 1 ; z2 =z2 - NGRID endIf EndIf ! Store results in buffers. ! four possible combinations: ! ii = 0 - particle is only in one region ! = 1 - in two regions ! = 2 - in four regions ! = 3 - in eight regions ! iFile = ind_x + (ind_y-1)*nx +(ind_z-1)*nx*ny ! write(*,'(/i10,3x,3f8.3,3x,3i4,3x,4i3,3x,3f8.3)') i,xx,yy,zz,ind_x,ind_y,ind_z,ii,i2,j2,k2,x2,y2,z2 !$OMP critical CALL BufferUpdate(ind_x,ind_y,ind_z,xx,yy,zz,i) !$OMP end critical If(ii == 1 ) Then ! two boxes If(i2 .ne. 0)Then !$OMP critical CALL BufferUpdate(i2,ind_y,ind_z,x2,yy,zz,i) !$OMP end critical Else If(j2.ne.0)Then !$OMP critical CALL BufferUpdate(ind_x,j2,ind_z,xx,y2,zz,i) !$OMP end critical Else !$OMP critical CALL BufferUpdate(ind_x,ind_y,k2,xx,yy,z2,i) !$OMP end critical EndIf ! i2 /=2 Else if(ii == 2 ) Then ! four boxes If(i2 == 0)Then !$OMP critical CALL BufferUpdate(ind_x,j2,ind_z, xx,y2,zz,i) CALL BufferUpdate(ind_x,ind_y,k2,xx,yy,z2,i) CALL BufferUpdate(ind_x,j2,k2, xx,y2,z2,i) !$OMP end critical Else If(j2 == 0)Then !$OMP critical CALL BufferUpdate(i2,ind_y,ind_z, x2,yy,zz,i) CALL BufferUpdate(ind_x,ind_y,k2,xx,yy,z2,i) CALL BufferUpdate(i2,ind_y,k2, x2,yy,z2,i) !$OMP end critical Else !$OMP critical CALL BufferUpdate(i2,ind_y,ind_z,x2,yy,zz,i) CALL BufferUpdate(ind_x,j2,ind_z,xx,y2,zz,i) CALL BufferUpdate(i2,j2,ind_z, x2,y2,zz,i) !$OMP end critical EndIf Else If(ii==3)Then ! ii == 3 eight boxes !!$OMP critical CALL BufferUpdate(i2,ind_y,ind_z, x2,yy,zz,i) CALL BufferUpdate(ind_x,j2,ind_z, xx,y2,zz,i) CALL BufferUpdate(i2,j2,ind_z, x2,y2,zz,i) CALL BufferUpdate(ind_x,ind_y,k2,xx,yy,z2,i) CALL BufferUpdate(ind_x,j2,k2, xx,y2,z2,i) CALL BufferUpdate(i2,ind_y,k2, x2,yy,z2,i) CALL BufferUpdate(i2,j2,k2, x2,y2,z2,i) !!$OMP end critical EndIf EndDo Do i=1,Nfiles If(iCount(i).ne.0)Then NN = iCount(i) write(100+iFile,'(i10)')NN write(100+iFile,200)(xb(jj,iFile),jj=1,NN) write(100+iFile,200)(yb(jj,iFile),jj=1,NN) write(100+iFile,200)(zb(jj,iFile),jj=1,NN) write(100+iFile,200)(vxb(jj,iFile),jj=1,NN) write(100+iFile,200)(vyb(jj,iFile),jj=1,NN) write(100+iFile,200)(vzb(jj,iFile),jj=1,NN) write(100+iFile,210)(id_part(jj,iFile),jj=1,NN) EndIf EndDo 200 format(16g12.4) 210 format(16i11) end SUBROUTINE Split !--------------------------------------------------------- ! Update buffers and write files SUBROUTINE BufferUpdate(i3,j3,k3,xx,yy,zz,i) !-------------------------------------------------------------- use SetArrs Integer*8 :: N,i,N_particles,Lsmall,Jpage,iL,IN iFile = i3 + (j3-1)*nx +(k3-1)*nx*ny jCount(iFile) =jCount(iFile) +1 iCount(iFile) =iCount(iFile) +1 xb(jCount(iFile) ,iFile) = xx yb(jCount(iFile) ,iFile) = yy zb(jCount(iFile) ,iFile) = zz vxb(jCount(iFile) ,iFile) = VXp(i) vyb(jCount(iFile) ,iFile) = VYp(i) vzb(jCount(iFile) ,iFile) = VZp(i) id_part(jCount(iFile) ,iFile) = i ! write(*,'(3x,i10,3f9.4,i10)')iFile,xb(jCount(iFile) ,iFile),yb(jCount(iFile) ,iFile), & ! zb(jCount(iFile) ,iFile), id_part(jCount(iFile) ,iFile) If(jCount(iFile) == Nrecord)Then write(100+iFile)Nrecord write(100+iFile)(xb(jj,iFile),jj=1,Nrecord) write(100+iFile)(yb(jj,iFile),jj=1,Nrecord) write(100+iFile)(zb(jj,iFile),jj=1,Nrecord) write(100+iFile)(vxb(jj,iFile),jj=1,Nrecord) write(100+iFile)(vyb(jj,iFile),jj=1,Nrecord) write(100+iFile)(vzb(jj,iFile),jj=1,Nrecord) write(100+iFile)(id_part(jj,iFile),jj=1,Nrecord) ! write(100+iFile,'(/a,i10)')' Npart=',Nrecord ! write(100+iFile,200)(xb(jj,iFile),jj=1,Nrecord) ! write(100+iFile,200)(yb(jj,iFile),jj=1,Nrecord) ! write(100+iFile,200)(zb(jj,iFile),jj=1,Nrecord) ! write(100+iFile,200)(vxb(jj,iFile),jj=1,Nrecord) ! write(100+iFile,200)(vyb(jj,iFile),jj=1,Nrecord) ! write(100+iFile,200)(vzb(jj,iFile),jj=1,Nrecord) ! write(100+iFile,210)(id_part(jj,iFile),jj=1,Nrecord) jCount(iFile) =0 EndIf 200 format(16g12.4) 210 format(16i11) end SUBROUTINE BufferUpdate !--------------------------------------------------------- ! Read particles SUBROUTINE ReadPntP(N) !-------------------------------------------------------------- use SetArrs Integer*8 :: N,i,N_particles,Lsmall,Jpage,iL,IN Real*8 :: xmin,xmax,ymin,ymax,zmin,zmax,xbig ! lspecies(1) = lspecies(1)/2_8 !!! ! write (*,*) ' Temp: N=',lspecies(1) N =0 Jpage = NPAGE xmin =1.e10 xmax =-1.e10 ymin =1.e10 ymax =-1.e10 zmin =1.e10 zmax =-1.e10 xzero = 1. xbig = float(NGRID)+1.d0-1.d-4 write (*,*) ' X small/big=',xzero,xbig Nsmall =0 Lsmall =lspecies(1) N_particles =lspecies(Nspecies) ! Total number of particles Npages = (N_particles -1)/Jpage +1 N_in_last=N_particles -Jpage*(Npages-1) write (*,*) ' Pages=',Npages,' Species=',Nspecies write (*,*) ' N_in_last=',N_in_last !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (IROW,In_page,iL,IN,j) DO IROW=1, Npages ! Loop over particle pages In_page =NPAGE If(IROW.eq.Npages)In_page =N_in_last iL = JPAGE*(IROW-1) if(mod(IROW,10).eq.0)write (*,*)' Read page=',IROW,' Np=',iL ! Loop over particles If(lspecies(Nspecies)>1024**3)Then j = (IROW-1)/64 +1 if(j>8.or.j<1)write(*,*) ' Error in file number: ',j read(20+j,REC=mod(IROW-1,64)+1)xpar,ypar,zpar,vx,vy,vz Else READ (21,REC=IROW) XPAR,YPAR,ZPAR,VX,VY,VZ EndIf DO IN=1,In_page Xp(iL+IN) =XPAR(IN) Yp(iL+IN) =YPAR(IN) Zp(iL+IN) =ZPAR(IN) VXp(iL+IN) =VX(IN) VYp(iL+IN) =VY(IN) VZp(iL+IN) =VZ(IN) EndDo ! IN EndDo ! IROW N = lspecies(Nspecies) write(*,*) ' Done Reading Points. N=',N !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (i) & !$OMP REDUCTION(MAX:xmax,ymax,zmax) & !$OMP REDUCTION(MIN:xmin,ymin,zmin) Do i =1,N If(Xp(i).gt.xbig)Xp(i)=Xp(i)-1.e-4 If(Yp(i).gt.xbig)Yp(i)=Yp(i)-1.e-4 If(Zp(i).gt.xbig)Zp(i)=Zp(i)-1.e-4 If(i.le.Lsmall)Then xmin = min(Xp(i),xmin) xmax =max(Xp(i),xmax) ymin = min(Yp(i),ymin) ymax =max(Yp(i),ymax) zmin = min(Zp(i),zmin) zmax =max(Zp(i),zmax) EndIf ENDDO write(*,*) ' Small Particles:',Lsmall write(*,*) ' x:',xmin,xmax write(*,*) ' y:',ymin,ymax write(*,*) ' z:',zmin,zmax write (*,*) ' Nmaxpart=',Nmaxpart,Np !do i=1,1000 !write (*,'(i10," Coords=",6g13.6)')i,Xp(i),Yp(i),Zp(i), !& VXp(i),VYp(i),VZp(i) !EndDo Return End SUBROUTINE ReadPntP