DISPATCH
parallel_io_mod.f90
1 !===============================================================================
2 !> Use MPI parallel I/O to write everything to a single file. No critical regions
3 !> should be needed here; all thread protection is done inside mpi_file_mod and
4 !> mpi_io_mod.
5 !>
6 !> There are two pairs of procedures; output/input and output_single/input_single.
7 !> The latter group supports writing each snapshot in a separate file, but only
8 !> for MPI libraries that support multiple thread, simultenous MPI calls.
9 !>
10 !> In each of the pairs, the io%format value chooses between an older storage
11 !> pattern, where the position of a patch in the disk image is always the same,
12 !> and independent of MPI configuration, and one where each rank writes a
13 !> contiguous part in the file, allowing a much faster output, and also a much
14 !> faster reading from IDL and Python.
15 !>
16 !> A compromise storage pattern would be one where the data is both stored in an
17 !> MPI-independent way, but still can be read in large chunks. A good arrangement
18 !> would have all patches in a variable in a contiguous disk space, arranged in
19 !> an order that does not depend on the MPI-configuration, nor depends on the
20 !> order of writing (there should be no need to get the order of writes from the
21 !> order of patches in the patch_rrrrr.nml files, for example).
22 !>
23 !> In the current task-ID scheme, task IDs are guaranteed to be unique, without
24 !> requiring MPI-communication when getting a new ID (which would have to be
25 !> arranged via MPI-RMA -- with indications of not being reliable). Hence, in
26 !> the current scheme, task-IDs are already MPI-dependent, and needs to be
27 !> regenerated when restarting with a different MPI-configuration. It would
28 !> thus be best if the position of a patch in the file only depende on the
29 !> position of the patch in space, and in such a way that a consequence was to
30 !> have most (if not all) patches written from a single rank written to the
31 !> same part of the file.
32 !>
33 !> A better approach, much easier to implement, is that the reader uses the
34 !> meta-data saved with the snapshot to deduce the arrangement in the file,
35 !> and reads the data based on that. It is then free to choose the pattern to
36 !> write data with based only on efficiency of writing, and speed of reading.
37 !===============================================================================
39  USE patch_mod
40  USE io_mod
41  USE io_unit_mod
42  USE omp_mod
43  USE omp_timer_mod
44  USE mpi_mod
45  USE mpi_coords_mod
46  USE mpi_io_mod
47  USE mpi_file_mod
48  USE trace_mod
49  USE time_slices_mod
50  USE bits_mod
51  USE kinds_mod
52  implicit none
53  private
54  type, public:: parallel_io_t
55  real, dimension(:,:,:,:,:,:), pointer:: buffer => null()
56  contains
57  procedure, nopass:: init
58  procedure, nopass:: input
59  procedure, nopass:: output
60  procedure, nopass:: close
61  procedure, nopass:: open_and_close
62  end type
63  integer, save:: verbose=0
64  logical:: check_files=.false.
65  type(mpi_file_t), pointer:: file_in=>null(), file_out=>null()
66  type(mpi_file_t), pointer, dimension(:):: files_out=>null()
67  type(parallel_io_t), public:: parallel_io
68  integer, save:: io0, io1
69  integer, save:: opened=-1, closed=-1
70 CONTAINS
71 
72 !===============================================================================
73 !> Initialize for MPI parallel I/O. This is an MPI collective operation, and
74 !> needs to be done at the start of the run.
75 !===============================================================================
76 SUBROUTINE init (patch)
77  class(patch_t):: patch
78  !.............................................................................
79  character(len=5):: snapname
80  character(len=120):: filename
81  logical:: exists
82  integer:: i, iostat
83  logical:: first_time=.true.
84  namelist /parallel_io_params/ verbose
85  character(len=120):: id = &
86  '$Id: d936aa6f2f05c80e7a6feced2da2dce6d4019b39 $ io/parallel_io_mod.f90'
87  !-----------------------------------------------------------------------------
88  call trace%print_id (id)
89  call trace%begin ('parallel_io_t%init')
90  !$omp critical (input_cr)
91  if (first_time) then
92  first_time = .false.
93  rewind(io%input); read(io%input, parallel_io_params, iostat=iostat)
94  write (io%output, parallel_io_params)
95  if (io%ntask==0) &
96  call io%abort ('parallel_io_t%init: io%ntask == 0')
97  call mpi_io%init
98  call mpi_io%set (nwrite=io%nwrite*(io%time_derivs+1))
99  if (io%method=='parallel') then
100  allocate (file_out)
101  write (filename,'(a,"/snapshots.dat")') trim(io%outputname)
102  call file_out%openw (filename)
103  end if
104  end if
105  !$omp end critical (input_cr)
106  call time_slices%init (patch%nt)
107  call trace%end()
108 END SUBROUTINE init
109 
110 !===============================================================================
111 !> Close the data I/O method
112 !===============================================================================
113 SUBROUTINE close
114  integer:: i
115  !.............................................................................
116  call trace%begin ('parallel_io_t%close')
117  if (io%method == 'snapshot') then
118  do i=io0,io1
119  call files_out(i)%close
120  end do
121  else
122  call file_out%close
123  end if
124  call file_in%close
125  call trace%end()
126 END SUBROUTINE close
127 
128 !===============================================================================
129 !> Adjust the open / closed status of the data I/O files
130 !===============================================================================
131 SUBROUTINE open_and_close ()
132  integer:: i
133  !-----------------------------------------------------------------------------
134  ! Only check files when a request to open or close has just been made. Only
135  ! one thread should do the checking.
136  !-----------------------------------------------------------------------------
137  if (check_files) then
138  !$omp critical (files_out_cr)
139  if (check_files) then
140  !$omp atomic write
141  check_files = .false.
142  !$omp end atomic
143  call trace%begin ('parallel_io_t%open_and_close')
144  if (associated(files_out)) then
145  do i=io0,io1
146  call files_out(i)%open_and_close()
147  end do
148  end if
149  call trace%end()
150  end if
151  !$omp end critical (files_out_cr)
152  if (verbose > 1) then
153  write (io_unit%log,*) wallclock(), 'parallel_io_t%open_and_close: exiting'
154  flush (io_unit%log)
155  end if
156  end if
157 END SUBROUTINE open_and_close
158 
159 !===============================================================================
160 !> Write output in increments of one patch data cube
161 !===============================================================================
162 SUBROUTINE output (patch, count)
163  class(patch_t):: patch
164  integer:: count
165  !.............................................................................
166  integer(8):: snapshot_words, buffer_words, pos
167  integer:: iv, n(3), jt(2), iw, nw
168  real(kind=KindScalarVar), dimension(:,:,:,:), pointer:: mem
169  real(kind=4), dimension(:,:,:,:), pointer:: buffer
170  character(len=120):: filename
171  real:: pt(2)
172  logical, save:: first_time=.true.
173  integer, save:: itimer=0
174  !-----------------------------------------------------------------------------
175  call trace%begin ('parallel_io_t%output', itimer=itimer)
176  if (io%format > 9 .and. io%format < 14) then
177  call output_single (patch, count)
178  call trace%end (itimer)
179  return
180  end if
181  !-----------------------------------------------------------------------------
182  ! Initialize file_out
183  !-----------------------------------------------------------------------------
184  if (.not.associated(file_out)) then
185  write (filename,'(a,"/snapshots.dat")') trim(io%outputname)
186  allocate (file_out)
187  call file_out%openw (filename)
188  end if
189  !-----------------------------------------------------------------------------
190  ! Set size parameters; a buffer is all variables in the patch, and the
191  ! snspshot offset is that times the number of patches per snapshot.
192  !-----------------------------------------------------------------------------
193  if (io%guard_zones) then
194  n = patch%gn
195  else
196  n = patch%n
197  end if
198  !-----------------------------------------------------------------------------
199  ! If 9 < io%format < 14 then slots in the file are ordered (id,iv,iout),
200  ! while otherwise they are ordered (iv,id,iout). The former is much faster
201  ! to read, while the latter is faster to write.
202  !-----------------------------------------------------------------------------
203  snapshot_words = product(n)*io%nv*io%ntotal
204  buffer_words = product(n)*io%nv
205  if (io%method == 'snapshot') then
206  call mpi_io%use (files_out(patch%iout))
207  else
208  call mpi_io%use (file_out)
209  end if
210  call mpi_io%set (snapshot_words, buffer_words, io%nwrite*(io%time_derivs+1))
211  !-----------------------------------------------------------------------------
212  ! Info
213  !-----------------------------------------------------------------------------
214  if (first_time) then
215  first_time = .false.
216  write(io%output,'(a,f8.3,a,i10,a)') &
217  ' output snapshot size =', snapshot_words*4d0/1024d0**3, ' GB, in', &
218  io%ntotal, ' patches'
219  write(io%output,'(a,4i6,2i12)') &
220  ' parallel_io_t%output: iout, n, snapshot_words, buffer_words =', &
221  patch%iout, n, snapshot_words, buffer_words
222  end if
223  !-----------------------------------------------------------------------------
224  ! Time interpolation to exact out_next
225  !-----------------------------------------------------------------------------
226  allocate (mem(n(1),n(2),n(3),patch%nt))
227  allocate (buffer(n(1),n(2),n(3),io%nv))
228  do iv=1,io%nv
229  call convert_variables ! copy `patch%mem` into `mem`, converting if desired
230  call time_slices%interpolate (patch, mem, buffer(:,:,:,iv))
231  end do
232  call io_write ((io%time_derivs+1)*patch%iout, 0)
233  !-----------------------------------------------------------------------------
234  ! On-the-fly time derivatives output
235  !-----------------------------------------------------------------------------
236  if (io%time_derivs>0) then
237  do iv=1,io%nv
238  call convert_variables ! copy `patch%mem` into `mem`, converting if desired
239  call time_slices%derivative1 (patch, mem, buffer(:,:,:,iv))
240  end do
241  call io_write ((io%time_derivs+1)*patch%iout, 1)
242  if (io%time_derivs==2) then
243  do iv=1,io%nv
244  call convert_variables ! copy `patch%mem` into `mem`, converting if desired
245  call time_slices%derivative2 (patch, mem, buffer(:,:,:,iv))
246  end do
247  call io_write ((io%time_derivs+1)*patch%iout, 2)
248  end if
249  end if
250  !-----------------------------------------------------------------------------
251  ! Deallocate mem, but do NOT deallocate buffer. Since io_write uses non-
252  ! blocking MPI_File_iwrite calls, the buffer must remain allocated, until a
253  ! later call to check_io finds that it is no longer needed, and deallocates it
254  !-----------------------------------------------------------------------------
255  deallocate (mem)
256  call trace%end (itimer)
257 contains
258  !-----------------------------------------------------------------------------
259  ! Copy variables from `patch%mem` to `mem` array, converting variables based
260  ! on the I/O format parameter. This procedure is internal to the 'output'
261  ! procedure, and as such only supports io%format in the range 6-9.
262  !-----------------------------------------------------------------------------
263  subroutine convert_variables
264  integer:: nt1, l(3), u(3), i, j, k
265  !...........................................................................
266  call trace%begin ('parallel_io-t%output::convert_variables')
267  if (io%guard_zones) then
268  l = patch%mesh%lb
269  u = max(patch%mesh%ub,patch%mesh%gn)
270  else
271  l = patch%mesh%li
272  u = patch%mesh%ui
273  end if
274  nt1 = patch%nt-1
275  if (io%format >= 8 .and. iv==patch%idx%d) then
276  ! --- Convert density to log density ---
277  do k=l(3),u(3); do j=l(2),u(2); do i=l(1),u(1)
278  mem(i-l(1)+1,j-l(2)+1,k-l(3)+1,1:nt1) &
279  = log(patch%mem(i,j,k,iv,patch%iit(1:nt1),1))
280  end do; end do; end do
281  else if ((io%format == 6 .or. io%format==8) .and. patch%solver_is('ramses') &
282  .and. iv >= patch%idx%px .and. iv <= patch%idx%pz) then
283  ! --- Convert momentum to velocity ---
284  do k=l(3),u(3); do j=l(2),u(2); do i=l(1),u(1)
285  mem(i-l(1)+1,j-l(2)+1,k-l(3)+1,1:nt1) &
286  = patch%mem(i,j,k,iv,patch%iit(1:nt1),1) &
287  / patch%mem(i,j,k,patch%idx%d,patch%iit(1:nt1),1)
288  end do; end do; end do
289  else
290  mem = patch%mem(l(1):u(1),l(2):u(2),l(3):u(3),iv,patch%iit(1:nt1),1)
291  end if
292  call trace%end()
293  end subroutine
294  !-----------------------------------------------------------------------------
295  subroutine io_write (slot, ider)
296  integer:: slot, ider
297  !---------------------------------------------------------------------------
298  ! This ordering has the derivative order added to the snapshot number, which
299  ! is incremented by (io%time_derivs+1), with patch%id starting with 1, and
300  ! ider and slot starting with 0. This make the time derivative results
301  ! appear to be extra snapshots (thus making accessing the snapshot metadata
302  ! difficult). Introduced for compatibility with Paolo's analysis scripts.
303  !---------------------------------------------------------------------------
304  if (io%method == 'pan') then
305  call mpi_io%iwrite (buffer, 1+ider+slot, patch%id)
306  !---------------------------------------------------------------------------
307  ! This ordering has the derivative order added to the patch%id, which is
308  ! in turn causing an increment by (io%time_derivs+1), so effectively adding
309  ! nv additional variables for each derivative order, making it easy to
310  ! address the derivatives, while using snapshot numbers incremented by 1.
311  !---------------------------------------------------------------------------
312  else
313  call mpi_io%iwrite (buffer, 1+slot, 1+ider+(patch%id-1)*(io%time_derivs+1))
314  end if
315  end subroutine
316 END SUBROUTINE output
317 
318 !===============================================================================
319 !> Write output as one single buffer
320 !===============================================================================
321 SUBROUTINE output_single (patch, count)
322  class(patch_t):: patch
323  integer:: count
324  !.............................................................................
325  integer(8):: snapshot_words, buffer_words, pos
326  integer:: iv, n(3), jt(2), iw, nw, nder, irec, ider, ip, i, iout
327  real(kind=KindScalarVar), dimension(:,:,:,:), pointer:: mem
328  real(kind=4), dimension(:,:,:), pointer:: buffer
329  real:: pt(2)
330  real(8):: start
331  logical:: first_time=.true.
332  character(len=120):: filename
333  integer, save:: itimer=0
334  !-----------------------------------------------------------------------------
335  call trace%begin ('parallel_io_t%output_single', itimer=itimer)
336  !-----------------------------------------------------------------------------
337  ! method='snapshot' => One file per snapshot
338  !-----------------------------------------------------------------------------
339  if (verbose > 1) then
340  write (io_unit%log,*) wallclock(), patch%id, patch%time, ' begin'
341  flush (io_unit%log)
342  end if
343  if (io%method == 'snapshot') then
344  if (.not. associated(files_out)) then
345  io0 = patch%iout
346  io1 = io0 + nint((io%end_time-patch%out_next)/io%out_time) + 1
347  closed = io0-1
348  if (verbose > 0 .and. omp%master) then
349  write (io_unit%log,'(a,i4," -- ",i4)') &
350  'parallel_io_t%output: initializing for snapshots', io0, io1
351  end if
352  allocate (files_out(io0:io1))
353  do i=io0,io0+1
354  write (filename,'(a,"/",i5.5,"/snapshot.dat")') trim(io%outputname), i
355  call files_out(i)%openw (filename)
356  end do
357  end if
358  !---------------------------------------------------------------------------
359  ! If new snapshot number, open a new file
360  !---------------------------------------------------------------------------
361  if (patch%iout+1 > opened) then
362  !$omp atomic write
363  opened = patch%iout+1
364  write (filename,'(a,"/",i5.5,"/snapshot.dat")') trim(io%outputname), opened
365  if (verbose > 1) then
366  write (io_unit%log,*) wallclock(), patch%id, patch%time, ' request_open'
367  flush (io_unit%log)
368  end if
369  call files_out(opened)%request_open (filename)
370  !$omp atomic write
371  check_files = .true.
372  end if
373  end if
374  !-----------------------------------------------------------------------------
375  ! Set size parameters; a buffer is all variables in the patch, and the
376  ! snspshot offset is that times the number of patches per snapshot.
377  !-----------------------------------------------------------------------------
378  if (verbose > 1) then
379  write (io_unit%log,*) wallclock(), patch%id, patch%time, ' set size'
380  flush (io_unit%log)
381  end if
382  if (io%guard_zones) then
383  n = patch%gn
384  else
385  n = patch%n
386  end if
387  !-----------------------------------------------------------------------------
388  ! Slots in the file are ordered (id,iv,iout).
389  !-----------------------------------------------------------------------------
390  if (verbose > 1) then
391  write (io_unit%log,*) wallclock(), patch%id, patch%time, ' call mpi_io'
392  flush (io_unit%log)
393  end if
394  if (io%ntask==0) &
395  call io%abort ('parallel_io_t%output_single: io%ntask == 0')
396  !-----------------------------------------------------------------------------
397  nder = io%time_derivs+1
398  snapshot_words = product(n)*io%ntotal*io%nv*nder
399  buffer_words = product(n)*io%ntask
400  if (io%method == 'snapshot') then
401  call mpi_io%use (files_out(patch%iout))
402  else
403  call mpi_io%use (file_out)
404  end if
405  !-----------------------------------------------------------------------------
406  ! Info
407  !-----------------------------------------------------------------------------
408  if (first_time) then
409  first_time = .false.
410  write(io%output,'(a,f8.3,a,i10,a)') &
411  ' output snapshot size =', snapshot_words*4d0/1024d0**3, ' GB, in', &
412  io%ntotal, ' patches'
413  write(io%output,'(a,4i6,2i12)') &
414  ' parallel_io_t%output: iout, n, snapshot_words, buffer_words =', &
415  patch%iout, n, snapshot_words, buffer_words
416  end if
417  !-----------------------------------------------------------------------------
418  call mpi_io%set (snapshot_words, buffer_words, 99)
419  !-----------------------------------------------------------------------------
420  if (.not.associated(parallel_io%buffer)) then
421  allocate (parallel_io%buffer(n(1),n(2),n(3),io%ntask,io%nv,nder))
422  if (verbose > 0) &
423  write (io%output,*) 'shape io%buffer =', shape(parallel_io%buffer)
424  end if
425  ! --- rank-local patch 1D index; cf. cartesian_mod ---
426  ip = patch%ip
427  if (ip < 0 .or. ip > io%ntask .and. patch%time < io%end_time) then
428  write (io_unit%output,*) &
429  'parallel_io_t%output_single: WARNING io, ntask =', ip, io%ntask
430  ip = max(1,min(ip,io%ntask))
431  end if
432  if (verbose > 0) then
433  write(io%output,'(a,6i6,2i12,l3)') &
434  ' parallel_io_t%output_single: id, ip, iout, n, snapshot_words, buffer_words =', &
435  patch%id, ip, patch%iout, n, snapshot_words, buffer_words, associated(parallel_io%buffer)
436  flush (io%output)
437  end if
438  if (verbose > 1) then
439  write (io_unit%log,*) wallclock(), patch%id, patch%time, ' allocate mem'
440  flush (io_unit%log)
441  end if
442  allocate (mem(n(1),n(2),n(3),patch%nt-1))
443  !-----------------------------------------------------------------------------
444  ! Time interpolation to exact time = out_next.
445  !-----------------------------------------------------------------------------
446  if (verbose > 1) then
447  write (io_unit%log,*) wallclock(), patch%id, patch%time, ' convert var, interpolate'
448  flush (io_unit%log)
449  end if
450  do iv=1,io%nv
451  call convert_variables ! copy `patch%mem` into `mem`, converting if desired
452  buffer => parallel_io%buffer(:,:,:,ip,iv,1)
453  call time_slices%interpolate (patch, mem, buffer)
454  end do
455  !-----------------------------------------------------------------------------
456  ! On-the-fly time derivatives output
457  !-----------------------------------------------------------------------------
458  if (io%time_derivs>0) then
459  do iv=1,io%nv
460  call convert_variables ! copy `patch%mem` into `mem`, converting if desired
461  buffer => parallel_io%buffer(:,:,:,ip,iv,2)
462  call time_slices%derivative1 (patch, mem, buffer)
463  end do
464  if (io%time_derivs==2) then
465  do iv=1,io%nv
466  call convert_variables ! copy `patch%mem` into `mem`, converting if desired
467  buffer => parallel_io%buffer(:,:,:,ip,iv,3)
468  call time_slices%derivative2 (patch, mem, buffer)
469  end do
470  end if
471  end if
472  deallocate (mem)
473  if (verbose > 0) then
474  write (io_unit%log,*) wallclock(), patch%id, patch%time, count, ' buffer'
475  flush (io_unit%log)
476  end if
477  if (count == 0) then
478  if (io%method == 'snapshot') then
479  iout = 0
480  file_out => files_out(patch%iout)
481  else
482  iout = patch%iout
483  end if
484  start = wallclock()
485  if (verbose > 0) then
486  write (io_unit%log,*) wallclock(), patch%id, patch%time, ' start:'//trim(file_out%filename)
487  flush (io_unit%log)
488  end if
489  do ider=1,nder
490  do iv=1,io%nv
491  !-------------------------------------------------------------------------
492  ! Each "record" is product(n)*io%ntasks words long, and contains values
493  ! for one variable, for all the patches in a rank. The ranks are written
494  ! out one after the other for one variable, followed by the next var
495  !-------------------------------------------------------------------------
496  irec = 1 + mpi%rank + mpi%size*(iv-1 + (ider-1)*io%nv)
497  if (verbose > 1) &
498  write (io_unit%output,*) 'write(1): irec, iout =', irec, iout
499  write (io_unit%log,*) 'mpi_io%write: iout, irec =', 1+iout, irec
500  call mpi_io%write (parallel_io%buffer(:,:,:,:,iv,ider), 1+iout, irec)
501  end do
502  end do
503  if (verbose > 0) then
504  write (io_unit%log,*) wallclock(), patch%id, patch%time, ' end:'//trim(file_out%filename)
505  flush (io_unit%log)
506  end if
507  write (io_unit%log,'(a,f7.3,1x,a)') &
508  'parallel_io_t%output_single: MPI_File_write_at I/O', wallclock()-start, 'sec'
509  if (io%method == 'snapshot') then
510  if (verbose > 0) &
511  write (io_unit%log,*) wallclock(), patch%id, patch%time, ' request_close'
512  if (patch%iout-1 > closed) then
513  !$omp atomic write
514  closed = patch%iout-1
515  call files_out(closed)%request_close ()
516  !$omp atomic write
517  check_files = .true.
518  end if
519  end if
520  end if
521  if (verbose > 1) then
522  write (io_unit%log,*) wallclock(), patch%id, patch%time, ' END'
523  flush (io_unit%log)
524  end if
525  call trace%end (itimer)
526 contains
527  !-----------------------------------------------------------------------------
528  ! Copy variables from `patch%mem` to `mem` array, converting variables based
529  ! on the I/O format parameter. This procedure is internal to the 'output_single'
530  ! procedure, and as such only supports io%format in the range 10-13. In the
531  ! process, the time slot ordering is rearranged into increasing time.
532  !-----------------------------------------------------------------------------
533  subroutine convert_variables
534  integer:: nt1, l(3), u(3), i, j, k
535  !...........................................................................
536  call trace%begin ('parallel_io-t%output_single::convert_variables')
537  if (io%guard_zones) then
538  l = patch%mesh%lb
539  u = max(patch%mesh%ub,patch%mesh%gn)
540  else
541  l = patch%mesh%li
542  u = patch%mesh%ui
543  end if
544  nt1 = patch%nt-1
545  if (io%format >= 12 .and. iv==patch%idx%d) then
546  ! --- Convert density to log density ---
547  do k=l(3),u(3); do j=l(2),u(2); do i=l(1),u(1)
548  mem(i-l(1)+1,j-l(2)+1,k-l(3)+1,1:nt1) &
549  = log(patch%mem(i,j,k,iv,patch%iit(1:nt1),1))
550  end do; end do; end do
551  else if ((io%format == 10 .or. io%format==12) .and. patch%solver_is('ramses') &
552  .and. iv >= patch%idx%px .and. iv <= patch%idx%pz) then
553  ! --- Convert momentum to velocity ---
554  do k=l(3),u(3); do j=l(2),u(2); do i=l(1),u(1)
555  mem(i-l(1)+1,j-l(2)+1,k-l(3)+1,1:nt1) &
556  = patch%mem(i,j,k,iv,patch%iit(1:nt1),1) &
557  / patch%mem(i,j,k,patch%idx%d,patch%iit(1:nt1),1)
558  end do; end do; end do
559  else
560  mem = patch%mem(l(1):u(1),l(2):u(2),l(3):u(3),iv,patch%iit(1:nt1),1)
561  end if
562  call trace%end()
563  end subroutine
564 END SUBROUTINE output_single
565 
566 !===============================================================================
567 !> Read input in increments of one variable data cube
568 !===============================================================================
569 SUBROUTINE input (patch, ok)
570  class(patch_t):: patch
571  logical:: ok
572  !.............................................................................
573  integer(8):: snapshot_words, buffer_words
574  integer:: iv, n(3), slot
575  real, dimension(:,:,:,:), pointer:: buffer
576  character(len=120):: filename
577  logical:: single_patch
578  logical, save:: first_time=.true.
579  integer, save:: itimer=0
580  !-----------------------------------------------------------------------------
581  call trace%begin ('parallel_io_t%input', itimer=itimer)
582  !-----------------------------------------------------------------------------
583  ! If the patch is a virtual patch, signal that it should be filled with IC
584  ! data at first -- it will become updated via MPI later. By setting time < 0
585  ! we allow restarting from the IC snapshot, for tests
586  !-----------------------------------------------------------------------------
587  if (patch%is_set(bits%virtual)) then
588  ok = .false.
589  patch%time = -1d0
590  call trace%end (itimer)
591  return
592  end if
593  !-----------------------------------------------------------------------------
594  ! Open the appropriate input file
595  !-----------------------------------------------------------------------------
596  if (io%method == 'snapshot') then
597  write (filename,'(a,"/",i5.5,"/snapshot.dat")') &
598  trim(io%outputname), io%restart
599  else
600  filename=trim(io%inputdir)//'snapshots.dat'
601  end if
602  allocate (file_in)
603  call file_in%openr (filename)
604  !-----------------------------------------------------------------------------
605  ! Choose input reader
606  !-----------------------------------------------------------------------------
607  if (io%format > 9 .and. io%format < 14) then
608  call input_single (patch, ok)
609  call trace%end (itimer)
610  return
611  end if
612  !-----------------------------------------------------------------------------
613  ! Set size parameters
614  !-----------------------------------------------------------------------------
615  if (file_in%handle==0) then
616  call init (patch)
617  end if
618  if (io%guard_zones) then
619  n = patch%gn
620  else
621  n = patch%n
622  end if
623  !-----------------------------------------------------------------------------
624  ! Catch reading from a snapshot with different mpi_size
625  !-----------------------------------------------------------------------------
626  single_patch = any(io%mpi_odims/=io%mpi_dims)
627  !-----------------------------------------------------------------------------
628  ! Compute offset size constants for the mpi_io reader
629  !-----------------------------------------------------------------------------
630  snapshot_words = product(n)*io%ntotal*io%nv*(io%time_derivs+1)
631  if (single_patch) then
632  buffer_words = product(n)
633  else
634  buffer_words = product(n)*io%ntask
635  end if
636  !-----------------------------------------------------------------------------
637  ! Info
638  !-----------------------------------------------------------------------------
639  !$omp critical (first_cr)
640  if (first_time) then
641  first_time = .false.
642  write(io%output,'(a,f8.3,a,i10,a)') &
643  ' output snapshot size =', snapshot_words*4d0/1024d0**3, ' GB, in', &
644  io%ntotal, ' patches'
645  write(io%output,'(a,4i6,2i12)') &
646  ' parallel_io_t%output: iout, n, snapshot_words, buffer_words =', &
647  patch%iout, n, snapshot_words, buffer_words
648  end if
649  !$omp end critical (first_cr)
650  !-----------------------------------------------------------------------------
651  call mpi_io%use (file_in)
652  call mpi_io%set (snapshot_words, buffer_words, io%nwrite*(io%time_derivs+1))
653  !-----------------------------------------------------------------------------
654  ! Read via buffer, in case guard zones not included. FIXME: read snapshot
655  ! guard zone parameter, to allow it to differ
656  !-----------------------------------------------------------------------------
657  allocate (buffer(n(1),n(2),n(3),io%nv))
658  patch%iout = io%restart
659  slot = (io%time_derivs+1)*patch%iout
660  call mpi_io%read (buffer, 1+slot, 1+(patch%id-1)*(io%time_derivs+1))
661  if (mpi_io%err /= 0) then
662  ok = .false.
663  call trace%end (itimer)
664  return
665  end if
666  if (io%guard_zones) then
667  patch%mem(:,:,:,1:io%nv,patch%it,1) = buffer
668  else
669  patch%mem(patch%li(1):patch%ui(1), &
670  patch%li(2):patch%ui(2), &
671  patch%li(3):patch%ui(3),1:io%nv,patch%it,1) = buffer
672  end if
673  if (verbose > 1) then
674  do iv=1,io%nv
675  write (io_unit%output,*) 'parallel_io_t%input: id, iv, iout, min, max =', &
676  patch%id, iv, patch%iout, patch%fminval(iv), patch%fmaxval(iv)
677  end do
678  else if (verbose > 0) then
679  write (io_unit%output,*) 'parallel_io_t%input: id, iout =', &
680  patch%id, patch%iout
681  end if
682  deallocate (buffer)
683  do iv=1,io%nv
684  call convert_variables ! copy `patch%mem` into `mem`, converting if desired
685  end do
686  !-----------------------------------------------------------------------------
687  ! Make sure all time slots have the correct time set
688  !-----------------------------------------------------------------------------
689  patch%t(patch%iit) = patch%time
690  !-----------------------------------------------------------------------------
691  ! Signal success
692  !-----------------------------------------------------------------------------
693  ok = .true.
694  if (verbose==1) &
695  write(io%output,'(a,2i6,1p,2e16.6)') ' parallel_io_t%input: iout, id, time, dtime =', &
696  patch%iout, patch%id, patch%time, patch%dtime
697  if (verbose>1) &
698  write(io%output,'(a,4i6,2i7)') &
699  ' parallel_io_t%input: iout, n, snapshot_words, buffer_words =', &
700  patch%iout, n, snapshot_words, buffer_words
701  call trace%end (itimer)
702 contains
703  !-----------------------------------------------------------------------------
704  ! Copy variables from `patch%mem` to `mem` array, converting variables based
705  ! on the I/O format parameter.
706  !-----------------------------------------------------------------------------
707  subroutine convert_variables
708  !-----------------------------------------------------------------------------
709  ! Using log(d) if io%format==8 or 9
710  !-----------------------------------------------------------------------------
711  if (io%format >= 8 .and. iv==patch%idx%d) then
712  ! --- Convert log density to density ---
713  patch%mem(:,:,:,iv,:,1) = exp(patch%mem(:,:,:,iv,:,1))
714  end if
715  !-----------------------------------------------------------------------------
716  ! Restore to momentum variables for RAMSES snapshots, if format=6,8
717  !-----------------------------------------------------------------------------
718  if ((io%format == 6 .or. io%format==8) .and. patch%solver_is('ramses') &
719  .and. iv >= patch%idx%px .and. iv <= patch%idx%pz) then
720  ! --- Convert velocity to momentum ---
721  patch%mem(:,:,:,iv,:,1) = patch%mem(:,:,:,iv,:,1)*patch%mem(:,:,:,patch%idx%d,:,1)
722  end if
723  end subroutine
724 END SUBROUTINE input
725 
726 !===============================================================================
727 !> Read input in just one read per variable, reading all patches for this rank
728 !> in one go. Use the same parallel_io%buffer as for writes.
729 !===============================================================================
730 SUBROUTINE input_single (patch, ok)
731  class(patch_t):: patch
732  logical:: ok
733  !.............................................................................
734  integer(8):: snapshot_words, buffer_words
735  integer:: iv, n(3), slot, nder, irec, ip, iout, rank, nrank
736  logical, save:: read_failed=.false., first_time=.true.
737  integer, save:: itimer=0
738  real, allocatable:: buffer(:,:,:)
739  logical:: single_patch
740  !-----------------------------------------------------------------------------
741  call trace%begin ('parallel_io_t%input_single', itimer=itimer)
742  !-----------------------------------------------------------------------------
743  ! If the patch is a virtual patch, read nothing, but return ok=.true. never-
744  ! theless, relying on the negative time to force nbor tasks to wait for the
745  ! first MPI refresh of the task.
746  !-----------------------------------------------------------------------------
747  if (patch%is_set(bits%virtual)) then
748  ok = .true.
749  patch%time = -1.0
750  call trace%end (itimer)
751  return
752  end if
753  !-----------------------------------------------------------------------------
754  ! Set size parameters
755  !-----------------------------------------------------------------------------
756  if (file_in%handle==0) then
757  call init (patch)
758  end if
759  if (io%guard_zones) then
760  n = patch%gn
761  else
762  n = patch%n
763  end if
764  !-----------------------------------------------------------------------------
765  ! Catch reading from a snapshot with different mpi_size
766  !-----------------------------------------------------------------------------
767  single_patch = any(io%mpi_odims/=io%mpi_dims)
768  !-----------------------------------------------------------------------------
769  ! Slots in the file are ordered (ip,rank,iv,iout), so all patches for one
770  ! variable that belong to this rank may be read at once
771  !-----------------------------------------------------------------------------
772  if (io%ntask==0) &
773  call io%abort ('parallel_io_t%input_single: io%ntask == 0')
774  !-----------------------------------------------------------------------------
775  nder = io%time_derivs+1
776  snapshot_words = product(n)*io%ntotal*io%nv*nder
777  if (single_patch) then
778  buffer_words = product(n)
779  else
780  buffer_words = product(n)*io%ntask
781  end if
782  !-----------------------------------------------------------------------------
783  !$omp critical (allocate_cr)
784  if (.not.associated(parallel_io%buffer)) then
785  allocate (parallel_io%buffer(n(1),n(2),n(3),io%ntask,io%nv,nder))
786  if (verbose > 0) &
787  write (io%output,*) 'shape io%buffer =', shape(parallel_io%buffer)
788  end if
789  !$omp end critical (allocate_cr)
790  !-----------------------------------------------------------------------------
791  call mpi_io%use (file_in)
792  call mpi_io%set (snapshot_words, buffer_words, 99)
793  !-----------------------------------------------------------------------------
794  ! Make sure only the first thread to come here allocates the buffer, if needed
795  !-----------------------------------------------------------------------------
796  patch%iout = io%restart
797  ip = patch%ip
798  if (single_patch) then
799  !-----------------------------------------------------------------------------
800  ! Info
801  !-----------------------------------------------------------------------------
802  if (first_time) then
803  !$omp critical (first_cr)
804  if (first_time) then
805  write(io%output,'(a,f8.3,a,i10,a)') &
806  ' output snapshot size =', snapshot_words*4d0/1024d0**3, ' GB, in', &
807  io%ntotal, ' patches'
808  write(io%output,'(a,4i6,2i12)') &
809  ' parallel_io_t%output: iout, n, snapshot_words, buffer_words =', &
810  patch%iout, n, snapshot_words, buffer_words
811  first_time = .false.
812  end if
813  !$omp end critical (first_cr)
814  end if
815  allocate (buffer(n(1),n(2),n(3)))
816  !-----------------------------------------------------------------------------
817  ! The previous number of ranks is equal to the product of the previous mpi_dims.
818  ! The previous rank is computable from what the old %ipos would have been.
819  !-----------------------------------------------------------------------------
820  nrank = product(io%mpi_odims)
821  rank = mpi_coords%coords_to_rank (patch%ipos*io%mpi_odims/io%dims)
822  if (verbose > 0) then
823  irec = 1 + rank
824  write (io_unit%output,'(a,i6,4(2x,3i5))') &
825  'single_patch: id, ip, rank, nrank, per_rank, ipos =', &
826  patch%id, ip, rank, nrank, io%dims/io%mpi_odims, patch%ipos, &
827  irec, patch%iout
828  end if
829  do iv=1,io%nv
830  irec = 1 + rank + nrank*(iv-1)
831  if (verbose > 1) &
832  write (io_unit%output,*) 'read(1): irec, iout =', irec, patch%iout
833  call mpi_io%read (buffer, 1+patch%iout, irec)
834  if (io%guard_zones) then
835  patch%mem(:,:,:,iv,patch%it,1) = buffer
836  else
837  patch%mem(patch%li(1):patch%ui(1), &
838  patch%li(2):patch%ui(2), &
839  patch%li(3):patch%ui(3),iv,patch%it,1) = buffer
840  end if
841  end do
842  deallocate (buffer)
843  else
844  if (first_time) then
845  !$omp critical (read_first_cr)
846  if (first_time) then
847  !---------------------------------------------------------------------------
848  ! Info
849  !---------------------------------------------------------------------------
850  write(io%output,'(a,f8.3,a,i10,a)') &
851  ' output snapshot size =', snapshot_words*4d0/1024d0**3, ' GB, in', &
852  io%ntotal, ' patches'
853  write(io%output,'(a,4i6,2i12)') &
854  ' parallel_io_t%output: iout, n, snapshot_words, buffer_words =', &
855  patch%iout, n, snapshot_words, buffer_words
856  !---------------------------------------------------------------------------
857  ! Read via buffer, in case guard zones not included. FIXME: read snapshot
858  ! guard zone parameter, to allow it to differ
859  !---------------------------------------------------------------------------
860  if (io%method == 'snapshot') then
861  iout = 0
862  else
863  iout = patch%iout
864  end if
865  do iv=1,io%nv
866  irec = 1 + mpi%rank + mpi%size*(iv-1)
867  if (verbose > 1) &
868  write (stdout,*) 'read(2): irec, iout =', irec, iout
869  call mpi_io%read (parallel_io%buffer(:,:,:,:,iv,1), 1+iout, irec)
870  if (mpi_io%err /= 0) then
871  read_failed = .true.
872  write(io%output,*) 'io%input_single: read failed', io%restart
873  exit
874  end if
875  end do
876  if (.not.read_failed) &
877  write(io%output,*) 'io%input_single: read succeeded, snapshot', io%restart
878  first_time = .false.
879  end if
880  !$omp end critical (read_first_cr)
881  end if
882  if (read_failed) then
883  ok = .false.
884  call trace%end (itimer)
885  return
886  end if
887  if (ip < 1 .or. ip > io%ntask) then
888  write (io_unit%output,*) &
889  'parallel_io_t%input_single: WARNING io, ntask =', ip, io%ntask
890  flush (io_unit%output)
891  ip = max(1,min(ip,io%ntask))
892  end if
893  if (io%guard_zones) then
894  patch%mem(:,:,:,:,patch%it,1) = parallel_io%buffer(:,:,:,ip,:,1)
895  else
896  patch%mem(patch%li(1):patch%ui(1), &
897  patch%li(2):patch%ui(2), &
898  patch%li(3):patch%ui(3),:,patch%it,1) = parallel_io%buffer(:,:,:,ip,:,1)
899  end if
900  end if
901  if (verbose > 1) then
902  do iv=1,io%nv
903  write (io_unit%output,'(a,i6,3i4,1p,2g14.6)') &
904  'parallel_io_t%input: id, ip, iv, iout, min, max =', &
905  patch%id, ip, iv, patch%iout, patch%fminval(iv), patch%fmaxval(iv)
906  end do
907  else if (verbose > 0) then
908  write (io_unit%output,*) 'parallel_io_t%input: id, ip, iout =', &
909  patch%id, ip, patch%iout
910  end if
911  do iv=1,io%nv
912  call convert_variables ! copy `patch%mem` into `mem`, converting if desired
913  end do
914  !-----------------------------------------------------------------------------
915  ! Make sure all time slots have the correct time set
916  !-----------------------------------------------------------------------------
917  patch%t(patch%iit) = patch%time
918  !-----------------------------------------------------------------------------
919  ! Signal success
920  !-----------------------------------------------------------------------------
921  ok = .true.
922  call trace%end (itimer)
923 contains
924  !-----------------------------------------------------------------------------
925  ! Copy variables from `patch%mem` to `mem` array, converting variables based
926  ! on the I/O format parameter.
927  !-----------------------------------------------------------------------------
928  subroutine convert_variables
929  integer:: l(3), u(3)
930  call trace%begin('input_single::convert_variables')
931  if (io%guard_zones) then
932  l = patch%mesh%lb
933  u = max(patch%mesh%ub,patch%mesh%gn)
934  else
935  l = patch%mesh%li
936  u = patch%mesh%ui
937  end if
938  !-----------------------------------------------------------------------------
939  ! Using log(d) if io%format==12 or 13
940  !-----------------------------------------------------------------------------
941  if (io%format >= 12 .and. iv==patch%idx%d) then
942  ! --- Convert log density to density ---
943  patch%mem(l(1):u(1),l(2):u(2),l(3):u(3),iv,:,1) = &
944  exp(patch%mem(l(1):u(1),l(2):u(2),l(3):u(3),iv,:,1))
945  !-----------------------------------------------------------------------------
946  ! Restore to momentum variables for RAMSES snapshots, if format=10,12
947  !-----------------------------------------------------------------------------
948  else if ((io%format == 10 .or. io%format==12) .and. patch%solver_is('ramses') &
949  .and. iv >= patch%idx%px .and. iv <= patch%idx%pz) then
950  ! --- Convert velocity to momentum ---
951  patch%mem(l(1):u(1),l(2):u(2),l(3):u(3),iv,:,1) = &
952  patch%mem(l(1):u(1),l(2):u(2),l(3):u(3),iv,:,1) * &
953  patch%mem(l(1):u(1),l(2):u(2),l(3):u(3),patch%idx%d,:,1)
954  end if
955  call trace%end()
956  end subroutine
957 END SUBROUTINE input_single
958 
959 END MODULE parallel_io_mod
Compute time derivative from a sequence of time slices, using Lagrange interpolation.
Module for handling blocking and non-blocking MPI parallel I/O to a single file. The module is initia...
Definition: mpi_io_mod.f90:31
Support tic/toc timing, as in MATLAB, and accurate wallclock() function. The timing is generally much...
Use MPI parallel I/O to write everything to a single file. No critical regions should be needed here;...
Template module for patches, which adds pointers to memory and mesh, and number of dimensions and var...
Definition: patch_mod.f90:6
Module for handling blocking and non-blocking MPI parallel I/O to a single file.
Definition: mpi_file_mod.f90:4
Definition: io_mod.f90:4