DISPATCH
direct_io_mod.f90
1 !===============================================================================
2 !> $Id: fcd9be11c70ddcd692251c2adb4f0d1686761af1 $
3 !===============================================================================
5  USE io_mod
6  USE mpi_mod
7  USE trace_mod
8  USE patch_mod
9  USE omp_timer_mod
10  implicit none
11  private
12  !-----------------------------------------------------------------------------
13  ! The output object contains procedures for opening, buffering, writing, and
14  ! closing buffered output to a direct access file.
15  !-----------------------------------------------------------------------------
16  type, public:: direct_io_t
17  integer:: id
18  integer:: dims(4)
19  integer:: count=0
20  integer:: nbuf
21  integer:: rec=0
22  real, dimension(:,:,:,:), allocatable:: buffer
23  character(len=128):: filename
24  logical:: locked=.false.
25  real(8):: position(3)=0d0
26  real:: gb=0.0, gbs=0.0
27  contains
28  procedure:: init
29  procedure:: output
30  procedure:: input
31  procedure:: out
32  procedure:: in
33  procedure:: close
34  end type
35  integer(8):: revision=1
36  integer:: verbose=1
37  type(direct_io_t), public:: direct_io
38 CONTAINS
39 
40 !===============================================================================
41 !> The unigrid size is available from the mesh box and patch sizes
42 !===============================================================================
43 SUBROUTINE init (self, patch, filename)
44  class(direct_io_t) :: self
45  class(patch_t):: patch
46  character(len=*), optional :: filename
47  integer:: mpatch(3)
48  !-----------------------------------------------------------------------------
49  call trace_begin('direct_io_t%init')
50  !$omp critical (direct_cr)
51  !
52  if (present(filename)) then
53  if (verbose>1) print *, 'direct_io_t: filename = ', trim(filename)
54  self%filename = filename
55  end if
56  !----------------------------------------------------------------------------
57  ! Cartesian dimensions of patch arrangement and unigrid dimensions self%dims
58  !----------------------------------------------------------------------------
59  mpatch = nint(patch%mesh%b/patch%mesh%s)
60  if (patch%no_mans_land) then
61  self%dims(1:3) = mpatch*patch%mesh%n
62  else
63  self%dims(1:3) = mpatch*(patch%mesh%n-1)+1
64  end if
65  self%id = patch%id
66  self%dims(4) = patch%nv
67  self%nbuf = product(mpatch)
68  !
69  self%count = 0
70  !$omp end critical (direct_cr)
71  call trace_end
72 END SUBROUTINE init
73 
74 !===============================================================================
75 !> Output a patch to a direct access unigrid data file
76 !===============================================================================
77 SUBROUTINE output (self, patch)
78  class(direct_io_t):: self
79  class(patch_t):: patch
80  !-----------------------------------------------------------------------------
81  character(len=64) :: filename
82  real, pointer :: buf(:,:,:,:)
83  integer :: n(3), l(3), u(3), ng(3)=0, jt(2), id=0, offset(4)=1
84  real :: pt(2)
85  !-----------------------------------------------------------------------------
86  call trace_begin('patch_t%output_direct')
87  !$omp critical (direct_cr)
88  !-----------------------------------------------------------------------------
89  ! Write a small text file with info when the last chunk is being entered
90  !-----------------------------------------------------------------------------
91  if (self%count==self%nbuf-1) then
92  write (filename,'(a,i5.5,"/",i5.5)') trim(io%outputname), patch%iout, id
93  if (patch%no_mans_land) then
94  io%format = 1
95  else
96  io%format = 2
97  end if
98  open (io%data_unit, file=trim(filename)//'.txt', form='formatted', status='unknown')
99  write (io%data_unit,'(i2.2)') io%format
100  write (io%data_unit,'(a)') trim(patch%kind)
101  write (io%data_unit,'(a)') trim(patch%eos)
102  write (io%data_unit,'(a)') trim(patch%opacity)
103  write (io%data_unit,'(1p,2e18.10,6(2x,3e18.10)1x,1x,2e10.3)') patch%out_next, patch%dtime, &
104  self%position, patch%ds, patch%velocity, patch%llc_nat, patch%llc_cart, &
105  patch%mesh%centre_nat, patch%quality, patch%gamma
106  l = 1
107  u = self%dims(1:3)
108  n = u-l+1
109  id = mpi%rank
110  write (io%data_unit,'(6(3i10.1,2x),3i2.1,1x,1i2.1,i3)') int(patch%box/patch%ds+0.5), &
111  n, l, u, n, n, ng, patch%nv, patch%level
112  write (io%data_unit,'(3(1i8.1,1x),i2.1)') id, patch%iout, patch%istep, patch%mesh_type
113  close (io%data_unit)
114  self%filename = trim(filename)//'.dat'
115  end if
116  !-----------------------------------------------------------------------------
117  ! Pass a chunk to the procedure that eventually writes a binary file
118  !-----------------------------------------------------------------------------
119  l = patch%mesh%li
120  u = patch%mesh%ui
121  n = u-l+1
122  if (verbose>1) print *, 'output:', patch%id, l, u
123  allocate (buf(n(1),n(2),n(3),patch%nv))
124  call patch%time_interval (patch%out_next, jt, pt)
125  buf = patch%mem(l(1):u(1),l(2):u(2),l(3):u(3),:,jt(1),1)*pt(1) &
126  + patch%mem(l(1):u(1),l(2):u(2),l(3):u(3),:,jt(2),1)*pt(2)
127  offset(1:3) = 1 + patch%ipos*patch%ncell
128  call self%out (buf, offset)
129  deallocate (buf)
130  !$omp end critical (direct_cr)
131  call trace_end
132 END SUBROUTINE output
133 
134 !===============================================================================
135 !> Output a patch to a direct access unigrid data file
136 !===============================================================================
137 SUBROUTINE input (self, patch, ok)
138  class(direct_io_t):: self
139  class(patch_t):: patch
140  logical:: ok
141  !-----------------------------------------------------------------------------
142  character(len=64) :: filename
143  real, pointer :: buf(:,:,:,:)
144  integer :: n(3), l(3), u(3), ng(3)=0, jt(2), id=0, offset(4)=1, iof
145  real :: pt(2)
146  logical :: exist
147  real(8), save :: time
148  !-----------------------------------------------------------------------------
149  call trace_begin('patch_t%input_direct')
150  !$omp critical (direct_in_cr)
151  !-----------------------------------------------------------------------------
152  ! Read a small text file with info, to get the snapshot time
153  !-----------------------------------------------------------------------------
154  write (filename,'(a,i5.5,"/",i5.5)') trim(io%outputname), patch%restart, id
155  inquire (file=trim(filename)//'.dat', exist=exist)
156  if (exist) then
157  if (self%count==0) then
158  open (io%data_unit, file=trim(filename)//'.txt', form='formatted', status='old')
159  read (io%data_unit,'(i2.2)') iof
160  read (io%data_unit,'(a)')
161  read (io%data_unit,'(a)')
162  read (io%data_unit,'(a)')
163  read (io%data_unit,'(1p,2e18.10,6(2x,3e18.10)1x,1x,2e10.3)') time
164  close (io%data_unit)
165  self%filename = trim(filename)//'.dat'
166  end if
167  patch%time = time
168  !-----------------------------------------------------------------------------
169  ! Get a chunk from the procedure that reads a whole piece
170  !-----------------------------------------------------------------------------
171  l = patch%mesh%li
172  u = patch%mesh%ui
173  n = u-l+1
174  allocate (buf(n(1),n(2),n(3),patch%nv))
175  !call patch%time_interval (patch%out_next, jt, pt)
176  offset(1:3) = 1 + patch%ipos*patch%ncell
177  call self%in (buf, offset)
178  patch%mem( : , : , : ,:,patch%it,1) = 0.0
179  patch%mem(l(1):u(1),l(2):u(2),l(3):u(3),:,patch%it,1) = buf
180  if (verbose>1) print *, 'id, dmin, dmax', patch%id, minval(buf(:,:,:,1)), maxval(buf(:,:,:,1)), l, u
181  if (verbose>1) print *, 'id, dmin, dmax', patch%id, patch%fminval(patch%idx%d), patch%fmaxval(patch%idx%d), shape(buf)
182  deallocate (buf)
183  ok = .true.
184  if (self%count==self%nbuf) then
185  call self%close
186  end if
187  else
188  ok = .false.
189  end if
190  if (io%verbose > 0) &
191  print *, 'input_direct:', trim(filename)//'.dat', exist, ok, self%count
192  !$omp end critical (direct_in_cr)
193  call trace_end
194 END SUBROUTINE input
195 
196 !===============================================================================
197 !> Interface to optionally enforce a critical region
198 !===============================================================================
199 SUBROUTINE out (self, data, offset)
200  class(direct_io_t):: self
201  real, dimension(:,:,:,:) :: data
202  integer, dimension(4) :: offset
203  !-----------------------------------------------------------------------------
204  call trace_begin('direct_io_t%out')
205  if (self%locked) then
206  call out_real (self, data, offset)
207  else
208  call out_real (self, data, offset)
209  end if
210  call trace_end
211 END SUBROUTINE out
212 
213 !===============================================================================
214 !> Buffer one patch, accumulating until the buffer is complete
215 !===============================================================================
216 SUBROUTINE out_real (self, data, offset)
217  class(direct_io_t):: self
218  integer, dimension(4) :: offset, n, l, u
219  real, dimension(:,:,:,:) :: data
220  !-----------------------------------------------------------------------------
221  call trace_begin('direct_io_t%out_real')
222  n = shape(data)
223  l = offset
224  u = l+n-1
225  if (io%debug(3)) &
226  print '("direct_io_t%out",4(4i4,2x))', l, u, self%dims, shape(data)
227  if (.not.allocated(self%buffer)) then
228  allocate (self%buffer(self%dims(1),self%dims(2),self%dims(3),self%dims(4)))
229  end if
230  self%buffer(l(1):u(1),l(2):u(2),l(3):u(3),l(4):u(4)) = data
231  self%count = self%count+1
232  if (self%count==self%nbuf) call write_out (self)
233  call trace_end
234 END SUBROUTINE out_real
235 
236 !===============================================================================
237 !> Interface to optionally enforce a critical region
238 !===============================================================================
239 SUBROUTINE in (self, data, offset)
240  class(direct_io_t):: self
241  real, dimension(:,:,:,:) :: data
242  integer, dimension(4) :: offset
243  !-----------------------------------------------------------------------------
244  call trace_begin('direct_io_t%in')
245  if (self%locked) then
246  call in_real (self, data, offset)
247  else
248  call in_real (self, data, offset)
249  end if
250  call trace_end
251 END SUBROUTINE in
252 
253 !===============================================================================
254 !> Buffer one patch, accumulating until the buffer is complete
255 !===============================================================================
256 SUBROUTINE in_real (self, data, offset)
257  class(direct_io_t):: self
258  integer, dimension(4) :: offset, n, l, u
259  real, dimension(:,:,:,:) :: data
260  !-----------------------------------------------------------------------------
261  call trace_begin('direct_io_t%in_real')
262  n = shape(data)
263  l = offset
264  u = l+n-1
265  if (io%debug(3)) &
266  print '("direct_io_t%in",4(4i4,2x))', l, u, self%dims, shape(data)
267  if (.not.allocated(self%buffer)) then
268  allocate (self%buffer(self%dims(1),self%dims(2),self%dims(3),self%dims(4)))
269  end if
270  if (self%count==0) call read_in (self)
271  data = self%buffer(l(1):u(1),l(2):u(2),l(3):u(3),l(4):u(4))
272  self%count = self%count+1
273  call trace_end
274 END SUBROUTINE in_real
275 
276 !===============================================================================
277 !> Close an output = deallocate buffer
278 !===============================================================================
279 SUBROUTINE close (self)
280  class(direct_io_t):: self
281  !-----------------------------------------------------------------------------
282  call trace_begin('direct_io_t%close')
283  if (allocated(self%buffer)) then
284  if (verbose>1) print *, 'deallocating buffer'
285  deallocate (self%buffer)
286  end if
287  call trace_end
288 END SUBROUTINE close
289 
290 !===============================================================================
291 !> Write out the buffer and reset the counter
292 !===============================================================================
293 SUBROUTINE write_out (self, rec)
294  class(direct_io_t) :: self
295  integer, optional :: rec
296  integer :: l_rec, recl, i, j, i_rec
297  integer(8) :: recl8
298  real(8) :: wt
299  !-----------------------------------------------------------------------------
300  call trace_begin('direct_io_t%write_out')
301  self%locked = .true.
302  if (present(rec)) then
303  l_rec = rec
304  else
305  l_rec = 1
306  end if
307  self%rec = l_rec
308  wt = wallclock()
309  if (sizeof(self%buffer) < 2e9) then
310  recl = 4*product(self%dims)
311  open (io_unit%direct, file=trim(self%filename), access='direct', &
312  status='unknown', recl=recl)
313  write (io_unit%direct, rec=l_rec) self%buffer
314  else
315  recl8 = 4_8*product(self%dims(1:3))
316  if (recl8 < 2e9) then
317  recl = recl8
318  open (io_unit%direct, file=trim(self%filename), access='direct', &
319  status='unknown', recl=recl)
320  do i=1,self%dims(4)
321  write (io_unit%direct, rec=i+(l_rec-1)*self%dims(4)) self%buffer(:,:,:,i)
322  end do
323  else
324  if (verbose>1) print *, 'DOUBLE LOOP'
325  recl = 4*product(self%dims(1:2))
326  open (io_unit%direct, file=trim(self%filename), access='direct', &
327  status='unknown', recl=recl)
328  do j=1,self%dims(4)
329  do i=1,self%dims(3)
330  i_rec = i + (j-1)*self%dims(3) + (l_rec-1)*self%dims(4)
331  write (io_unit%direct, rec=i_rec) self%buffer(:,:,i,j)
332  end do
333  end do
334  end if
335  end if
336  wt = max(wallclock()-wt,1d-9)
337  self%gb = 4.0*product(self%dims)/1024.**3
338  self%gbs = self%gb/wt
339  self%count = 0
340  close (io_unit%direct)
341  self%locked = .false.
342  call trace_end
343 END SUBROUTINE write_out
344 
345 !===============================================================================
346 !> Read in the buffer and reset the counter
347 !===============================================================================
348 SUBROUTINE read_in (self, rec)
349  class(direct_io_t) :: self
350  integer, optional :: rec
351  integer :: l_rec, recl, i, j, i_rec
352  integer(8) :: recl8
353  real(8) :: wt
354  !-----------------------------------------------------------------------------
355  call trace_begin('direct_io_t%read_in')
356  self%locked = .true.
357  if (present(rec)) then
358  l_rec = rec
359  else
360  l_rec = 1
361  end if
362  self%rec = l_rec
363  wt = wallclock()
364  if (sizeof(self%buffer) < 2e9) then
365  recl = 4*product(self%dims)
366  open (io_unit%direct, file=trim(self%filename), access='direct', &
367  status='old', recl=recl)
368  read (io_unit%direct, rec=l_rec) self%buffer
369  else
370  recl8 = 4_8*product(self%dims(1:3))
371  if (recl8 < 2e9) then
372  recl = recl8
373  open (io_unit%direct, file=trim(self%filename), access='direct', &
374  status='old', recl=recl)
375  do i=1,self%dims(4)
376  read (io_unit%direct, rec=i+(l_rec-1)*self%dims(4)) self%buffer(:,:,:,i)
377  end do
378  else
379  if (verbose>1) print *, 'DOUBLE LOOP'
380  recl = 4*product(self%dims(1:2))
381  open (io_unit%direct, file=trim(self%filename), access='direct', &
382  status='old', recl=recl)
383  do j=1,self%dims(4)
384  do i=1,self%dims(3)
385  i_rec = i + (j-1)*self%dims(3) + (l_rec-1)*self%dims(4)
386  read (io_unit%direct, rec=i_rec) self%buffer(:,:,i,j)
387  end do
388  end do
389  end if
390  end if
391  print *, 'read_in:', shape(self%buffer), minval(self%buffer), maxval(self%buffer)
392  wt = max(wallclock()-wt,1d-9)
393  self%gb = 4.0*product(self%dims)/1024.**3
394  self%gbs = self%gb/wt
395  self%count = 0
396  close (io_unit%direct)
397  self%locked = .false.
398  call trace_end
399 END SUBROUTINE read_in
400 
401 END MODULE direct_io_mod
402 
Support tic/toc timing, as in MATLAB, and accurate wallclock() function. The timing is generally much...
Template module for patches, which adds pointers to memory and mesh, and number of dimensions and var...
Definition: patch_mod.f90:6
Definition: io_mod.f90:4