DISPATCH
pdf_io_mod.f90
1 !===============================================================================
2 !> Module to periodically output a density probability density function (PDF).
3 !> The procedure assumes that 'out_time' is larger than the largest time step.
4 !===============================================================================
5 MODULE pdf_io_mod
6  USE io_mod
7  USE io_unit_mod
8  USE trace_mod
9  USE patch_mod
10  USE link_mod
11  USE task_mod
12  USE counters_mod
13  USE mpi_mod
14  USE bits_mod
15  implicit none
16  private
17  !.............................................................................
18  integer, parameter:: n_rms=6
19  type, public:: pdf_io_t
20  logical:: on=.false.
21  integer:: iout=0
22  real(8):: time, out_time
23  real:: rms(n_rms)=0.0
24  type(counters_t):: counters
25  contains
26  procedure:: init
27  procedure:: update
28  procedure:: output
29  end type
30  !.............................................................................
31  real, dimension(:), allocatable:: bins, counts
32  real(8):: out_time=0.1_8
33  real:: bin_min=-12.0, bin_max=10.0, bin_delta
34  integer:: nbins=111
35  integer:: verbose=0
36  logical:: on=.false.
37  type(pdf_io_t), public:: pdf_io
38 CONTAINS
39 
40 !===============================================================================
41 !> Initialize parameters for PDF output
42 !===============================================================================
43 SUBROUTINE init (self, patch)
44  class(pdf_io_t):: self
45  class(patch_t):: patch
46  !.............................................................................
47  integer:: i, iostat
48  logical:: first_time=.true.
49  namelist /pdf_io_params/ on, verbose, out_time, nbins, bin_min, bin_max
50  !-----------------------------------------------------------------------------
51  call trace%begin ('pdf_io_t%init')
52  if (first_time) then
53  !$omp critical (input_cr)
54  if (first_time) then
55  rewind(io%input)
56  read (io%input, pdf_io_params, iostat=iostat)
57  write (io%output, pdf_io_params)
58  if (on) then
59  allocate (bins(nbins), counts(nbins))
60  bin_delta = (bin_max-bin_min)/(nbins-1.0)
61  do i=1,nbins
62  bins(i) = bin_min + (i-1)*bin_delta
63  counts(i) = 0.0
64  end do
65  end if
66  first_time = .false.
67  end if
68  !$omp end critical (input_cr)
69  end if
70  self%on = on
71  if (on) then
72  self%out_time = out_time
73  if (patch%pdf_next==0d0) &
74  patch%pdf_next = out_time
75  call self%counters%init
76  end if
77  call trace%end()
78 END SUBROUTINE init
79 
80 !===============================================================================
81 !> PDF update, to be called when a patch first passes out_next
82 !===============================================================================
83 SUBROUTINE update (self, patch)
84  class(pdf_io_t):: self
85  class(patch_t):: patch
86  !.............................................................................
87  character(len=80):: filename
88  integer:: ix, iy, iz, l(3), u(3), i, ibin, jt(2), iout, count
89  real:: dv, bin, pt(2), d, px, py, pz, ux, uy, uz, u2, rms(n_rms)
90  real(8):: v_tot, d_tot, px_tot, py_tot, pz_tot, u2_tot
91  logical:: deallocated
92  !-----------------------------------------------------------------------------
93  if (.not.on) return
94  if (patch%is_set(bits%remove)) then
95  write (stdout,*) 'WARNING: tried to count PDF of removed', patch%id
96  return
97  end if
98  if (patch%is_set(bits%virtual)) then
99  write (stderr,*) 'ERROR: tried to count PDF of virtual', patch%id
100  return
101  end if
102  call trace%begin ('pdf_io_t%update')
103  if (patch%time > patch%pdf_next) then
104  call patch%time_interval (patch%pdf_next, jt, pt)
105  deallocated = .not.allocated(patch%irefine)
106  if (deallocated) &
107  allocate (patch%irefine(patch%gn(1),patch%gn(2),patch%gn(3)))
108  call make_irefine (patch)
109  dv = product(patch%mesh%d)
110  d_tot = 0.0
111  px_tot = 0.0
112  py_tot = 0.0
113  pz_tot = 0.0
114  u2_tot = 0.0
115  v_tot = 0.0
116  l = patch%mesh%li
117  u = patch%mesh%ui
118  do iz=l(3),u(3)
119  do iy=l(2),u(2)
120  do ix=l(1),u(1)
121  if (patch%irefine(ix,iy,iz) == 0) then
122  d = patch%mem(ix,iy,iz,patch%idx%d ,jt(1),1)*pt(1) &
123  + patch%mem(ix,iy,iz,patch%idx%d ,jt(2),1)*pt(2)
124  px = patch%mem(ix,iy,iz,patch%idx%px,jt(1),1)*pt(1) &
125  + patch%mem(ix,iy,iz,patch%idx%px,jt(2),1)*pt(2)
126  py = patch%mem(ix,iy,iz,patch%idx%py,jt(1),1)*pt(1) &
127  + patch%mem(ix,iy,iz,patch%idx%py,jt(2),1)*pt(2)
128  pz = patch%mem(ix,iy,iz,patch%idx%pz,jt(1),1)*pt(1) &
129  + patch%mem(ix,iy,iz,patch%idx%pz,jt(2),1)*pt(2)
130  ux = px/d
131  uy = py/d
132  uz = pz/d
133  u2 = ux**2 + uy**2 + uz**2
134  px_tot = px_tot + px*dv
135  py_tot = py_tot + py*dv
136  pz_tot = pz_tot + pz*dv
137  u2_tot = u2_tot + u2*dv
138  d_tot = d_tot + d*dv
139  v_tot = v_tot + dv
140  bin = (alog(d) - bin_min)/bin_delta
141  ibin = 1.0+bin
142  ibin = max(1,min(ibin,nbins))
143  counts(ibin) = counts(ibin) + dv
144  end if
145  end do
146  end do
147  end do
148  rms = [v_tot, d_tot, px_tot, py_tot, pz_tot, u2_tot]
149  do i=1,n_rms
150  !$omp atomic update
151  self%rms(i) = self%rms(i) + rms(i)
152  end do
153  if (deallocated) &
154  deallocate (patch%irefine)
155  !---------------------------------------------------------------------------
156  ! Decrement the count of tasks still remaining to be accumulated. iout is
157  ! the output index (0 at time=0.0) and the counter index, which needs to be
158  ! non-zero, is iout+1. The first task that passes a give pdf_next refers to
159  ! a non-existing counter, which causes the counter to be created, with initial
160  ! count = io%nwrite-1.
161  !---------------------------------------------------------------------------
162  !$omp critical (update_counters_cr)
163  iout = nint(patch%pdf_next/pdf_io%out_time)
164  call self%counters%update (iout+1, io%nwrite, -1, count)
165  if (verbose > 0) &
166  write(stdout,'(a,i6,3f12.6,4i4,i5,2x,a)') &
167  'pdf_io_t%update: id, time, pdf_next =', patch%id, patch%time, &
168  patch%time-patch%dtime, patch%pdf_next, count, self%iout, iout, &
169  patch%level, patch%n_nbors, 'output'
170  !---------------------------------------------------------------------------
171  ! Call the output procedure when all tasks have been accumulated. All tasks
172  ! have their pdf_next incremented after their data are accumulated.
173  !---------------------------------------------------------------------------
174  if (count==0) then
175  self%time = patch%pdf_next
176  call self%output
177  end if
178  patch%pdf_next = (iout+1)*out_time
179  !$omp end critical (update_counters_cr)
180  end if
181  call trace%end()
182 END SUBROUTINE update
183 
184 !===============================================================================
185 !> PDF output, delayed until all patches have been processed. The PDF snapshot
186 !> number is kept as pdf_io%iout, and since we only use one counts(:) array, we
187 !> need to assume that this is the same for all patches.
188 !===============================================================================
189 SUBROUTINE output (self)
190  class(pdf_io_t):: self
191  !.............................................................................
192  character(len=80):: file
193  integer:: pdf_iout
194  integer, parameter:: ioformat=1
195  !-----------------------------------------------------------------------------
196  call trace%begin ('pdf_io_t%output')
197  !---------------------------------------------------------------------------
198  ! Increment the output counter, keeping track of the no of snapshots written.
199  ! The counter corresponding to the pending snapshot is (or has been) created
200  ! when the first task sets pdf_next to that output time.
201  !---------------------------------------------------------------------------
202  self%iout = self%iout+1
203  write (file,'(a,"pdf_",i5.5,"_",i5.5,".dat")') &
204  trim(io%outputname), self%iout, mpi%rank
205  open (io_unit%tmp, file=file, form='unformatted', status='unknown')
206  write (io_unit%tmp) ioformat, nbins
207  write (io_unit%tmp) self%time
208  write (io_unit%tmp) bins
209  write (io_unit%tmp) counts
210  write (io_unit%tmp) self%rms
211  close (io_unit%tmp)
212  if (verbose > -1) &
213  write (stdout,*) trim(file)
214  counts(:) = 0.0
215  self%rms(:) = 0.0
216  call trace%end()
217 END SUBROUTINE output
218 
219 !===============================================================================
220 !> Make a map of where the patch is refined
221 !===============================================================================
222 SUBROUTINE make_irefine (patch)
223  class(patch_t):: patch
224  class(link_t), pointer:: nbor
225  class(task_t), pointer:: nbpatch
226  integer:: l(3), u(3)
227  real(8):: dist(3)
228  !-----------------------------------------------------------------------------
229  patch%irefine = 0
230  nbor => patch%link%nbor
231  do while (associated(nbor))
232  nbpatch => nbor%task
233  if (nbpatch%level > patch%level+1) then
234  write (stderr,*) 'WARNING: nbpatch w/o support =', nbpatch%id, patch%id
235  else if (nbpatch%level > patch%level) then
236  select type (nbpatch)
237  class is (patch_t)
238  if (patch%contains(nbpatch)) then
239  dist = (nbpatch%size-nbpatch%ds)/2d0
240  l = patch%index_only (nbpatch%position-dist)
241  u = patch%index_only (nbpatch%position+dist)
242  if (verbose > 1) &
243  write(stdout,*) 'irefine: nbpatch =', &
244  nbpatch%id, nbpatch%level, l, u, product(u-l+1)
245  patch%irefine(l(1):u(1),l(2):u(2),l(3):u(3)) = 1
246  end if
247  end select
248  end if
249  nbor => nbor%next
250  end do
251  if (verbose > 0) then
252  l = patch%li
253  u = patch%ui
254  write(stdout,*) 'make_irefine: unrefined =', &
255  sum(1-patch%irefine(l(1):u(1),l(2):u(2),l(3):u(3))), product(u-l+1)
256  end if
257 END SUBROUTINE make_irefine
258 
259 END MODULE pdf_io_mod
Help keep track of when all patches have passed some counter, by decrementing a counter, from start to 0. Typical use:
Template module for patches, which adds pointers to memory and mesh, and number of dimensions and var...
Definition: patch_mod.f90:6
Module to periodically output a density probability density function (PDF). The procedure assumes tha...
Definition: pdf_io_mod.f90:5
Definition: io_mod.f90:4
Template module for tasks.
Definition: task_mod.f90:4