-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathpowerspectrum.f90
119 lines (103 loc) · 3.45 KB
/
powerspectrum.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
!
!> Compute the power spectrum from an L-Gadget snapshot file.
!
! Created by Michael Schneider on 2010-04-09
!
program powerspectrum
use types_nsample
use parameters
use io_nsample
use fft_nsample
use gridtools, only: part2grid, print_gridarray
use resample, only: npkbins, power_spectrum_estimator
implicit none
character(len=200) :: arg, outdir, outfile
real(dp), dimension(:,:), allocatable :: coor, incoor
integer(i8b), dimension(npart) :: pid
real(dp) :: scale_factor
integer(i4b) :: ng, i,j, isnap
real(dp), dimension(npk,2) :: pk
real(dp), dimension(:), allocatable :: ps
real(dp), dimension(:,:,:), allocatable :: delta
character(len=10) :: isnaplab
logical :: icfile
#ifdef MPI
! ----- Initialize MPI
call mpi_init(ierr)
call mpi_comm_rank(mpi_comm_world,rank,ierr)
call mpi_comm_size(mpi_comm_world,numtasks,ierr)
print *, "Process ",rank," of ",numtasks," is alive"
#else
rank = 0
numtasks = 1
#endif MPI
icfile = .false.
! ----- Parse command line arguments
if (rank .eq. 0) then
if (command_argument_count() /= 3) then
write(6,*) "Usage: powerspectrum.x ng isnap outdir"
stop
end if
call get_command_argument(1,arg)
read(arg,*) ng
call get_command_argument(2,arg)
read(arg,*) isnap
call get_command_argument(3,outdir)
end if
#ifdef MPI
call mpi_bcast(ng, 1, mpi_integer, 0, mpi_comm_world, ierr)
call mpi_bcast(isnap, 1, mpi_integer, 0, mpi_comm_world, ierr)
call mpi_bcast(outdir, 200, mpi_character, 0, mpi_comm_world, ierr)
#endif
! ----- Allocate output power spectrum array
allocate(ps(0:npkbins(ng)))
! ----- Initialize FFTW plans and MPI slab decomposition
#ifdef MPI
call init_fftw(rank, numtasks, ng, ng)
if (rank .eq. 0) print *, "Finished FFTW init"
#else !MPI
local_nlast = ng
local_last_start = 0
local_nlastL = ng
local_last_startL = 0
#endif MPI
call parameters_init()
! ----- Read snapshot and distribute particles
allocate(coor(3, npart/numtasks))
if (rank .eq. 0) then
allocate(incoor(3,npart))
call read_snapshot(outdir, isnap, incoor, pid, scale_factor, icfile=icfile)
print *, "scale factor of input snapshot: ",scale_factor
end if
#ifdef MPI
call mpi_bcast(scale_factor, 1, mpi_double_precision, 0, mpi_comm_world, ierr)
call mpi_scatter(incoor, 3*npart/numtasks, mpi_double_precision, &
& coor, 3*npart/numtasks, mpi_double_precision, &
& 0, mpi_comm_world, ierr)
call mpi_barrier(mpi_comm_world, ierr)
#else
coor = incoor
#endif MPI
if (rank .eq. 0) deallocate(incoor)
allocate(delta(ng+fftpad, ng, local_nlast))
call part2grid(ng, npart, boxsize, coor(:,1:npart/numtasks), delta)
deallocate(coor)
call print_gridarray(delta, "snapshot density")
! ----- Compute power spectrum
write(isnaplab,'(i3.3)') isnap
if (pkbintype .eq. 1) then
outfile = trim(outdir)//"/snapdir_"//trim(adjustl(isnaplab))// &
& "/ps_takbins_"//trim(adjustl(isnaplab))//".txt"
else
if (icfile) then
outfile = trim(outdir)// &
& "/ps_"//trim(adjustl(isnaplab))//".txt"
else
outfile = trim(outdir)//"/snapdir_"//trim(adjustl(isnaplab))// &
& "/ps_"//trim(adjustl(isnaplab))//".txt"
end if
end if
call power_spectrum_estimator(ng, ps, delta, 1, 1, outfile)
deallocate(delta, ps)
stop
end program powerspectrum