! _______________________ START 3-D SIMULATIONS ! ! Klypin, August 1993, March 2010 -MPI ! extras: 1-10 = wspecies ! 11 -32 = lspecies (8bytes/number, 10 numbers) ! 40 = NPAGE ! 41 = Npages = number of pages in one file ! 80 = luxury switch (=1 to use luxury) ! 81 = level of integers for luxury ! 100 = Box size in Mpc !------------------------------------------------------------- MODULE Randoma Contains !-------------------------------------- ! normal random numbers FUNCTION GAUSS (M) !-------------------------------------- X = 0. DO I = 1, 5 X = X + RANDd (M) EndDo X2 = 1.5491933 * (X - 2.5) GAUSS = X2 * (1. - 0.01 * (3. - X2 * X2) ) RETURN END FUNCTION GAUSS !------------------------------------------------ ! random number generator FUNCTION RANDd (M) !------------------------------------------------ DATA LC, AM, KI, K1, K2, K3, K4, L1, L2, L3, L4 / 453815927, & 2147483648., 2147483647, 536870912, 131072, 256, 16777216, 4, & 16384, 8388608, 128 / ML = M / K1 * K1 M1 = (M - ML) * L1 ML = M / K2 * K2 M2 = (M - ML) * L2 ML = M / K3 * K3 M3 = (M - ML) * L3 ML = M / K4 * K4 M4 = (M - ML) * L4 M5 = KI - M IF (M1.GE.M5) M1 = M1 - KI - 1 ML = M + M1 M5 = KI - ML IF (M2.GE.M5) M2 = M2 - KI - 1 ML = ML + M2 M5 = KI - ML IF (M3.GE.M5) M3 = M3 - KI - 1 ML = ML + M3 M5 = KI - ML IF (M4.GE.M5) M4 = M4 - KI - 1 ML = ML + M4 M5 = KI - ML IF (LC.GE.M5) ML = ML - KI - 1 M = ML + LC RANDd = M / AM RETURN END FUNCTION RANDd end module Randoma ! ! !------------------------------------------------------------------ MODULE setInitialConditions PUBLIC INCLUDE 'PMinitialp.h' INTEGER, PARAMETER :: NPAGE = min(1024**2,NROW**2), & ! # particles in a record Npages = 64, & ! maximum # of pages per file PMcrsXXX Nslaves =32, & ! # of mpi tasks for IO NMAX = NGRID/2, & NRECL = NPAGE*6, & ! # particles in a record NARR = MAX(2*NROW+1,NGRID+1), & ! FFT NF67 = MAX(NROW,NGRID/2), & Nblocks=NROW/Lblock, & NSPEC = NROW/2, & ! No waves shorter than Ny marr = NROW+1, & ! Arrays for FFT mf67 = NROW/2 REAL, DIMENSION(10) :: timing =0. REAL :: alpha =0. 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) COMMON / HEADDR/ HEADER CHARACTER*45 HEADER COMMON /FOURAR/Zf(NARR),Yf(NARR) COMMON /F67COM/ & IBC, IP, ISL, L1, N2, & N3, N4, N7, & SI(NF67), INDX(NF67) Real*4, ALLOCATABLE, dimension(:) :: XPAR,YPAR,ZPAR,& VX,VY,VZ,YPt,ZPt,VYt,VZt Real*4, ALLOCATABLE, dimension(:,:,:) :: Buff Real*4 :: wspecies(10) Integer*8, DIMENSION(10) :: lspecies(10), nLevel (10) Integer*4 :: rank EQUIVALENCE & (wspecies(1),extras(1)), & (lspecies(1),extras(11)) integer*4 :: iMask(Nblocks,Nblocks) ! Refinement Mask integer*8 :: Allspecies(10,NROW) REAL, ALLOCATABLE, DIMENSION (:,:) :: GR,GRbuff Real*4, DIMENSION(NPAGE) :: XPt,VXt INTEGER*4, DIMENSION(NPAGE) :: MPt COMMON /SeedsPage/ NseedP(NROW) COMMON / LevelsM/ Max_lev_abs,Nmin_lev,Nmax_lev Contains !-------------------------------------------------- ! : Store and retrieve seeds for parallelization ! of random number generator !-------------------------------------------------- SUBROUTINE SetRandom(iCont) use Randoma use mpi INCLUDE 'luxuryp.h' Character (LEN=70) :: FileName Logical :: FileExists =.false., RNGluxury =.false. Logical :: iCont REAL :: t0=0,t1=0. ! counters for timing CALL CPU_TIME(t0) ! initialize time counter If(abs(extras(80)-1.) < 1.e-5)RNGluxury =.true. If(rank == 0)Then write(FileName,'(a,i1.1,".",i4.4,".",i8.8,".dat")') & 'Seeds.',INT(extras(80)),NROW,Nseed INQUIRE(file=FileName,EXIST=FileExists) If(FileExists)Then ! read the data write(*,'(10x,2a)')'Open Existing FileName=',FileName Open(62,file=FileName,form='unformatted') read(62)NROWdummy,Nseeddummy if(NROWdummy /= NROW .or. Nseeddummy /= Nseed)Then write(*,*) ' File with seeds does not match parameters' write(*,*) ' Read Nrow =',NROWdummy,' Internal=',NROW write(*,*) ' Read Nseed=',Nseeddummy,' Internal=',Nseed STOP endIf If(RNGluxury)Then ! read(62)i24A,j24A,in24A,kountA,carryA,gsetA,iFlagA read(62)SeedsPage Close(62) Else read(62)i24A EndIf Close(62) Else ! open file with seeds write(*,'(10x,2a)')'Create FileName=',FileName Open(62,file=FileName,form='unformatted') ! generate seeds If(RNGluxury)Then ! Detect random generator: luxury or nrandd iFlag = 0 gSet = 0. Ns = Nseed ; lux = 2 Do k=1,NROW ! luxury if(k/10*10==k)write(*,*) ' page number=',k,iFlag i24A(k) = i24 j24A(k) = j24 in24A(k) = in24 kountA(k) = kount carryA(k) = carry gsetA(k) = gSet iFlagA(k) = iFlag Do m=1,24 SeedsPage(m,k) = seeds(m) EndDo If(NROW>4096)Then dummy = RANDd(Ns) CALL rluxgo (lux, Ns, 0, 0) ! initialize luxury Else Do j=1,NROW Do i=1,NROW ! If(i+j+k.ne.3)Then x = GAUSS3fake(gSet,iFlag) ! EndIf EndDo EndDo EndIf EndDo write(62)NROW,Nseed write(62)i24A,j24A,in24A,kountA,carryA,gsetA,iFlagA write(62)SeedsPage Close(62) Else write(*,*) ' Nseed=',Nseed NRAND = Nseed Do k=1,NROW ! gauss and randd if(k/10*10==k)write(*,*) ' page number=',k i24A(k) = NRAND Do j=1,NROW Do i=1,NROW If(i+j+k.ne.3)Then x = GAUSS(NRAND) EndIf EndDo EndDo EndDo write(62)NROW,Nseed write(62)i24A Close(62) EndIf ! end luxury/randd switch write(*,'(/2(20("-"),a/))')'File with seeds was created. STOP', & 'Do not change anything. Just restart the code to produce Initial Conditions' iCont = .false. EndIf ! file exists switch CALL CPU_TIME(t1) timing(2) = t1-t0 ! time for setting random numbers EndIf ! rank == 0 CALL MPI_Barrier(MPI_COMM_WORLD,ierr) CALL MPI_Bcast(RNGluxury,1,MPI_LOGICAL,0,MPI_COMM_WORLD,ierr) CALL MPI_Bcast(i24A,NROW,MPI_INTEGER,0,MPI_COMM_WORLD,ierr) If(RNGluxury)Then ! CALL MPI_Bcast(j24A,NROW,MPI_INTEGER,0,MPI_COMM_WORLD,ierr) CALL MPI_Bcast(in24A,NROW,MPI_INTEGER,0,MPI_COMM_WORLD,ierr) CALL MPI_Bcast(kountA,NROW,MPI_INTEGER,0,MPI_COMM_WORLD,ierr) CALL MPI_Bcast(iFlagA,NROW,MPI_INTEGER,0,MPI_COMM_WORLD,ierr) CALL MPI_Bcast(carryA,NROW,MPI_REAL,0,MPI_COMM_WORLD,ierr) CALL MPI_Bcast(gsetA,NROW,MPI_REAL,0,MPI_COMM_WORLD,ierr) CALL MPI_Bcast(iFlagA,NROW,MPI_INTEGER,0,MPI_COMM_WORLD,ierr) CALL MPI_Bcast(SeedsPage,24*NROW,MPI_REAL,0,MPI_COMM_WORLD,ierr) EndIf If(rank /= 0)Then If(RNGluxury)Then ! Detect random generator: luxury or nrandd Ns = Nseed lux = extras(81) CALL rluxgo (lux, Nseed, 0, 0) ! initialize luxury EndIf EndIf Return End SUBROUTINE SetRandom !-------------------------------------------------- ! sqrt(Power spectrum) ! k = (2pi/L) !-------------------------------------------------- FUNCTION TRUNF (WK) !------------------------------------------------- PARAMETER (NtabM = 100000) COMMON / PTABLE / xkt (0:NtabM), Pkt (0:NtabM), Ntab, StepK, & alog0, iFlag COMMON / TRUNCOM / Om, Omb, Omc, Omnu, Par (6), ns, qqscaleb, & QSCALE Real ns, k IF (WK.GE.FLOAT (NSPEC) ) THEN TRUNF = 0. RETURN ENDIF k = QSCALE * wk ! use approximations If (iFlag.gt.0) Then ! Holtzman approx If (Par (6) .ne.0.) Then sk = sqrt (k) TRUNF = k** (ns / 2.) / (1. + sk * (Par (2) + sk * (Par (3) & + sk * (Par (4) + sk * (Par (5) ) ) ) ) ) ** (Par (6) ) ! BBKS + Sugiyama approx Else ! Gammaeff =Om*hsmall/exp(Omb*(1.+sqrt(hsmall/0.5)/Om)) ! Q = wk/hsmall/Gammaeff Q = k * qqscaleb TRUNF = k** (ns / 2.) * LOG (1. + Par (1) * Q) / (Par (1) * Q) & / sqrt (sqrt (1. + Q * (Par (2) + Q * (Par (3) **2 + Q * (Par (4) & **3 + Q * (Par (5) **4) ) ) ) ) ) EndIf ! use table for P(k) Else TRUNF = sqrt (Ppk (k) ) EndIf RETURN END FUNCTION TRUNF !--------------------------------------- FUNCTION Ppk (x) ! interpolate table with p(k) ! x is in real 1/Mpc !--------------------------------------- PARAMETER (NtabM = 100000) COMMON / PTABLE / xkt (0:NtabM), Pkt (0:NtabM), Ntab, StepK, & alog0, iFlag ! slope is ns =-3 If (x.ge.xkt (Ntab) ) Then Ppk = Pkt (Ntab) / (x / xkt (Ntab) ) **3 Return EndIf If (x.lt.xkt (1) ) Then ! slope is ns=1 Ppk = Pkt (1) * (x / xkt (1) ) Return EndIf ind = INT ( (log10 (x) - alog0) / StepK) + 1 dk = xkt (ind+1) - xkt (ind) Ppk = (Pkt (ind) * (xkt (ind+1) - x) + Pkt (ind+1) * (x - xkt ( & ind) ) ) / dk ! write (*,'(5x," k=",g12.4," table=",2g12.4," P=",3g12.4)') ! & x,xkt(ind),xkt(ind+1), ! & Ppk,Pkt(ind),Pkt(ind+1) Return END FUNCTION Ppk !--------------------------------------- SUBROUTINE ReadTable (hsmall) ! interpolate table with p(k) !--------------------------------------- PARAMETER (NtabM = 100000) COMMON / PTABLE / xkt (0:NtabM), Pkt (0:NtabM), Ntab, StepK, & alog0, iFlag CHARACTER(80) Line If (iFlag.eq. - 1) Then ! W.Hu table twop = 2. * 3.14145926**2 Ntab = 0 10 READ (50, *, end = 30, err = 30) xx, pp Ntab = Ntab + 1 xkt (Ntab) = xx * hsmall ! scale to real Mpc Pkt (Ntab) = pp / xx**3 * twop ! k^3Pk -> Pk GOTO 10 30 Write (*,*) Write (*,*) ' Read ', Ntab, ' lines from WHUpk.dat table' Write (*,*) ' hsmall =', hsmall Else ! read header Do i = 1, 9 READ (50, '(a)', end = 100, err = 100) Line EndDo Ntab = 0 12 READ (50, *, end = 32, err = 32) xx, pp Ntab = Ntab + 1 xkt (Ntab) = xx * hsmall ! scale to real Mpc Pkt (Ntab) = pp ! Pk GOTO 12 32 Write (*,*) Write (*,*) ' Read ', Ntab, ' lines from P(k) table ' Write (*,*) ' hsmall =', hsmall EndIf ! end case iFlag If (Ntab.le.1) stop 'wrong table for p(k)' StepK = log10 (xkt (2) / xkt (1) ) alog0 = log10 (xkt (1) ) ! test that spacing is the same Do k = 2, Ntab - 1 ss = log10 (xkt (k + 1) / xkt (k) ) If (abs (ss / StepK - 1.) .gt.2.e-2) Then Write (*,*) ' error in K spacing. k=', k, xkt (k + 1),xkt (k) STOP EndIf EndDo Return 100 Write (*,*) ' Error reading table of P(k)' STOP END SUBROUTINE ReadTable !-------------------------------------------------------------------- ! INITIALIZE CONSTANTS: ! Scalel =length of square in MPC ! AEXPN =expansion factor (current) ! ASTEP =increment of expansion factor: AEXPN =AEXP_0+N*ASTEP ! PARTW =weight of a particle: total 'mass' per cell ! even for LCDM or OCDM ! AMPLT =amplitude of perturbations inside the simulation box ! SUBROUTINE InitValues (SCALEL) !-------------------------------------------- PARAMETER (NtabM = 100000) COMMON / PTABLE / xkt (0:NtabM), Pkt (0:NtabM), Ntab, StepK, & alog0, iFlag COMMON / TRUNCOM / Om, Omb, Omc, Omnu, Par (6), ns, qqscaleb, & QSCALE COMMON / FERMI / Vscale Real ns Character Answer*1 !Real rINPUT DATA PI / 3.1415926535 / Do i = 1, 100 extras (i) = 0. enddo Write (*,*) 'Would you like to use a model provided in cdm.fit ' Write (*,*) 'or your own model?' Write (*,'(A,$)') ' Enter Yes for cdm.fit; No for your model =' Write (*,*) Read (*,'(A)') Answer HEADER = 'N=128x256 L=20h-1CDMsig=0.239 z=30-----------' Write (*,*) '------ Enter Header for the run up to 45 characters' Write (*,*) ' Example:' Write (*,'(A)') HEADER Write (*,*) READ (*,'(A)') HEADER Write (*,'(A)') HEADER AEXPN = rINPUT (' Initial expansion parameter (0-1)=') ASTEP = rINPUT (' Step of the expansion parameter =') AMPLT = rINPUT (' Amplitude of density fluctuations=') ns = rINPUT (' Slope of the Power spectrum n=') SCALEL= rINPUT (' Box size (Mpc/h) =') Ocurv = rINPUT (' Omega_curvature at present =') Nseed = rINPUT (' Random seed (Integer 1-2^31) =') If (Nseed.le.0) Then lux = 2 ! level for luxury extras (80) = 1. extras (81) = lux If (Nseed.eq.0) Nseed = 121071 Nseed = ABS (Nseed) CALL rluxgo (lux, Nseed, 0, 0) ! initialize luxury Else ! use randd and gauss extras (80) = 0. EndIf ! extras(80) is used as switch for ! random numbers generators ! = 0 - old Randd+Gauss ! = 1 - luxury+Gauss3 ! extras(81) = lux - level of luxury ! store the box size extras (100) = SCALEL If (Answer.eq.'Y'.or.Answer.eq.'y') Then CALL MODEL (hsmall) hubble = hsmall Om0 = Om Oml0 = 1. - Om0 - Ocurv Else Write (*,*) ' You use your cosmological model.' Write (*,*) ' Be sure you provide routine TRUNF' hubble = rINPUT (' The Hubble constant (h=H/100) =') Om0 = rINPUT (' Omega_0 matter at present =') Oml0 = rINPUT (' Omega_lambda at present =') EndIf ! scale it to real megaparsecs SCALEL = SCALEL / hubble NBYTE = NPAGE * 6 * 4 Write (*, '(A,$)') ' Enter Min and Max number of mass species =' READ ( *, * ) Nmin_lev, Nmax_lev W = (FLOAT (NGRID) / FLOAT (NROW) ) **3 PARTW = W ! multiple mass resolution If (Nmin_lev.lt.Nmax_lev) Then Write ( *, * ) Write (*,*) ' You use multiple mass resolution' ENDIF ! change the order of levels If(Nmin_lev.gt.Nmax_lev)Then i = Nmax_lev Nmax_lev = Nmin_lev Nmin_lev = i endif ISTEP = 0 TINTG = 0. AU0 = 0. AEU0 = 0. EKIN = 0. EKIN1 = 0. EKIN2 = 0. NROWC = NROW NGRIDC= NGRID QSCALE = 2.*PI/SCALEL QS = hubble**(-2) write (16,'(/3x,3(a,g12.3))') ' Qscale=',QSCALE, & ' Box=',SCALEL,' slope=',ns RETURN END SUBROUTINE InitValues !-------------------------------------------------------------------- ! Broadcast Initial values ! SUBROUTINE Bcast_InitValues(SCALEL) !-------------------------------------------- use mpi PARAMETER (NtabM = 100000) COMMON / PTABLE / xkt (0:NtabM), Pkt (0:NtabM), Ntab, StepK, & alog0, iFlag COMMON / TRUNCOM / Om, Omb, Omc, Omnu, Par (6), ns, qqscaleb, & QSCALE COMMON / FERMI / Vscale Real ns if(rank < 33)write(16,*) ' Inside Bcast_Init. AEXPN=',AEXPN CALL MPI_Bcast(AEXPN,1,MPI_REAL,0,MPI_COMM_WORLD,ierr) CALL MPI_Bcast(ASTEP,1,MPI_REAL,0,MPI_COMM_WORLD,ierr) CALL MPI_Bcast(AMPLT,1,MPI_REAL,0,MPI_COMM_WORLD,ierr) CALL MPI_Bcast(SCALEL,1,MPI_REAL,0,MPI_COMM_WORLD,ierr) CALL MPI_Bcast(Nseed,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr) CALL MPI_Bcast(AMPLT,1,MPI_REAL,0,MPI_COMM_WORLD,ierr) CALL MPI_Bcast(Nmax_lev,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr) CALL MPI_Bcast(Nmin_lev,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr) CALL MPI_Bcast(QSCALE,1,MPI_REAL,0,MPI_COMM_WORLD,ierr) CALL MPI_Bcast(hubble,1,MPI_REAL,0,MPI_COMM_WORLD,ierr) CALL MPI_Bcast(Om0,1,MPI_REAL,0,MPI_COMM_WORLD,ierr) CALL MPI_Bcast(Oml0,1,MPI_REAL,0,MPI_COMM_WORLD,ierr) CALL MPI_Bcast(Ocurv,1,MPI_REAL,0,MPI_COMM_WORLD,ierr) CALL MPI_Bcast(PARTW,1,MPI_REAL,0,MPI_COMM_WORLD,ierr) CALL MPI_Bcast(EKIN,1,MPI_REAL,0,MPI_COMM_WORLD,ierr) CALL MPI_Bcast(extras,100,MPI_REAL,0,MPI_COMM_WORLD,ierr) CALL MPI_Bcast(xkt,1+NtabM,MPI_REAL,0,MPI_COMM_WORLD,ierr) CALL MPI_Bcast(Pkt,1+NtabM,MPI_REAL,0,MPI_COMM_WORLD,ierr) CALL MPI_Bcast(Ntab,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr) CALL MPI_Bcast(Ntab,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr) CALL MPI_Bcast(StepK,1,MPI_REAL,0,MPI_COMM_WORLD,ierr) CALL MPI_Bcast(alog0,1,MPI_REAL,0,MPI_COMM_WORLD,ierr) if(rank < 33)write(16,*) ' Inside Bcast_Init. rank =',rank if(rank < 33)write(16,*) ' Inside Bcast_Init. AEXPN=',AEXPN if(rank < 33)write(16,*) ' Inside Bcast_Init. ScaleL=',SCALEL End SUBROUTINE BCAST_INITVALUES !________________________________________Read parameters of the model ! Om = Omega_0 ! Omb = Omega_baryon ! Omnu = Omega_neutrino ! hsmall = hubble constant ! Par = fitting parameters ! Par(6) = 0 --- bbks-style+Hu&Sugiyama ! Par(6) ne 0 --- holtzman-style SUBROUTINE MODEL (hsmall) !--------------------------------------- COMMON / TRUNCOM / Om, Omb, Omc, Omnu, Par (6), ns, qqscaleb, & QSCALE PARAMETER (NtabM = 100000) COMMON / PTABLE / xkt (0:NtabM), Pkt (0:NtabM), Ntab, StepK, & alog0, iFlag Real ns !, rINPUT Character Header1*79 Line1 = rINPUT (' Enter Line Number in cdm.fit =') If (Line1.gt.0) Then ! use analytical approximation iFlag = Line1 OPEN (2, file = 'cdm.fit') Read (2, '(A)') Header1 If (Line1.gt.2) Then Do i = 1, Line1 - 2 Read (2, * ) a EndDo EndIf Read (2, * ) Om, Omb, Omc, Omnu, hsmall, Par CLOSE (2) ! = T_cmb/2.7 theta = 2.726 / 2.7 ! Hu&Sugiyama fitting Ob0 = Omb / Om Omh2 = Om * hsmall**2 a1 = (46.9 * Omh2) **0.670 * (1. + (32.1 * Omh2) ** ( - 0.532) ) a2 = (12.0 * Omh2) **0.424 * (1. + (45. * Omh2) ** ( - 0.582) ) alph = 1. / a1**Ob0 / a2** (Ob0**3) qqscaleb = theta**2 / (1. - Ob0) **0.60 / sqrt (alph) / Omh2 Write (*,20) Om, Omb, Omc, Omnu, hsmall, Par ! Write (16,20)Om,Omb,Omc,Omnu,hsmall,Par 20 Format (' Model: Om_0=', F5.3, ' Om_baryon=', F5.3, ' Om_cold=', & F5.3, ' O_nu=', F5.3, ' hsmall=', F4.2, / 8x, 'Parameters=', & (6G9.3) ) Else ! read spectrum from table iFlag = Line1 iFlag0 = -iFlag If(iFlag0 == 1)Then Call ReadPkTable(hsmall) Om0 = Omc Else hsmall = ParSet(1,iFlag0) Omb = ParSet(2,iFlag0) Om = ParSet(3,iFlag0) Omc = Om Open(50,file=TableNames(iFlag0)) CALL ReadTable (hsmall) Close (50) End If Write ( *, 22) Om, Omb, hsmall,TableNames(iFlag0) Write (1, 22) Om, Omb, hsmall,TableNames(iFlag0) 22 Format (' Model: Om_0=', F5.3, ' Om_baryon=', F5.3, ' hsmall=', & F4.2, ' Using table for P(k)=>',a) EndIf Return END SUBROUTINE MODEL !--------------------------------------- SUBROUTINE ReadPkTable(hsmall) ! interpolate table with p(k) !--------------------------------------- PARAMETER (NtabM = 100000) COMMON / PTABLE / xkt (0:NtabM), Pkt (0:NtabM), Ntab, StepK, & alog0, iFlag Character*80 pkTable,Line Logical :: exst Inquire(file='PkTable.dat',exist=exst) If(exst)Then Open(50,file='PkTable.dat') Else write(*,*) ' Table with the power', & ' spectrum PkTable.dat is not found' stop End If Call ReadHeader(hsmall) ! read header, get Om0, sigma8 Ntab = 0 12 READ (50, *, end = 32, err = 32) xx, pp Ntab = Ntab + 1 xkt (Ntab) = xx*hsmall Pkt (Ntab) = pp ! Pk GOTO 12 32 Write (*,*) Write (*,*) ' Read ', Ntab, ' lines from P(k) table ' If (Ntab.le.1) stop 'wrong table for p(k)' StepK = log10 (xkt (Ntab) / xkt (1) )/(Ntab-1) alog0 = log10 (xkt (1) ) ! test that spacing is the same Do k = 2, Ntab - 1 ss = log10 (xkt (k + 1) / xkt (k) ) If (abs (ss / StepK - 1.) .gt.2.e-2) Then Write (*,*) ' error in K spacing. k=', k, xkt (k + 1),xkt (k) STOP EndIf EndDo END SUBROUTINE ReadPkTable !------------------- Subroutine ReadHeader(hsmall) COMMON / TRUNCOM / Om, Omb, Omc, Omnu, Par (6), ns, qqscaleb, & QSCALE Real ns, k Omb = ParseLine() Omc = ParseLine() OmL = ParseLine() Om0 = ParseLine() sigma8 = ParseLine() hsmall =0.70 Ocurv =0. Om =Om0 ns =0.96 write(*,'(3(a,ES12.3))') & ' Model: Ombar =',Omb,' Om_matter =',Om0,' Sigma8 =',sigma8 end Subroutine ReadHeader !------------------- ! read line from Header of Pk ! Function ParseLine() COMMON / TRUNCOM / Om, Omb, Omc, Omnu, Par (6), ns, qqscaleb, & QSCALE Real ns, k Character :: Line*120,Line2*120,Line3(120) Read(50,'(a)')Line Ieq =INDEX(Line,'=',BACK=.TRUE.) !write(*,*) ' Ieq =',Ieq backspace (50) write(Line2,'(a1,i2,a)') '(',Ieq,'a1,g12.5)' !write(*,'(a)') Line2 Read(50,Line2)(Line3(i),i=1,Ieq),dummy ParseLine = dummy !write(*,'(a,ES12.3)') ' Result =',ParseLine end Function ParseLine !------------------------------------------------ ! Clear MASK for multiple mass resolution ! Initiat to level Nmin_lev SUBROUTINE ZeroMask !------------------------------------------------ DO MK2 = 1, Nblocks DO MK1 = 1, Nblocks iMask (MK1, MK2) = Nmin_lev ENDDO ENDDO RETURN END SUBROUTINE ZeroMask !------------------------------------------------ ! Set MASK for multiple mass resolution ! Read particles with lagrangian coordinat ! to set the mask ! Block = 1 = low resolution ! = 2 = higher resolution ! = 3 = even higher resolution ! After the mask is set, make maximum resolution Nmax_lev ! SUBROUTINE SetMask !------------------------------------------------ use mpi Dimension nLevel (10) Character *70 Line Integer*4, ALLOCATABLE, DIMENSION(:,:,:) :: iMs Integer*4 :: status(MPI_STATUS_SIZE) Integer*8 :: icurrent If(rank == 0)Then ALLOCATE(iMs(Nblocks,Nblocks,Nblocks)) DO MK3 = 1, Nblocks DO MK2 = 1, Nblocks DO MK1 = 1, Nblocks iMs (MK1, MK2,MK3) = Nmin_lev ENDDO ENDDO ENDDO Do i = 1, 6 ! count lines in the file Read (30, '(A)') Line Write (*, '(A)') Line EnDDo Nlines = 0 DO Read (30, *,iostat=ist) i, xw, yw, zw, vwx, vwy, vwz, im, jm, km !write (*,'(i10,6g12.4,5x,3i5)') i, xw, yw, zw, vwx, vwy, vwz, im, jm, km if(ist/=0)exit Nlines = Nlines + 1 If (im.lt.0.or.im.gt.Nblocks) Then Write (16,*) ' wrong x_block: ', im, jm, km, ' Line=', Nlines CALL MPI_FINALIZE(ierr) Stop EndIf If (jm.lt.0.or.jm.gt.Nblocks) Then Write (16,*) ' wrong y_block: ', im, jm, km, ' Line=', Nlines CALL MPI_FINALIZE(ierr) Stop EndIf If (km.lt.0.or.km.gt.Nblocks) Then Write (16,*) ' wrong z_block: ', im, jm, km, ' Line=', Nlines CALL MPI_FINALIZE(ierr) Stop EndIf ! mark high resolution block iMs (im, jm, km) = Nmax_lev Do ik = - 1, 1 k = km + ik If (k.le.0) k = k + Nblocks If (k.gt.Nblocks) k = k - Nblocks Do ij = - 1, 1 j = jm + ij If (j.le.0) j = j + Nblocks If (j.gt.Nblocks) j = j - Nblocks Do ii = - 1, 1 i = im + ii If (i.le.0) i = i + Nblocks If (i.gt.Nblocks) i = i - Nblocks ind = abs (ik) + abs (ij) + abs (ii) ! only seven blocks ! are marked out If (ind.le.1) Then iMs (i, j, k) = Nmax_lev ! change mask Endif ENDDO ENDDO ENDDO EndDo ! end reading Write (*,*) ' Total Lines read=', Nlines ! Change the mask: iteratevely add layers of lower resolution. ! Find blocks adjasent to filled blocks and ! mark them with kSpec =Nspecies-1, -2, ! DO kSpec = Nmax_lev - 1, 2, - 1 DO MK3 = 1, Nblocks DO MK2 = 1, Nblocks DO MK1 = 1, Nblocks iM = iMs (MK1, MK2, MK3) If (iM.eq.kSpec + 1) Then ! change to kSpec if not high res Do ik = - 1, 1 k = MK3 + ik If (k.le.0) k = k + Nblocks If (k.gt.Nblocks) k = k - Nblocks Do ij = - 1, 1 j = MK2 + ij If (j.le.0) j = j + Nblocks If (j.gt.Nblocks) j = j - Nblocks Do ii = - 1, 1 i = MK1 + ii If (i.le.0) i = i + Nblocks If (i.gt.Nblocks) i = i - Nblocks If (iMs (i, j, k) .lt.kSpec) iMs (i, j, k) = kSpec ENDDO ENDDO ENDDO EndIf ! end iM=kSpec+1 test ENDDO ! end MK1 ENDDO ! end MK2 ENDDO ! end MK3 Enddo ! end kSpec Do i = 1, Max_lev_abs ! count marked blocks nLevel (i) = 0 Enddo DO MK3 = 1, Nblocks DO MK2 = 1, Nblocks DO MK1 = 1, Nblocks iM = iMs (MK1, MK2, MK3) If (iM.le.Nmax_lev.and.iM.ge.Nmin_lev) Then nLevel (iM) = nLevel (iM) + 1 Else Write (*,*) 'wrong mask: ', iM, ' block=', MK1, MK2, MK3 Write (*, * ) ' it should be within limits:', Nmin_lev, Nmax_lev Stop EndIf ENDDO ENDDO ENDDO Ntot = 0 Do i = Max_lev_abs, 1, - 1 icurrent = nLevel (i) * 8** (i - 1) Write (*,*) ' Level=', i, ' Number of Blocks=', nLevel (i) , & ' N_Particles=', icurrent Ntot = Ntot + icurrent Enddo Write (*,*) ' Total expected particles =', Ntot EndIf ! rank == 0 If(rank == 0)Then Do k =2,NROW kk = (k-1)/Lblock+1 iMask = iMs(:,:,kk) CALL MPI_Send(iMask,Nblocks**2,MPI_INTEGER,k-1,k, & MPI_COMM_WORLD,ierr) EndDo iMask = iMs(:,:,1) write(*,*) ' Root have sent iMask to all nodes' DEALLOCATE(iMs) Else CALL MPI_Recv(iMask,Nblocks**2,MPI_INTEGER,0,MPI_ANY_TAG, & MPI_COMM_WORLD,status,ierr) write(16,*) ' rank=',rank,' got iMask:',iMask(1,1) endIf !Do j=1,Nblocks ! write(16,'(32i3)')(iMask(i,j),i=1,Nblocks) !EndDO RETURN END SUBROUTINE SetMask !------------------------------------------------ ! Make a realization of spectrum of perturbations ! ALPHA = normalization factor for displacement ! NRAND = seed for random numbers SUBROUTINE SPECTR (NRAND, iDirection) !------------------------------------------------ use Randoma use mpi INCLUDE 'luxuryp.h' REAL(8) :: SUMM=0., Wi3, Wj3, Wk3, WD, TS, TRX REAL :: t0=0,t1=0. ! counters for timing INTEGER*4, DIMENSION(NROW) :: mapz, map3 Real*4, DIMENSION(NROW) :: bbuff_send,bbuff_rec INTEGER*4 :: status(MPI_STATUS_SIZE) k = rank +1 ! current page CALL CPU_TIME(t0) ! initialize time counter iRandom = INT (extras (80) + 1.e-5) ! set type of Random Numbers ! 0 - old GAUSS 1-GAUSS3(ranlux) SUMM = 0. !CALL MPI_Finalize(ierr) !write(16,*) ' In Spectr. SUM=',SUMM,' rank=',rank !stop ' --------------' Do j=1,NROW jsign = - 1 jz = j + NSPEC j3 = j - 1 If (j3.GT.NSPEC) j3 = j3 - NSPEC IF (jz.gt.NROW) THEN jz = jz - NROW jsign = 1 ENDIF mapz(j) = jz map3(j) = j3*jsign EndDo GR = 0. If (iRandom.eq.0) Then NRAND = i24A(k) Else i24 = i24A(k) j24 = j24A(k) in24 = in24A(k) kount = kountA(k) mkount = 0 carry = carryA(k) iflag = iFlagA(k) gset = gSetA(k) Do m=1,24 seeds(m) =SeedsPage(m,k) EndDo EndIf ! end iRandom Wk3 = map3(k)**2 DO j = 1, NROW Wj3 = map3(j)**2 DO i = 1, NROW Wi3 = map3(i)**2 If (iRandom.eq.0) Then TS = GAUSS (NRAND) Else TS = GAUSS3 (gSet,iFlag) EndIf IF (i + j + k.EQ.3) THEN GR (1, 1) = 0. ELSE WD = Wi3 + Wj3 + Wk3 Wk = SQRT (WD) TS = TRUNF (Wk) * TS TRX = TS / WD ! start switch If (iDirection.eq.1) Then GR (mapz(i), j) = TRX * map3(i) ElseIF (iDirection.eq.2) Then GR (i, mapz(j)) = TRX * map3(j) Else !GR (i, j, mapz(k)) = TRX * map3(k) GR (i, j) = TRX * map3(k) EndIf ! end iDirection SUMM = SUMM + TS**2 ENDIF ! end kx=ky=kz=0 switch ENDDO ! i ENDDO ! j TS = SUMM ! make global summ CALL MPI_Reduce(TS,SUMM,1,MPI_DOUBLE_PRECISION, & MPI_SUM,0,MPI_COMM_WORLD,ierr) CALL MPI_Bcast(SUMM,1,MPI_DOUBLE_PRECISION, 0,& MPI_COMM_WORLD,ierr) If (iDirection.eq.3) Then ! swap pages i_send = mapz(k) -1 !write(16,*) ' Exchange pages=',i_send+1,k !CALL MPI_Sendrecv(GR, NROW**2,MPI_REAL,i_send,99, & ! GRbuff, NROW**2,MPI_REAL,i_send,99, & ! MPI_COMM_WORLD,status,ierr) !GR = GRbuff Do j=1,NROW bbuff_send = GR(:,j) CALL MPI_Sendrecv(bbuff_send, NROW,MPI_REAL,i_send,99, & bbuff_rec, NROW,MPI_REAL,i_send,99, & MPI_COMM_WORLD,status,ierr) GR(:,j) = bbuff_rec EndDo !Do j=1,min(32,NROW) ! write(16,'(1P,16g12.4)') (GR(i,j),i=1,16) !EndDo EndIf IF (SUMM.LE.0.) Write (*, * ) ' Error!!! Summ over spectrum = 0' ALPHA = AMPLT / SQRT (SUMM) * sqrt (8.) if(rank < 33)Write (16,'(10x,a20,3g13.6)') 'SPECTR ==> SUMM=', SUMM, ALPHA CALL CPU_TIME(t1) timing(3) = timing(3) +t1-t0 ! time for setting SPECTR !write(16,*) ' ---------- k=',k ! Do j=1,32 ! write(16,'(32g12.3)') (GR(i,j),i=1,16) ! EndDo RETURN END SUBROUTINE SPECTR !------------------------------------------------ ! Define coordinates and velocities for given plane (i,j) ! particles with resolution = iMs(i,j) ! icurrent = number of particles SUBROUTINE BLOCKS (XCONS, VCONS, iDirection) !------------------------------------------------ use mpi Real(8) :: sDispl =0. REAL :: t0=0,t1=0. ! counters for timing Integer*8 :: icurrent Character*80 :: FileName CALL CPU_TIME(t0) ! initialize time counter nLevel =0 Boff = Lblock / 2 + 0.5 sDispl = 0. xShift = 0.5 /2**LevelParticles krank = rank +1 Mk3 = (krank-1)/Lblock+1 M3 = (Mk3 -1)*Lblock icurrent = 0 iFlag = 0 XMAX = FLOAT (NGRID) + 1. XSHF = FLOAT (NGRID) MPt(:) = 0 ; XPt(:) = 0. ; VXt(:) = 0. CALL MPI_COMM_SIZE(MPI_COMM_WORLD, numprocs,ierr) lchunk =min(numprocs,256) ! length of a chunk of MPI tasks for disk IO nchunks = (numprocs-1)/lchunk+1 mychunk = (krank-1)/lchunk+1 Do m=1,nchunks ! create output file for every thread CALL MPI_Barrier(MPI_COMM_WORLD,ierr) If(mychunk == m)Then write(FileName,'(a,i1.1,a,i4.4,a)') 'PMc/PMcoord.',iDirection,'.',krank,'.DAT' QFACT = FLOAT (NGRID) / FLOAT (NROW) !write(16,*) ' Blocks: Mk3=',Mk3,M3,icurrent !write(16,*) ' Blocks: xShift=',xShift,' iMask=',iMask(1,1),iMask(2,2) !write(16,*) ' Blocks: Nblocks=',Nblocks !CALL MPI_Finalize(ierr) !stop ' BL--------------' DO Mk2 = 1, Nblocks M2 = (Mk2 - 1) * Lblock ! part of offset DO Mk1 = 1, Nblocks M1 = (Mk1 - 1) * Lblock ! part of offset iM = iMask (Mk1, Mk2) ! KBlock = 2**(iM-1) = number of particles in 1 ! Nbcells= Lblock/KBlock = number of cells per particl ! (MK-1)*Lblock +(i-1)*Nbcells = offset of cells in 1D ! (MK-1)*Lblock +(i-1)*Nbcells +Nbcells/2 +1/2 = lagra ! coordi ! M1 = (Mk1-1)*Lbloc KBlock = 2** (iM - 1) ! number of particles in 1D per block Nbcells = Lblock / KBlock ! number of cells per particle in 1D Xlag = Nbcells / 2 + 0.5 ! offset for lagrangian coordinates Nbcells3 = Nbcells**3 ! total number of cells per particle Do k = 1, KBlock k1 = M3 + (k - 1) * Nbcells ! offset for cells in Z Q3 = QFACT * (k1 + Xlag - 1.) + 1. kp = k1 + (Nbcells + 1) / 2 !If(icurrent == 100)Then ! write(16,*)icurrent ,k,kp ! CALL MPI_Finalize(ierr) ! stop ' BL 2--------------' !EndIf If(kp == krank)Then ! If this is my plane, get particles Do i = 1, KBlock ! loop over particles i1 = M1 + (i - 1) * Nbcells ! offset for cells in X ip = i1 + (Nbcells + 1) / 2 Q1 = QFACT * (i1 + Xlag - 1.) + 1. ! scale it to box size Do j = 1, KBlock j1 = M2 + (j - 1) * Nbcells ! offset for cells in Y jp = j1 + (Nbcells + 1) / 2 Q2 = QFACT * (j1 + Xlag - 1.) + 1. DX = GR (ip, jp) icurrent = icurrent + 1 If(iFlag == 0)Then Open(30,file=TRIM(FileName),form='unformatted') iFlag = 1 !write(16,*)icurrent ,Q1,Q2,Q3 EndIf nLevel (iM) = nLevel (iM) + 1 If (iDirection.eq.1) Then Q = Q1 ElseIF (iDirection.eq.2) Then Q = Q2 Else Q = Q3 EndIf XPt (icurrent) = Q - XCONS * DX + xShift IF (XPt(icurrent) .GT.XMAX) XPt(icurrent) = XPt(icurrent) - XSHF IF (XPt(icurrent) .LE.1.) XPt(icurrent) = XPt(icurrent) + XSHF If (XPt(icurrent) .ge.XMAX)XPt(icurrent) = XPt(icurrent) - 1.e-4 !write(16,'(3i8,8g13.5)')ip,jp,icurrent,XPt (icurrent),Q,DX,xShift,XCONS VXt (icurrent) = VCONS * DX MPt (icurrent) = iM sDispl = sDispl + DX**2 If (icurrent.eq.NPAGE) Then if(rank < 33)Write (16,*) ' Write a page of particles: N=', icurrent Do kk =1,icurrent If(XPt(kk)<1.)& write(16,*) ' Error coordinates:',kk,XPt(kk),VXt(kk),MPt(kk) EndDo write(30)XPt,VXt,MPt MPt(:) = 0 ; XPt(:) = 0. ; VXt(:) = 0. icurrent = 0 EndIf EndDo ! end j EndDo ! end j EndIf ! kp ==k EndDo ! end k EndDo ! end Mk1 =1,Nblocks EndDo ! end Mk2 =1,Nblocks If(iFlag==1)Then ! there are particles on this task If(icurrent/=0)Then ! flash the buffer Write (16,*) ' Write a page of particles:', icurrent Do kk =1,icurrent If(XPt(kk)<1.)& write(16,*) ' Error coordinates:',kk,XPt(kk),VXt(kk),MPt(kk) EndDo write(30)XPt,VXt,MPt icurrent = 0 EndIf close (30) icurrent = SUM(nLevel) if(rank < 33)Write (16, '(10x,a20,2g13.5,i10)') ' RMS =',XCONS,sDispl if(rank < 33)Write (16, '(10x,a,g12.4)') ' RMS 3d diplacement=', & XCONS*sqrt (sDispl/icurrent ) if(rank < 33)write (16,'(10x,a,i11,a,i10,a,6i10)') & ' Current Particles=',icurrent,' on all levels:', & (nLevel(i),i=1,Max_lev_abs ) EndIf EndIf ! mychunk EndDo ! m CALL MPI_Gather(nLevel,10,MPI_LONG,Allspecies,10,MPI_LONG,0,MPI_COMM_WORLD,ierr) ! root gets all information on nLevel on all tasks If(rank ==0)then nLevel = 0 Do i =1,10 Do j =1,NROW nLevel(i) = nLevel(i) + Allspecies(i,j) ! final table of particles EndDo EndDo Nspecies = 0 Do i =10,1,-1 ! find Nspecies If(nLevel(i) /= 0)Then Nspecies = i; exit EndIf EndDo write (16,'(10x,a20,i4,a,10i10)') ' Species=',Nspecies, & ' on all levels:', & (nLevel(i),i=1,Nspecies) If(iDirection == 3)Then open(70,file='Allspecies.DAT') write (70,'(i4,a)') Nspecies,' <= Nspecies. On all different levels=>' write(70,'(10i12)')(nLevel(i),i=1,Nspecies) Do j=1,NROW write(70,'(10i12)')(Allspecies(i,j),i=1,Nspecies) EndDo close (70) EndIf EndIf CALL CPU_TIME(t1) timing(4) = timing(4) +t1-t0 ! time for setting SPECTR RETURN END SUBROUTINE BLOCKS !------------------------------------------------ ! FFT of the spectrum SUBROUTINE VECTOR (IFOUR) !------------------------------------------------ use mpi DIMENSION Uf (marr), Vf (marr) DIMENSION Qi (mf67), jndx (mf67) REAL :: t0=0,t1=0. ! counters for timing Real*4, DIMENSION(NROW) :: bbuff_send,bbuff_rec CALL CPU_TIME(t0) jb1 = 4 jq = IFOUR CALL SETF67 (jb1, jq, jbc, jp, jsl, ll1, k2, k3, k4, k7, Qi, jndx,& Uf, Vf) k = rank +1 ! !.... Open_MP -1- !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE ( k,j,i,Uf,Vf,jp,jsl,ll1,k2,k4 ) DO J = 1, NROW DO I = 1, NROW Uf (I) = GR (I, J) ENDDO CALL FOUR67 (jb1, jq, jp, jsl, ll1, k2, k7, Qi, jndx, Uf, Vf) DO I = 1, NROW GR (I, J) = Vf (I) ENDDO ENDDO !Write (16,'(10x,a)') 'Xfft is done' !.... Open_MP -2- !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE ( k,j,i,Uf,Vf,jp,jsl,ll1,k2,k4 ) DO I = 1, NROW DO J = 1, NROW Uf (J) = GR (I, J) ENDDO CALL FOUR67 (jb1, jq, jp, jsl, ll1, k2, k7, Qi, jndx, Uf, Vf) DO J = 1, NROW GR (I, J) = Vf (J) ENDDO ENDDO !Write (16,'(10x,a)') 'Yfft is done' !write (16,'(10x,a)') 'Swap k<->j' !CALL MPI_ALLtoALL(GR,NROW,MPI_INTEGER, & ! GRbuff,NROW,MPI_INTEGER,MPI_COMM_WORLD,ierr) !GR = GRbuff !write(16,*) Do i=1,NROW bbuff_send = GR(i,:) CALL MPI_ALLtoALL(bbuff_send,1,MPI_INTEGER, & bbuff_rec,1,MPI_INTEGER,MPI_COMM_WORLD,ierr) GR(i,:) = bbuff_rec EndDo !.... Open_MP -3- !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE ( k,j,i,Uf,Vf,jp,jsl,ll1,k2,k4 ) DO I = 1, NROW DO J = 1, NROW Uf (J) = GR (I, J) ENDDO CALL FOUR67 (jb1, jq, jp, jsl, ll1, k2, k7, Qi, jndx, Uf, Vf) DO J = 1, NROW GR (I, J) = Vf (J) / 8. ENDDO ENDDO !write (16,'(10x,a)') 'Swap j<->k' !CALL MPI_ALLtoALL(GR,NROW,MPI_INTEGER, & ! GRbuff,NROW,MPI_INTEGER,MPI_COMM_WORLD,ierr) !GR = GRbuff Do i=1,NROW bbuff_send = GR(i,:) CALL MPI_ALLtoALL(bbuff_send,1,MPI_INTEGER, & bbuff_rec,1,MPI_INTEGER,MPI_COMM_WORLD,ierr) GR(i,:) = bbuff_rec EndDo !write(16,*) ' --- page ---' !Do j=1,min(32,NROW) ! write(16,'(1P,16g12.4)') (GR(i,j),i=1,16) !EndDo if(rank < 33)Write (16,'(10x,a)') ' FFT is done' CALL CPU_TIME(t1) timing(5) = timing(5) +t1-t0 ! time for VECTOR RETURN END SUBROUTINE VECTOR !--------------------------------------------- ! Read files of individual pages ! Write final data to disk ! SUBROUTINE FilesWrite !---------------------------------------------- ! write header and control data use mpi REAL :: t0=0,t1=0. ! counters for timing Integer*8 :: j0,jfirst,jlast,NinPage,JPAGE Character*80 :: fname integer*4 :: mrank, ranks(0:Nslaves) Integer*4 :: map(NROW),mapz(NROW) Integer*4 :: status(MPI_STATUS_SIZE),ogroup,gr_slave,com_slave Integer*4 :: request(Nslaves),status_all(MPI_STATUS_SIZE,Nslaves),mes(3) Integer*4:: irequest(7),NpS(Nslaves),brequest(Nslaves),status_buf(MPI_STATUS_SIZE,7) CALL CPU_TIME(t0) ! initialize time counter Do i=0,Nslaves ! create new communicator group for parallel IO ranks(i) =i endDo CALL MPI_COMM_GROUP(MPI_COMM_WORLD,ogroup,ierr) write(16,*) ' ogroup= ',ogroup,Nslaves+1,ranks CALL MPI_GROUP_INCL(ogroup,Nslaves+1,ranks,gr_slave,ierr) write(16,*) ' gr_slave= ',gr_slave CALL MPI_COMM_CREATE(MPI_COMM_WORLD,gr_slave,com_slave,ierr) CALL MPI_GROUP_RANK(gr_slave,mrank,ierr) If(mrank /= MPI_UNDEFINED)Then CALL MPI_Bcast(Nspecies,1,MPI_INTEGER,0,com_slave,ierr) ALLOCATE (XPAR(NPAGE),YPAR(NPAGE),ZPAR(NPAGE)) ALLOCATE(VX(NPAGE),VY(NPAGE),VZ(NPAGE)) write(*,'(2(a,i4))') ' my rank=',mrank,' Nspecies =',Nspecies ! open files and write control information------------------------- If(mrank ==0) Then ! ------------------- root --------------- OPEN (UNIT = 9, FILE = 'PMcrd.DAT', form = 'unformatted') ! this clears old header file Write (9) CLOSE (9) OPEN (9, FILE = 'PMcrd.DAT', FORM = 'UNFORMATTED', STATUS = & 'UNKNOWN') OPEN (60, file = 'pt.dat', form = 'unformatted') Write (60) astep/2**LevelParticles ! write pt.dat file: time-step for particles close (60) js = Nspecies ! find first non-zero species Do i=1,Nspecies If(nLevel(i)/=0)Then js = i ; exit EndIf EndDo j0 = 0 Do i=Nspecies,js,-1 ! count particles and reverse the order of species j0 = j0 + nLevel(i) lspecies(Nspecies-i+1) = j0 wspecies (Nspecies-i+1) = PARTW * Lblock**3 / 8** (i - 1) EndDo extras(40) = NPAGE extras(41) = Npages Write (9) HEADER, AEXPN, AEXP0, AMPLT, ASTEP, ISTEP, PARTW, TINTG,& EKIN, EKIN1, EKIN2, AU0, AEU0, NROWC, NGRIDC, Nspecies-js+1, Nseed, & Om0, Oml0, hubble, Wp5, Ocurv, extras close (9) NBYTE = NRECL*4 NACCES = NBYTE / nbyteword icurrent = 0 ! number of particles in current page lpage = 0 ! current page number JPAGE = (lspecies(Nspecies-js+1)-1)/NPAGE+1 ! number of pages Nfiles = (JPAGE-1)/Npages +1 ! number of files write(*,'(5x,a,i5)') ' Nfiles to open =',Nfiles write(*,'(5x,a,4i10)') ' Jpage,NPAGEs, Npage =',JPAGE,Npages,NPAGE write(*,'(5x,a,i5)') ' Nspecies =',Nspecies-js+1 write(*,'(5x,a,i11)') ' Ntotal particles = ',lspecies(Nspecies-js+1) write(*,'(7x,a)') 'Specie N(>) dN' Do j=1,Nspecies write(*,'(10x,i3,2i11)') j,lspecies(j),nlevel(Nspecies-j+1) EndDo Do i =1,Nfiles write(fname,'(a,i2.2,a)')'PMcrs',i-1,'.DAT' OPEN(UNIT=20+i,FILE=TRIM(fname),ACCESS='DIRECT', & FORM='unformatted',STATUS='UNKNOWN',RECL=NACCES) EndDo ALLOCATE(Buff(NPAGE,6,Nslaves)) Else ! ----------------- slaves write(*,'(5x,a,i5)') ' Slave =',INT(mrank) write(fname,'(a,i3.3,a)')'Slave.',mrank,'.DAT' open(136,file=TRIM(fname)) write(136,'(5x,a,i5)') ' Slave =',INT(mrank) ALLOCATE(YPt(NPAGE),ZPt(NPAGE),VYt(NPAGE),VZt(NPAGE)) EndIf ! CALL MPI_Finalize(ierr) ! CALL MPI_Abort(MPI_COMM_WORLD,iierr,ierr) ! stop ' --------------' CALL MPI_Bcast(js,1,MPI_INTEGER,0,com_slave,ierr) write(136,'(5x,a,i11)') ' js = ',js Do jSpecies=js,Nspecies ! write one species at a time If(mrank ==0) Then ! ------------------- root --------------- nread = 0; map =0; mapz =0 Do i=1,NROW If(Allspecies(Nspecies-jSpecies+js,i)/=0)Then ! there are particles nread = nread+1 ! counter for pages with particles map(i) = Allspecies(Nspecies-jSpecies+js,i) ! number of particles in this page mapz(nread) = i EndIf EndDo nChunks = (nRead-1)/Nslaves +1 write(*,'(5x,a,T45,i5)') ' Number of pages to read from disk=',nRead write(*,'(5x,a,T45,i5)') ' Number of slaves to handle IO =',Nslaves write(*,'(5x,a,T45,i5)') ' Number of chunks =',nChunks write(*,*)' Number of particles on different pages for current species=' write(*,'(20i8)') map write(*,*)' List of pages, which have some particles for writing to files=' write(*,'(20i6)') (mapz(i),i=1,nRead) write(*,*)' Species number = ',jSpecies CALL MPI_Bcast(nRead,1,MPI_INTEGER,0,com_slave,ierr) CALL MPI_Bcast(nChunks,1,MPI_INTEGER,0,com_slave,ierr) j0 = 0 Do ichunk= 1,nChunks ! Deal with only Nslaves pages at a time imax = 0 ! max number of particles in this chunk Do j=1,Nslaves j1 = j+(ichunk-1)*Nslaves ! global page pointer If(j1.le.nRead)Then i1 = mapz(j1) ! page number in NROW set of pages if(map(i1)>imax)imax =map(i1) EndIf EndDo Nsections = (imax-1)/NPAGE+1 ! receive this number of messages from each slave mes(3) = Nsections Do j=1,Nslaves ! send message to slaves j0 =j0+1 ! global counter for pages If(j0.le.nRead)Then mes(1) = mapz(j0) ! page to read k = j+(ichunk-1)*Nslaves mes(2) = map(mapz(k)) ! number of particles in this page at this level Else mes(1) = 0 mes(2) = 0 k = 0 EndIf write(*,'(5x,2(a,i4),a,i9,9i8)') ' line=',j0,' send to slave=',j,' NpartMax= ',imax,mes CALL MPI_Send(mes,3,MPI_INTEGER,j,j,com_slave,ierr) EndDo ! now receive all information Do k=1,Nsections NpS = 0 Do j=1,Nslaves ! get number of particles in each section CALL MPI_Recv(NpS(j),1,MPI_INTEGER,j,7,com_slave,status,ierr) EndDo ! j write(*,'(5x,a,64i8)') ' Number of particles received from slaves = ',NpS ! get all coords and Vels for a section Do j=1,Nslaves Do kk =1,6 ! loop for x,y,z,vx,vy,vz CALL MPI_Recv(Buff(1,kk,j),NPAGE,MPI_REAL,j,kk,com_slave,status,ierr) EndDo ! kk EndDo ! j Do j=1,Nslaves If(NpS(j)/=0)Then Do i=1,NpS(j) If(Buff(i,1,j)<1..or.Buff(i,2,j)<1..or.Buff(i,3,j)<1.) & write(*,'(a,2i5,3g13.5)') ' Error in root:',i,j,(Buff(i,kk,j),kk=1,3) icurrent = icurrent +1 XPAR(icurrent) = Buff(i,1,j) YPAR(icurrent) = Buff(i,2,j) ZPAR(icurrent) = Buff(i,3,j) VX(icurrent) = Buff(i,4,j) VY(icurrent) = Buff(i,5,j) VZ(icurrent) = Buff(i,6,j) If(icurrent == NPAGE)Then ! write to file lpage = lpage +1 ! current global page jfile = (lpage-1)/Npages +1 ! current file number ipage = lpage -(jfile-1)*Npages ! page in current file !If(jfile>6)& write(20+jfile,REC=ipage) XPAR,YPAR,ZPAR,VX,VY,VZ write(*,*) ' write to file=',20+jfile,' global=',lpage,' Local=',ipage write(*,'(1p,3g13.4,3x,0p,i9)') (XPAR(mm),YPAR(mm),ZPAR(mm),mm,mm=1,16) icurrent = 0 EndIf EndDo ! i EndIf ! Nps/=0 EndDo ! j =Nslaves EndDo ! k=1,Nsections EndDo ! chunks write(16,*) ' Npages=',Npages Else !-------------------------------- SLAVE ------------- CALL FilesSlave(nChunks,nRead,Nspecies-jSpecies+js,com_slave,mrank) EndIf ! mrank EndDo ! jSpecies=1,Nspecies If(mrank ==0)Then If(icurrent /= 0)Then ! write to file last set of particles lpage = lpage +1 ! current global page jfile = (lpage-1)/Npages +1 ! current file number ipage = lpage -(jfile-1)*Npages ! page in current file write(20+jfile,REC=ipage) XPAR,YPAR,ZPAR,VX,VY,VZ write(*,*) ' write to file=',20+jfile,' global=',lpage,' Local=',ipage !write(*,'(1p,3g13.4,3x,0p,i9)') (XPAR(mm),YPAR(mm),ZPAR(mm),mm,mm=1,NPAGE) icurrent = 0 EndIf end If CALL CPU_TIME(t1) timing(7) = timing(7) +t1-t0 ! EndIf ! mrank /= undefined !CALL MPI_Finalize(ierr) !CALL MPI_Abort(MPI_COMM_WORLD,iierr,ierr) ! stop ' ====-----------' CALL MPI_Barrier(MPI_COMM_WORLD,ierr) RETURN END SUBROUTINE FilesWrite !--------------------------------------------- ! Slave's part of FilesWrite: ! read and parse data, ! then send it to root SUBROUTINE FilesSlave(nChunks,nRead,jSpecies,com_slave,mrank) !---------------------------------------------- ! write header and control data use mpi REAL :: t0=0,t1=0. ! counters for timing Integer*8 :: j0,jfirst,jlast,NinPage,JPAGE Character*80 :: fname integer*4 :: mrank Integer*4 :: status(MPI_STATUS_SIZE),ogroup,gr_slave,com_slave Integer*4 :: request(Nslaves),mes(3) Integer*4 :: irequest(7),status_buf(MPI_STATUS_SIZE,7) Integer*4 :: nin(0:10) XMAX = FLOAT (NGRID) + 1. XSHF = FLOAT (NGRID) ibuff = 0 CALL MPI_Bcast(nRead,1,MPI_INTEGER,0,com_slave,ierr) CALL MPI_Bcast(nChunks,1,MPI_INTEGER,0,com_slave,ierr) write(136,'(/a,4i9)') ' Slave: Rank=',mrank,nRead,Nchunks,jSpecies !write(*,*) ' Rank=',mrank,nRead,Nchunks!,Nspecies ! Return !!! Do i =1,nChunks CALL MPI_Recv(mes,3,MPI_INTEGER,0,MPI_ANY_TAG,com_slave,status,ierr) !CALL MPI_iRecv(mes,3,MPI_INTEGER,0,MPI_ANY_TAG,com_slave,request(mrank),ierr) !CALL MPI_Wait(request(mrank),status,ierr) write(136,*) ' Slave: message= ',mes,' rank=',mrank k = mes(1) ! page to read nn = mes(2) ! number of particles to read Nsections = mes(3) ! number of sections to send back write(fname,'(a,i1.1,a,i4.4,a)') 'PMc/PMcoord.',1,'.',k,'.DAT' Open(31,file=TRIM(fName),form='unformatted') write(fname,'(a,i1.1,a,i4.4,a)') 'PMc/PMcoord.',2,'.',k,'.DAT' Open(32,file=TRIM(fName),form='unformatted') write(fname,'(a,i1.1,a,i4.4,a)') 'PMc/PMcoord.',3,'.',k,'.DAT' Open(33,file=TRIM(fName),form='unformatted') write(136,'(5x,2a)') 'Open file = ',TRIM(fName) iSec = 0 ibuff =0 XPAR = 0.; YPAR =0.; ZPAR =0. n = 0 icount = 0 If(nn>0)Then DO ! loop over records read(31,iostat=io) XPt,VXt,MPt read(32,iostat=io) YPt,VYt,MPt read(33,iostat=io) ZPt,VZt,MPt If(io /= 0)exit n=n+1 nin = 0 ! count number of particles in different species !!! test do j =1,NPAGE kk =max(min(Mpt(j),20),0) nin(kk) = nin(kk) +1 enddo write(136,'(5x,a,20i8)') 'Number of particles on different Species = ', & (nin(kk),kk=Nspecies,0,-1) !write(136,'(5x,a,30i8)') ' Nspecies-jSpecies+1 =', Nspecies-jSpecies+1, & ! (MPt(j),j=1,16) LL = ibuff Do j =1,NPAGE If(MPt(j) == jSpecies.and.icountXMAX) XPAR(ibuff) = XPAR(ibuff) - XSHF IF(YPAR(ibuff)>XMAX) YPAR(ibuff) = YPAR(ibuff) - XSHF IF(ZPAR(ibuff)>XMAX) ZPAR(ibuff) = ZPAR(ibuff) - XSHF If(XPAR(ibuff)<1..or.YPAR(ibuff)<1..or.ZPAR(ibuff)<1.)& write(136,*) ' Error coordinates:',ibuff,XPAR(ibuff),YPAR(ibuff),ZPAR(ibuff) !write(136,'(15x,1p,3g12.4,0p,i4)')XPt(j),YPt(j),ZPt(j),MPt(j) If(ibuff == NPAGE)Then CALL MPI_Send(ibuff,1,MPI_INTEGER,0,7,com_slave,ierr) CALL MPI_Send(XPAR,NPAGE,MPI_REAL,0,1,com_slave,ierr) CALL MPI_Send(YPAR,NPAGE,MPI_REAL,0,2,com_slave,ierr) CALL MPI_Send(ZPAR,NPAGE,MPI_REAL,0,3,com_slave,ierr) CALL MPI_Send(VX,NPAGE,MPI_REAL,0,4,com_slave,ierr) CALL MPI_Send(VY,NPAGE,MPI_REAL,0,5,com_slave,ierr) CALL MPI_Send(VZ,NPAGE,MPI_REAL,0,6,com_slave,ierr) write(136,*) ' sent to root N particles= ',ibuff ibuff = 0 iSec = iSec + 1 XPAR = 0.; YPAR =0.; ZPAR =0. EndIf ! ibuff EndIf ! MPt EndDo ! j EndDo ! end reading files write(136,'(5x,4(a,i7))')'read page =',n,' Ncurrent=',ibuff,' Need to get=',nn,' rank=',mrank If(ibuff /= 0)Then Do ib=1,ibuff If(XPAR(ib)<1..or.YPAR(ib)<1..or.ZPAR(ib)<1.)& write(136,'(a,i9,1p,3g13.4)') ' Error in coordinates:', & ib,XPAR(ib),YPAR(ib),ZPAR(ib) End Do CALL MPI_Send(ibuff,1,MPI_INTEGER,0,7,com_slave,ierr) CALL MPI_Send(XPAR,NPAGE,MPI_REAL,0,1,com_slave,ierr) CALL MPI_Send(YPAR,NPAGE,MPI_REAL,0,2,com_slave,ierr) CALL MPI_Send(ZPAR,NPAGE,MPI_REAL,0,3,com_slave,ierr) CALL MPI_Send(VX,NPAGE,MPI_REAL,0,4,com_slave,ierr) CALL MPI_Send(VY,NPAGE,MPI_REAL,0,5,com_slave,ierr) CALL MPI_Send(VZ,NPAGE,MPI_REAL,0,6,com_slave,ierr) ibuff = 0 iSec = iSec + 1 EndIf EndIf ! nn>0 write(136,*) ' iSec = ',iSec,' Nsections = ',Nsections If(iSec .lt.Nsections)Then ! make fake send XPAR = 0.; YPAR =0.; ZPAR =0. ; Vx =0.; VY =0.; VZ =0. Do j= iSec+1,Nsections CALL MPI_Send(0,1,MPI_INTEGER,0,7,com_slave,ierr) CALL MPI_Send(XPAR,NPAGE,MPI_REAL,0,1,com_slave,ierr) CALL MPI_Send(YPAR,NPAGE,MPI_REAL,0,2,com_slave,ierr) CALL MPI_Send(ZPAR,NPAGE,MPI_REAL,0,3,com_slave,ierr) CALL MPI_Send(VX,NPAGE,MPI_REAL,0,4,com_slave,ierr) CALL MPI_Send(VY,NPAGE,MPI_REAL,0,5,com_slave,ierr) CALL MPI_Send(VZ,NPAGE,MPI_REAL,0,6,com_slave,ierr) EndDo write(136,*) ' sent to root empty particles= ',iSec+1,Nsections EndIf ! iSec 100 continue EndDo ! chunks write(16,*) ' Npages=',Npages write(136,'(5x,4(a,i7))')' Need to get=',nn,' Sent total=',icount,' rank=',mrank End SUBROUTINE FILESSLAVE !-------------------------------------- ! normal random numbers FUNCTION GAUSS3 (gSet,iFlag) ! Uses ranlux for homogenous rand numbers ! Uses Box-Muller inversion for gaussian ! Excellent quality and speed ! N= 100000000 GAUSS3+ranlux ------ ! sigma= 0.99997 mean= 0.18692E-03 ! sigma Frac(>sig) Frac(<-sig) True n(>) n(<) ! 1.00 0.1587 0.1586 0.1587 15869095 15857724 ! 1.50 0.6682E-01 0.6680E-01 0.6681E-01 6681652 6679768 ! 2.00 0.2277E-01 0.2273E-01 0.2275E-01 2276651 2273197 ! 2.50 0.6222E-02 0.6206E-02 0.6210E-02 622190 620610 ! 3.00 0.1356E-02 0.1347E-02 0.1350E-02 135628 134674 ! 4.00 0.3242E-04 0.3145E-04 0.3170E-04 3242 3145 ! !-------------------------------------- DIMENSION RanNum (2) If (iFlag.eq.0) Then 1 CALL ranlux (RanNum, 2 ) x1 = 2. * RanNum (1) - 1. x2 = 2. * RanNum (2) - 1. R = x1**2 + x2**2 If (R.ge.1.) GoTo 1 F = sqrt ( - 2. * LOG (R) / R) gSet = x1 * F GAUSS3 = x2 * F iFlag = 1 Else GAUSS3 = gSet iFlag = 0 EndIf RETURN END FUNCTION GAUSS3 !-------------------------------------- ! FUNCTION GAUSS3fake (gSet,iFlag) ! !-------------------------------------- DIMENSION RanNum (2) If (iFlag.eq.0) Then 1 CALL ranlux (RanNum, 2 ) x1 = 2. * RanNum (1) - 1. x2 = 2. * RanNum (2) - 1. R = x1**2 + x2**2 If (R.ge.1.) GoTo 1 iFlag = 1 Else iFlag = 0 EndIf GAUSS3fake =0. RETURN END FUNCTION GAUSS3fake !-------------------------------------------------------------------- ! LUXURY LEVELS. ! level 0 (p=24): equivalent to the original RCARRY of Marsaglia ! and Zaman, very long period, but fails many tests. ! level 1 (p=48): considerable improvement in quality over level 0, ! now passes the gap test, but still fails spectral test. ! level 2 (p=97): passes all known tests, but theoretically still ! defective. ! level 3 (p=223): DEFAULT VALUE. Any theoretically possible ! correlations have very small chance of being observed. ! level 4 (p=389): highest possible luxury, all 24 bits chaotic. !--------------------------------------------- SUBROUTINE ranlux (rvec, lenv) ! returns a vector RVEC of LEN ! 32-bit random floating point numbers between ! zero (not included) and one (also not incl.). ! Next call to ranlux gives next LEN of random numbers !--------------------------------------------- INTEGER lenv REAL rvec (lenv) INTEGER iseeds (24) INCLUDE 'luxuryp.h' ! DATA ndskip / 0, 24, 73, 199, 365 / ! DATA i24, j24, luxlev / 24, 10, lxdflt / ! DATA notyet / .true. / ! DATA in24, kount, mkount, carry / 0, 0, 0, 0. / ! NOTYET is .TRUE. if no initialization has been performed yet. ! Default Initialization by Multiplicative Congruential IF (notyet) THEN stop ' Access to runlux before initialization STOP' ENDIF ! The Generator proper: "Subtract-with-borrow", ! as proposed by Marsaglia and Zaman, ! Florida State University, March, 1989 DO ivec = 1, lenv uni = seeds (j24) - seeds (i24) - carry IF (uni.LT.0.) THEN uni = uni + 1.0 carry = twom24 ELSE carry = 0. ENDIF seeds (i24) = uni i24 = next (i24) j24 = next (j24) rvec (ivec) = uni ! small numbers (with less than 12 "significant" bits) are "padde IF (uni.LT.twom12) THEN rvec (ivec) = rvec (ivec) + twom24 * seeds (j24) ! and zero is forbidden in case someone takes a logarithm IF (rvec (ivec) .EQ.0.) rvec (ivec) = twom24 * twom24 ENDIF ! Skipping to luxury. As proposed by Martin Luscher. in24 = in24 + 1 IF (in24.EQ.24) THEN in24 = 0 kount = kount + nskip DO isk = 1, nskip uni = seeds (j24) - seeds (i24) - carry IF (uni.LT.0.) THEN uni = uni + 1.0 carry = twom24 ELSE carry = 0. ENDIF seeds (i24) = uni i24 = next (i24) j24 = next (j24) ENDDO ENDIF ENDDO kount = kount + lenv IF (kount.GE.igiga) THEN mkount = mkount + 1 kount = kount - igiga ENDIF RETURN ! SUBROUTINE ranlux END SUBROUTINE ranlux !--------------------------------------------- ! Subroutine to initialize from one or three integers SUBROUTINE rluxgo (lux, ins, k1, k2) ! initializes the generator from ! one 32-bit integer INT and sets Luxury Level LUX ! which is integer between zero and MAXLEV, or if ! LUX .GT. 24, it sets p=LUX directly. K1 and K2 ! should be set to zero unless restarting at a break ! point given by output of RLUXAT !--------------------------------------------- INCLUDE 'luxuryp.h' INTEGER iseeds(24) IF (lux.LT.0) THEN luxlev = lxdflt ELSEIF (lux.LE.maxlev) THEN luxlev = lux ELSEIF (lux.LT.24.OR.lux.GT.2000) THEN luxlev = maxlev Write (6, '(A,I7)') ' RANLUX ILLEGAL LUXURY RLUXGO: ', lux ELSE luxlev = lux DO ilx = 0, maxlev IF (lux.EQ.ndskip (ilx) + 24) luxlev = ilx ENDDO ENDIF IF (luxlev.LE.maxlev) THEN nskip = ndskip (luxlev) !Write (6, '(A,I2,A,I4)') ' RANLUX LUXURY LEVEL SET BY RLUXGO :', & ! luxlev, ' P=', nskip + 24 ELSE nskip = luxlev - 24 !Write (6, '(A,I5)') ' RANLUX P-VALUE SET BY RLUXGO TO:', & !luxlev ENDIF in24 = 0 IF (ins.LT.0) Write (6,'(A)')' Illegal initialization by RLUXGO ',& 'negative input seed' IF (ins.GT.0) THEN jseed = ins !Write (6, '(A,3I12)') ' RANLUX INITIALIZED BY RLUXGO FROM SEEDS', & ! jseed, k1, k2 ELSE jseed = jsdflt !Write (6, '(A)') ' RANLUX INITIALIZED BY RLUXGO FROM DEFAULT SEED' ENDIF inseed = jseed notyet = .false. twom24 = 1. DO i = 1, 24 twom24 = twom24 * 0.5 k = jseed / 53668 jseed = 40014 * (jseed-k * 53668) - k * 12211 IF (jseed.LT.0) jseed = jseed+icons iseeds (i) = MOD (jseed, itwo24) ENDDO twom12 = twom24 * 4096. DO i = 1, 24 seeds (i) = REAL (iseeds (i) ) * twom24 next (i) = i - 1 ENDDO next (1) = 24 i24 = 24 j24 = 10 carry = 0. IF (seeds (24) .EQ.0.) carry = twom24 ! If restarting at a break point, skip K1 + IGIGA*K2 ! Note that this is the number of numbers delivered to ! the user PLUS the number skipped (if luxury .GT. 0). kount = k1 mkount = k2 IF (k1 + k2.NE.0) THEN DO iouter = 1, k2 + 1 inner = igiga IF (iouter.EQ.k2 + 1) inner = k1 DO isk = 1, inner uni = seeds (j24) - seeds (i24) - carry IF (uni.LT.0.) THEN uni = uni + 1.0 carry = twom24 ELSE carry = 0. ENDIF seeds (i24) = uni i24 = next (i24) j24 = next (j24) ENDDO ENDDO ! Get the right value of IN24 by direct calculation in24 = MOD (kount, nskip + 24) IF (mkount.GT.0) THEN izip = MOD (igiga, nskip + 24) izip2 = mkount * izip + in24 in24 = MOD (izip2, nskip + 24) ENDIF ! Now IN24 had better be between zero and 23 inclusive IF (in24.GT.23) THEN Write (6, '(A/A,3I11,A,I5)') ' Error in RESTARTING with RLUXGO:',& ' The values', ins, k1, k2, ' cannot occur at luxury level', luxlev in24 = 0 ENDIF ENDIF RETURN ! SUBROUTINE rluxgo END SUBROUTINE rluxgo !---------------------------------- Read in variables FUNCTION rINPUT (text) !------------------------------------------------ Character text * (*) Write (*,'(A,$)') text READ (*,*) x rINPUT = x Return END FUNCTION rINPUT !-------------------------------------- Fourier Transform SUBROUTINE SETF67 (JBC1, JQ1, jbc, jp, jsl, ll1, k2, k3, k4, k7, & Qi, jndx, Uf, Vf) ! ----------------------------------------------------- DIMENSION Uf (marr), Vf (marr) DIMENSION Qi (mf67), jndx (mf67) PI = DATAN (1.D0) * 4. JBC = JBC1 JQ = JQ1 IF (JBC.LT.3) GOTO 101 JQ = JQ - 1 101 K3 = 2**JQ K7 = K3 / 2 N5 = K3 / 4 I = 1 JNDX (I) = N5 QI (I) = 0.5 * SQRT (2.) K = I I = I + 1 102 IL = I IF (I.EQ.K7) GOTO 104 103 K1 = JNDX (K) / 2 JNDX (I) = K1 QI (I) = SIN (PI * FLOAT (K1) / FLOAT (K3) ) K1 = K7 - K1 I = I + 1 JNDX (I) = K1 QI (I) = SIN (PI * FLOAT (K1) / FLOAT (K3) ) K = K + 1 I = I + 1 IF (K.EQ.IL) GOTO 102 GOTO 103 104 RETURN END SUBROUTINE SETF67 !----------------------------------------- SUBROUTINE TFOLD (IS, L, ZZZ, k2) DIMENSION ZZZ (MARR) IH2 = K2 / 2 - 1 DO 100 I = IS, IH2 I1 = I + L I2 = K2 - I + L A = ZZZ (I1) ZZZ (I1) = A - ZZZ (I2) ZZZ (I2) = A + ZZZ (I2) 100 END DO RETURN END SUBROUTINE TFOLD !------------------------------------ SUBROUTINE NEG (I1, I3, I2, jbc, jp, jsl, ll1, k2, k3, k4, k7, Qi,& jndx, Uf, Vf) DIMENSION Uf (marr), Vf (marr) DIMENSION Qi (mf67), jndx (mf67) DO K = I1, I3, I2 Vf (K) = - Vf (K) END DO RETURN END SUBROUTINE NEG !------------------------------------- SUBROUTINE REVNEG (jbc, jp, jsl, ll1, k2, k3, k4, k7, Qi, jndx, & Uf, Vf) DIMENSION Uf (marr), Vf (marr) DIMENSION Qi (mf67), jndx (mf67) INTEGER jbc, jp, jsl, ll1, k2, k3, k4, k7 DO I = 1, K7 J = K3 + 1 + I K = K4 + 1 - I A = Vf (J) Vf (J) = - Vf (K) Vf (K) = - A END DO RETURN END SUBROUTINE REVNEG !------------------------------------- SUBROUTINE ZEERO (L, jbc, jp, jsl, ll1, k2, k3, k4, k7, Qi, jndx, & Uf, Vf) DIMENSION Uf (marr), Vf (marr) DIMENSION Qi (mf67), jndx (mf67) INTEGER jbc, jp, jsl, ll1, k2, k3, k4, k7 DO I = 1, K2 LI = L + I Uf (LI - 1) = 0.0 END DO RETURN END SUBROUTINE ZEERO !------------------------------------------ SUBROUTINE TFOLD1 (jbc, jp, jsl, ll1, k2, k3, k4, k7, Qi, jndx, & Uf, Vf) DIMENSION Uf (marr), Vf (marr) DIMENSION Qi (mf67), jndx (mf67) INTEGER jbc, jp, jsl, ll1, k2, k3, k4, k7 II = K2 - 1 DO I = 1, II I1 = JSL + I I2 = LL1 - I A = Uf (I1) Uf (I1) = A + Uf (I2) Uf (I2) = A - Uf (I2) END DO RETURN END SUBROUTINE TFOLD1 !------------------------------------- SUBROUTINE KFOLD (jbc, jp, jsl, ll1, k2, k3, k4, k7, Qi, jndx, Uf,& Vf) DIMENSION Uf (marr), Vf (marr) DIMENSION Qi (mf67), jndx (mf67) INTEGER jbc, jp, jsl, ll1, k2, k3, k4, k7 JS1 = K2 I = 1 J5 = JSL + K2 IS1 = JSL IC1 = LL1 JS1 = JS1 / 2 IF (JS1.NE.1) GOTO 200 K1 = JNDX (I) SN = QI (I) IS0 = IS1 IS1 = IS1 + JS1 IC0 = IC1 IC1 = IC1 + JS1 ODD1 = SN * (Uf (IC1) - Uf (IS1) ) Vf (K1 + 1) = Uf (IC0) + ODD1 K3MK1 = K3 - K1 Vf (K3MK1 + 1) = Uf (IC0) - ODD1 IF (JBC.LT.3) GOTO 110 ODD2 = SN * (Uf (IC1) + Uf (IS1) ) K3PK1 = K3 + K1 Vf (K3PK1 + 1) = Uf (IS0) + ODD2 K4MK1 = K4 - K1 Vf (K4MK1 + 1) = - Uf (IS0) + ODD2 110 RETURN 200 SN = QI (I) IS1 = IS1 + JS1 IC1 = IC1 + JS1 J3 = IS1 + JS1 210 IS0 = IS1 - JS1 IC0 = IC1 - JS1 ODD1 = SN * (Uf (IC1) - Uf (IS1) ) ODD2 = SN * (Uf (IC1) + Uf (IS1) ) Uf (IC1) = Uf (IC0) - ODD1 Uf (IS1) = - Uf (IS0) + ODD2 Uf (IC0) = Uf (IC0) + ODD1 Uf (IS0) = Uf (IS0) + ODD2 IS1 = IS1 + 1 IC1 = IC1 + 1 IF (IS1.NE.J3) GOTO 210 I = I + 1 300 IS1 = JSL IC1 = LL1 JS1 = JS1 / 2 IF (JS1.EQ.1) GOTO 400 310 SN = QI (I) I = I + 1 CS = QI (I) IS1 = IS1 + JS1 IC1 = IC1 + JS1 J3 = IS1 + JS1 320 IS0 = IS1 - JS1 IC0 = IC1 - JS1 ODD1 = CS * Uf (IC1) - SN * Uf (IS1) ODD2 = SN * Uf (IC1) + CS * Uf (IS1) Uf (IC1) = Uf (IC0) - ODD1 Uf (IC0) = Uf (IC0) + ODD1 Uf (IS1) = - Uf (IS0) + ODD2 Uf (IS0) = Uf (IS0) + ODD2 IS1 = IS1 + 1 IC1 = IC1 + 1 IF (IS1.NE.J3) GOTO 320 IS1 = IS1 + JS1 IC1 = IC1 + JS1 J3 = IS1 + JS1 330 IS0 = IS1 - JS1 IC0 = IC1 - JS1 ODD1 = SN * Uf (IC1) - CS * Uf (IS1) ODD2 = CS * Uf (IC1) + SN * Uf (IS1) Uf (IC1) = Uf (IC0) - ODD1 Uf (IC0) = Uf (IC0) + ODD1 Uf (IS1) = - Uf (IS0) + ODD2 Uf (IS0) = Uf (IS0) + ODD2 IS1 = IS1 + 1 IC1 = IC1 + 1 IF (IS1.NE.J3) GOTO 330 I = I + 1 IF (IS1.EQ.J5) GOTO 300 GOTO 310 400 K1 = JNDX (I) SN = QI (I) I = I + 1 CS = QI (I) IS0 = IS1 IS1 = IS1 + JS1 IC0 = IC1 IC1 = IC1 + JS1 ODD1 = CS * Uf (IC1) - SN * Uf (IS1) Vf (K1 + 1) = Uf (IC0) + ODD1 K3MK1 = K3 - K1 Vf (K3MK1 + 1) = Uf (IC0) - ODD1 IF (JBC.LT.3) GOTO 410 ODD2 = SN * Uf (IC1) + CS * Uf (IS1) K3PK1 = K3 + K1 Vf (K3PK1 + 1) = Uf (IS0) + ODD2 K4MK1 = K4 - K1 Vf (K4MK1 + 1) = - Uf (IS0) + ODD2 410 IS1 = IS1 + 1 IC1 = IC1 + 1 K1 = JNDX (I) IS0 = IS1 IS1 = IS1 + JS1 IC0 = IC1 IC1 = IC1 + JS1 ODD1 = SN * Uf (IC1) - CS * Uf (IS1) Vf (K1 + 1) = Uf (IC0) + ODD1 K3MK1 = K3 - K1 Vf (K3MK1 + 1) = Uf (IC0) - ODD1 IF (JBC.LT.3) GOTO 420 ODD2 = CS * Uf (IC1) + SN * Uf (IS1) K3PK1 = K3 + K1 Vf (K3PK1 + 1) = Uf (IS0) + ODD2 K4MK1 = K4 - K1 Vf (K4MK1 + 1) = - Uf (IS0) + ODD2 420 IS1 = IS1 + 1 IC1 = IC1 + 1 I = I + 1 IF (IS1.NE.J5) GOTO 400 RETURN END SUBROUTINE KFOLD ! ------------------------------------------------------ SUBROUTINE FOUR67 (JBC1, JQ1, jp1, jsl1, ll11, k21, k71, Qi, jndx,& Uf, Vf) ! ------------------------------------------------------ DIMENSION Uf (marr), Vf (marr) DIMENSION Qi (mf67), jndx (mf67) INTEGER jbc, jp, jsl, ll1, k2, k3, k4, k7 JBC = JBC1 JQ = JQ1 jp = jp1 jsl = jsl1 ll1 = ll11 k2 = k21 k7 = k71 A5 = 0.5 * SQRT (2.0) K4 = 2**JQ K3 = K4 GOTO (103, 103, 101, 102), JBC 101 Uf (1) = Uf (1) / 2.0 Uf (K3 + 1) = Uf (1) K2 = K3 CALL TFOLD (0, 1, Uf, k2) K3 = K3 / 2 JQ = JQ - 1 GOTO 103 102 K3 = K3 / 2 JQ = JQ - 1 103 N5 = K3 / 4 K7 = K3 / 2 N11 = 3 * K7 K31 = K3 + 1 GOTO (300, 400, 500, 600), JBC 300 Uf (K31) = 0.0 Uf (1) = 0.0 K2 = K3 DO 301 I = 2, JQ CALL TFOLD (1, 1, Uf, k2) 301 K2 = K2 / 2 Vf (K7 + 1) = Uf (2) JF = N5 JSL = 1 DO 302 JP = 2, JQ LL1 = K2 + 1 CALL ZEERO (1, jbc, jp, jsl, ll1, k2, k3, k4, k7, Qi, jndx, Uf,& Vf) CALL KFOLD (jbc, jp, jsl, ll1, k2, k3, k4, k7, Qi, jndx, Uf, & Vf) I1 = 3 * JF + 1 I2 = 4 * JF I3 = I1 + (K2 / 2 - 1) * I2 CALL NEG (I1, I3, I2, jbc, jp, jsl, ll1, k2, k3, k4, k7, Qi, & jndx, Uf, Vf) K2 = K2 + K2 302 JF = JF / 2 RETURN 400 Uf (1) = Uf (1) / 2.0 Uf (K31) = Uf (K31) / 2.0 K2 = K3 DO 401 I = 2, JQ CALL TFOLD (0, K31 - K2, Uf, k2) 401 K2 = K2 / 2 LL1 = K31 - K2 A = Uf (LL1) + Uf (LL1 + 2) Vf (1) = - A - Uf (LL1 + 1) Vf (K31) = - A + Uf (LL1 + 1) Vf (K7 + 1) = Uf (LL1) - Uf (LL1 + 2) DO 402 JP = 2, JQ JSL = K31 - K2 LL1 = JSL - K2 idum = jsl CALL ZEERO (idum, jbc, jp, jsl, ll1, k2, k3, k4, k7, Qi, jndx, & Uf, Vf) CALL KFOLD (jbc, jp, jsl, ll1, k2, k3, k4, k7, Qi, jndx, Uf, & Vf) 402 K2 = K2 + K2 CALL NEG (1, K31, 2, jbc, jp, jsl, ll1, k2, k3, k4, k7, Qi, jndx, & Uf, Vf) RETURN 500 K2 = K3 L2 = K4 DO 501 JP = 2, JQ CALL TFOLD (1, 1, Uf, k2) CALL TFOLD (0, L2 - K2 + 1, Uf, k2) 501 K2 = K2 / 2 LL1 = L2 - K2 + 1 A = Uf (LL1) + Uf (LL1 + 2) Vf (K7 + 1) = 2.0 * ( - Uf (LL1) + Uf (LL1 + 2) ) Vf (1) = 2.0 * (A + Uf (LL1 + 1) ) Vf (K31) = 2.0 * (A - Uf (LL1 + 1) ) Vf (N11 + 1) = 2.0 * Uf (2) DO 502 JP = 2, JQ Uf (K2 + 1) = 2.0 * Uf (K2 + 1) JSL = K2 + 1 CALL TFOLD1 (jbc, jp, jsl, ll1, k2, k3, k4, k7, Qi, jndx, Uf, & Vf) LL1 = LL1 - K2 Uf (LL1) = - 2.0 * Uf (LL1) CALL KFOLD (jbc, jp, jsl, ll1, k2, k3, k4, k7, Qi, jndx, Uf, & Vf) 502 K2 = K2 + K2 Vf (1) = Vf (1) * A5 Vf (K31) = Vf (K31) * A5 RETURN 600 Uf (1) = Uf (1) * A5 Uf (K31) = Uf (K31) * A5 K2 = K3 DO 601 JP = 2, JQ CALL TFOLD (0, K31 - K2, Uf, k2) CALL TFOLD (1, K31, Uf, k2) 601 K2 = K2 / 2 LL1 = K31 - K2 A = Uf (LL1) + Uf (LL1 + 2) Vf (1) = 2.0 * (A + Uf (LL1 + 1) ) Vf (K31) = 2.0 * (A - Uf (LL1 + 1) ) Vf (K7 + 1) = 2.0 * ( - Uf (LL1) + Uf (LL1 + 2) ) Vf (N11 + 1) = 2.0 * Uf (K3 + 2) DO 602 JP = 2, JQ JSL = K31 + K2 Uf (JSL) = 2.0 * Uf (JSL) CALL TFOLD1 (jbc, jp, jsl, ll1, k2, k3, k4, k7, Qi, jndx, Uf, & Vf) LL1 = LL1 - K2 Uf (LL1) = - 2.0 * Uf (LL1) CALL KFOLD (jbc, jp, jsl, ll1, k2, k3, k4, k7, Qi, jndx, Uf, & Vf) 602 K2 = K2 + K2 CALL REVNEG (jbc, jp, jsl, ll1, k2, k3, k4, k7, Qi, jndx, Uf, Vf) idum = k3 CALL NEG (2, idum, 2, jbc, jp, jsl, ll1, k2, k3, k4, k7, Qi, jndx,& Uf, Vf) K2 = K4 CALL TFOLD (1, 1, Vf, k2) Vf (K4 + 1) = 0.0 RETURN END SUBROUTINE FOUR67 end module setInitialConditions ! ! BLOCKDATA PARAMETER ( maxlev = 4, lxdflt = 3) LOGICAL notyet COMMON /LUX1/ & i24, j24, in24, kount, mkount, carry, seeds(24) COMMON /LUX2/ notyet, twom24, twom12, & nskip, ndskip(0:maxlev), next(24), inseed, luxlev DATA ndskip / 0, 24, 73, 199, 365 / DATA i24, j24, luxlev / 24, 10, lxdflt / DATA notyet / .true. / DATA in24, kount, mkount, carry / 0, 0, 0, 0. / End BLOCKDATA !------------------------------------------------------------------------- PROGRAM PMstartMp use setInitialConditions use mpi ! NROW = number of particles in 1 dimension ! NGRID= number of cells in 1 dimension ! Npage = number of particles per page = NROW**2 ! COMMON / TRUNCOM / Om, Omb, Omc, Omnu, Par (6), ns, qqscaleb,QSCALE DOUBLE PRECISION :: Wtotal=0., SKINE=0. REAL, PARAMETER :: PI=3.1415926535 REAL :: t0=0,t1=0. ! counters for timing Logical :: iCont =.true., ext ! continue (true) after setting seeds Integer*8 :: icurrent character*80 fName CALL MPI_INIT(ierr) ! Initialize MPI CALL MPI_COMM_RANK(MPI_COMM_WORLD,rank,ierr) CALL MPI_COMM_SIZE(MPI_COMM_WORLD, numprocs,ierr) myrank = rank +1 lchunk =8 ! length of a chunk of MPI tasks for disk IO nchunks = (numprocs-1)/lchunk+1 mychunk = (myrank-1)/lchunk+1 AEXPN = 0. ! Do m=1,nchunks ! create output file for every thread ! CALL MPI_Barrier(MPI_COMM_WORLD,ierr) ! If(mychunk == m)Then If(myrank < 33)Then write(fName,'(a,i4.4,a)') 'OutputRank.',myrank,'.dat' open(16,file=fName) write(16,*)' Init:',myrank EndIf ! EndDo ! Initialize ------------------------------------------ CALL CPU_TIME(t0) If(rank == 0) CALL InitValues (SCALEL) CALL MPI_Barrier(MPI_COMM_WORLD,ierr) ! CALL MPI_Finalize(ierr) ! stop if(myrank < 33)write(16,*) ' go to BcastInit A=',AEXPN CALL Bcast_InitValues(SCALEL) ! CALL MPI_Finalize(ierr) ! stop if(myrank < 33)Write (16,*) ' InitValues is done' Max_lev_abs = INT (ALOG (REAL (Lblock) ) / ALOG (2.) + 1.1) if(myrank < 33)Write (16,'(10x,a,2i4)') 'Absolute Limits on mass resolution levels: ', 1,& Max_lev_abs if(myrank < 33)Write (16,'(10x,2(a,i4))') 'Actual Max resolution Level=', Nmax_lev, & ' Min Level=', Nmin_lev ! If (Nmax_lev.gt.Max_lev_abs) Then ! Write (16,*) ' Wrong number of Levels for mass refinement (', & ! Nmax_lev, ') Block size =', Lblock, ' Use Absolute value=', & ! Max_lev_abs ! Nmax_lev = Max_lev_abs ! Endif ! Write (16, * ) ' Seed=', Nseed AEXP0 = AEXPN AEXPV = AEXPN - ASTEP / 2./2**LevelParticles INQUIRE(FILE='Allspecies.DAT',EXIST=ext) If(ext)Then ! files exist. Go to gathering part. Skip initialization If(rank == 0)Then open(70,file='Allspecies.DAT') read (70,*) Nspecies read(70,*)(nLevel(i),i=1,Nspecies) Do j=1,NROW read(70,*)(Allspecies(i,j),i=1,Nspecies) EndDo write(16,*) ' ------------- Reading Allspecies file. Skip Spectrum and FFT ---' close (70) EndIf ! rank ==0 go to 10 EndIf ! exist CALL SetRandom(iCont) ! set seeds for ramdom numbers ! CALL MPI_Finalize(ierr) ! stop !' --------------' CALL ZeroMASK ! set mask for multiple mass resolution If (Nmax_lev.ne.Nmin_lev) Then Open (30, file = 'SelCoords.DAT', status = 'unknown') CALL SetMASK EndIf ! CALL MPI_Finalize(ierr) ! stop ' --------------' ! set scaling constants for coordinates and velocities IFOUR = INT (ALOG (FLOAT (NROW) ) / ALOG (2.) + 0.5) !Write (*,'(3x,a12,i4,a,i4)') ' Ifour =', IFOUR, ' Nrow=', NROW !Write(16,'(a,2i8)' ) ' Before Allocate ',NROW,NPAGE !CALL MPI_Finalize(ierr) !stop ' --------------' Fact = sqrt (Om0 + Oml0 * AEXPV**3 + Ocurv * AEXPV) QFACT = FLOAT (NGRID) / FLOAT (NROW) Vscale = ScaleL * 100. * hubble / NGRID CALL CPU_TIME(t1) timing(1) = t1-t0 ! time for intialization If(.NOT.iCont)Then write(*,'(2(10x,a,f10.2/))') & 'cpu_time for init =',timing(1), & 'cpu_time for seeds=',timing(2) Stop EndIf !Write(16,'(a,2i8)' ) ' Before Allocate ',NROW,NPAGE !CALL MPI_Finalize(ierr) !stop ' --------------' ALLOCATE(GR(NROW,NROW)) ! allocate memory for spectrum Do iDirection = 1, 3 ! ----------------- loop over x,y,z directions NRAND = Nseed if(myrank < 33)Write(16,'(//3x,a12,i3,15("-"))' )'- Direction=',iDirection !If(iDirection == 3)Then ! CALL MPI_Finalize(ierr) ! stop ' --------------' !EndIf CALL SPECTR (NRAND, iDirection) !If(iDirection == 3)Then ! CALL MPI_Finalize(ierr) ! stop 'MM --------------' !EndIf VCONS = - ALPHA/(2.*PI/NGRID)*(AEXPV/AEXP0)*SQRT(AEXPV)*Fact XCONS = ALPHA / (2. * PI / NGRID) * (AEXPN / AEXP0) if(myrank < 33)Write(16,'(3x,a12,g12.4,a,g12.4)' ) 'Scaling:(x)=', XCONS, ' (v)=', VCONS CALL VECTOR (IFOUR) ! get the displacement vector by FFT if(myrank < 33)Write (16,'(3x,a12,i4,a,g12.4)' ) 'SPECTR done:', IFOUR, ' Alpha=', ALPHA if(myrank < 33)write(16,'(2(a,f9.1))')' Spectr t=',timing(3),' FFT=',timing(5) ! find x,v and write to file, page by page !CALL MPI_Finalize(ierr) !stop ' 2--------------' lspecies = 0 CALL BLOCKS (XCONS, VCONS, iDirection) !CALL MPI_Finalize(ierr) !stop ' --------------' EndDo DEALLOCATE (GR) EKIN = 0.5 * SKINE / AEXPV**2 !Write (16, '('' Ekin='',E12.4,'' Weght per cell='',g12.5,'' (must be 1)'')') EKIN,Wtotal/NGRID**3 ! CALL MPI_Finalize(ierr) ! stop ' --------------' 10 CALL FilesWrite ! gather particles and write PM files !CALL CPU_TIME(t1) ! write(*,'(10x,a/,(10x,a15,f10.2))') & ! 'Wall cpu_time (seconds) ', & ! 'Init =',timing(1), & ! 'Seeds =',timing(2), & ! 'Spectr =',timing(3), & ! 'Coords =',timing(4), & ! 'FFT =',timing(5), & ! 'Write temp =',timing(6), & ! 'Write final =',timing(7), & ! 'Total =',t1-t0 CALL MPI_Finalize(ierr) STOP END PROGRAM PMstartMp