KTH framework for Nek5000 toolboxes; testing version  0.0.1
lpm_io.f
Go to the documentation of this file.
1 !-----------------------------------------------------------------------
2  subroutine lpm_io_read(filein1,npart)
3  include 'SIZE'
4  include 'SOLN'
5  include 'INPUT'
6  include 'MASS'
7  include 'GEOM'
8  include 'CTIMER'
9  include 'TSTEP'
10  include 'PARALLEL'
11 # include "LPM"
12 
13  real*4 rout_pos(3 *lpm_lpart)
14  > ,rout_sln(lpm_lrs*lpm_lpart)
15  > ,rout_lrp(lpm_lrp*lpm_lpart)
16  > ,rout_lip(3 *lpm_lpart)
17 
18  common /nekmpi/ mid,mp,nekcomm,nekgroup,nekreal
19 
20  real tcoef(3,3),dt_cmt,time_cmt
21  common /timestepcoef/ tcoef,dt_cmt,time_cmt
22 
23  character*5 sprop1
24  character*9 rprop1
25 
26  character (len = *) filein1
27  character*3 filein
28  character*12 vtufile
29  character*13 vtufile1
30  character*50 dumstr
31  character*1 dum_read
32 
33  integer icalld1
34  save icalld1
35  data icalld1 /0/
36 
37  logical partl
38  integer vtu,vtu1,prevs(2,np)
39  integer*4 iint
40  integer*8 idisp_pos,idisp_sln,idisp_lrp,idisp_lip,disp
41  integer*8 stride_len
42  real*4 rnptot
43 
44  icalld1 = icalld1+1
45 
46  lpm_restart = .true.
47 
48  nnp = np
49  nxx = lpm_npart
50  npt_total = iglsum(nxx,1)
51 
52  jx = 1
53  jy = 2
54  jz = 3
55 
56 ! --------------------------------------------------
57 ! COPY PARTICLES TO OUTPUT ARRAY
58 ! --------------------------------------------------
59 
60  iadd = 0
61  if_pos = 3 *isize*npt_total
62  if_sln = lpm_lrs*isize*npt_total
63  if_lrp = lpm_lrp*isize*npt_total
64  if_lip = 3 *isize*npt_total
65 
66  if (nid .eq. 0) then
67 
68  vtu=867+nid
69  open(unit=vtu,file=filein1,access='stream',form="unformatted")
70 
71  ivtu_size = -1
72  ifound = 0
73  do i=1,1000000
74  read(vtu) dum_read
75  if (dum_read == '_') ifound = ifound + 1
76  if (ifound .eq. 2) then
77  ivtu_size = i
78  exit
79  endif
80  enddo
81  read(vtu) npt_total
82  close(vtu)
83 
84  npt_total = npt_total/isize/3
85 
86  endif
87 
88  call bcast(npt_total,isize)
89 
90  idisp_pos = ivtu_size + isize*(1)
91  idisp_sln = ivtu_size + isize*(3 *npt_total + 2)
92  idisp_lrp = ivtu_size + isize*(3 *npt_total + lpm_lrs*npt_total
93  > + 3)
94  idisp_lip = ivtu_size + isize*(3 *npt_total + lpm_lrs*npt_total
95  > + lpm_lrp*npt_total + 4)
96 
97  npmax = min(npt_total/lpm_lpart+1,mp)
98  stride_len = 0
99  if (nid .le. npmax-1 .and. nid. ne. 0) stride_len = nid*lpm_lpart
100 
101  npart = lpm_lpart
102  if (nid .gt. npmax-1) npart = 0
103 
104  ndiff = npt_total - (npmax-1)*lpm_lpart
105  if (nid .eq. npmax-1) npart = ndiff
106 
107  idisp_pos = idisp_pos + isize*3 *stride_len
108  idisp_sln = idisp_sln + isize*lpm_lrs*stride_len
109  idisp_lrp = idisp_lrp + isize*lpm_lrp*stride_len
110  idisp_lip = idisp_lip + isize*3 *stride_len
111 
112  icount_pos = npart*3
113  icount_sln = npart*lpm_lrs
114  icount_lrp = npart*lpm_lrp
115  icount_lip = npart*3
116 
117  iorank = -1
118 
119  call byte_open_mpi(filein1,pth,.true.,ierr)
120 
121 
122  call byte_set_view(idisp_pos,pth)
123  call byte_read_mpi(rout_pos,icount_pos,iorank,pth,ierr)
124  call byte_set_view(idisp_sln,pth)
125  call byte_read_mpi(rout_sln,icount_sln,iorank,pth,ierr)
126  call byte_set_view(idisp_lrp,pth)
127  call byte_read_mpi(rout_lrp,icount_lrp,iorank,pth,ierr)
128  call byte_set_view(idisp_lip,pth)
129  call byte_read_mpi(rout_lip,icount_lip,iorank,pth,ierr)
130 
131 
132  call byte_close_mpi(pth,ierr)
133 
134  ic_pos = 0
135  ic_sln = 0
136  ic_lrp = 0
137  ic_lip = 0
138  do i=1,npart
139 
140  ic_pos = ic_pos + 1
141  lpm_y(jx,i) = rout_pos(ic_pos)
142  ic_pos = ic_pos + 1
143  lpm_y(jy,i) = rout_pos(ic_pos)
144  if (if3d) then
145  ic_pos = ic_pos + 1
146  lpm_y(jz,i) = rout_pos(ic_pos)
147  else
148  ic_pos = ic_pos + 1
149  lpm_y(jz,i) = 0.0
150  endif
151 
152  do j=1,lpm_lrs
153  ic_sln = ic_sln + 1
154  lpm_y(j,i) = rout_sln(ic_sln)
155  enddo
156 
157  do j=1,lpm_lrp
158  ic_lrp = ic_lrp + 1
159  lpm_rprop(j,i) = rout_lrp(ic_lrp)
160  enddo
161 
162  ic_lip = ic_lip + 1
163  lpm_iprop(5,i) = int(rout_lip(ic_lip))
164  ic_lip = ic_lip + 1
165  lpm_iprop(6,i) = int(rout_lip(ic_lip))
166  ic_lip = ic_lip + 1
167  lpm_iprop(7,i) = int(rout_lip(ic_lip))
168 
169  enddo
170 
171  return
172  end
173 !-----------------------------------------------------------------------
174  subroutine lpm_io_write(filein1,iobig)
175  include 'SIZE'
176  include 'SOLN'
177  include 'INPUT'
178  include 'MASS'
179  include 'GEOM'
180  include 'CTIMER'
181  include 'TSTEP'
182  include 'PARALLEL'
183 # include "LPM"
184 
185  real*4 rout_pos(3 *lpm_lpart)
186  > ,rout_sln(lpm_lrs*lpm_lpart)
187  > ,rout_lrp(lpm_lrp*lpm_lpart)
188  > ,rout_lip(3 *lpm_lpart)
189 
190  real tcoef(3,3),dt_cmt,time_cmt
191  common /timestepcoef/ tcoef,dt_cmt,time_cmt
192 
193  character*5 sprop1
194  character*9 rprop1
195 
196  character filein1*(*)
197  character*3 filein
198  character*12 vtufile
199  character*13 vtufile1
200  character*50 dumstr
201 
202  integer icalld1
203  save icalld1
204  data icalld1 /0/
205 
206  logical partl
207  integer vtu,vtu1,prevs(2,np)
208  integer*4 iint
209  integer*8 idisp_pos,idisp_sln,idisp_lrp,idisp_lip
210  integer*8 stride_len
211 
212  icalld1 = icalld1+1
213 
214  nnp = np
215  nxx = lpm_npart
216  npt_total = iglsum(nxx,1)
217 
218  jx = 1
219  jy = 2
220  jz = 3
221 
222  if_sz = len(filein1)
223  if (if_sz .lt. 3) then
224  filein = 'par'
225  else
226  write(filein,'(A3)') filein1
227  endif
228 
229 ! --------------------------------------------------
230 ! COPY PARTICLES TO OUTPUT ARRAY
231 ! --------------------------------------------------
232 
233  iadd = 0
234  if_pos = 3 *isize*npt_total
235  if_sln = lpm_lrs*isize*npt_total
236  if_lrp = lpm_lrp*isize*npt_total
237  if_lip = 3 *isize*npt_total
238 
239  ic_pos = iadd
240  ic_sln = iadd
241  ic_lrp = iadd
242  ic_lip = iadd
243  do i=1,nxx
244 
245  ic_pos = ic_pos + 1
246  rout_pos(ic_pos) = lpm_y(jx,i)
247  ic_pos = ic_pos + 1
248  rout_pos(ic_pos) = lpm_y(jy,i)
249  if (if3d) then
250  ic_pos = ic_pos + 1
251  rout_pos(ic_pos) = lpm_y(jz,i)
252  else
253  ic_pos = ic_pos + 1
254  rout_pos(ic_pos) = 0.0
255  endif
256 
257  do j=1,lpm_lrs
258  ic_sln = ic_sln + 1
259  rout_sln(ic_sln) = lpm_y(j,i)
260  enddo
261 
262  do j=1,lpm_lrp
263  ic_lrp = ic_lrp + 1
264  rout_lrp(ic_lrp) = lpm_rprop(j,i)
265  enddo
266 
267  ic_lip = ic_lip + 1
268  rout_lip(ic_lip) = lpm_iprop(5,i)
269  ic_lip = ic_lip + 1
270  rout_lip(ic_lip) = lpm_iprop(6,i)
271  ic_lip = ic_lip + 1
272  rout_lip(ic_lip) = lpm_iprop(7,i)
273 
274  enddo
275 
276 ! --------------------------------------------------
277 ! FIRST GET HOW MANY PARTICLES WERE BEFORE THIS RANK
278 ! --------------------------------------------------
279  do i=1,nnp
280  prevs(1,i) = i-1
281  prevs(2,i) = nxx
282  enddo
283 
284  nps = 1 ! index of new proc for doing stuff
285  nglob = 1 ! unique key to sort by
286  nkey = 1 ! number of keys (just 1 here)
287  ndum = 2
288  call fgslib_crystal_ituple_transfer(i_cr_hndl,prevs,
289  > ndum,nnp,nnp,nps)
290  call fgslib_crystal_ituple_sort(i_cr_hndl,prevs,
291  > ndum,nnp,nglob,nkey)
292 
293  stride_len = 0
294  if (nid .ne. 0) then
295  do i=1,nid
296  stride_len = stride_len + prevs(2,i)
297  enddo
298  endif
299 
300 ! ----------------------------------------------------
301 ! WRITE EACH INDIVIDUAL COMPONENT OF A BINARY VTU FILE
302 ! ----------------------------------------------------
303  write(vtufile,'(A3,I5.5,A4)') filein,icalld1,'.vtu'
304 
305  if (nid .eq. 0) then
306 
307  vtu=867+nid
308  open(unit=vtu,file=vtufile,status='replace')
309 
310 ! ------------
311 ! FRONT MATTER
312 ! ------------
313  write(vtu,'(A)',advance='no') '<VTKFile '
314  write(vtu,'(A)',advance='no') 'type="UnstructuredGrid" '
315  write(vtu,'(A)',advance='no') 'version="1.0" '
316  if (iobig .eq. 0) then
317  write(vtu,'(A)',advance='yes') 'byte_order="LittleEndian">'
318  elseif (iobig .eq. 1) then
319  write(vtu,'(A)',advance='yes') 'byte_order="BigEndian">'
320  endif
321 
322  write(vtu,'(A)',advance='yes') ' <UnstructuredGrid>'
323 
324  write(vtu,'(A)',advance='yes') ' <FieldData>'
325  write(vtu,'(A)',advance='no') ' <DataArray ' ! time
326  write(vtu,'(A)',advance='no') 'type="Float32" '
327  write(vtu,'(A)',advance='no') 'Name="TIME" '
328  write(vtu,'(A)',advance='no') 'NumberOfTuples="1" '
329  write(vtu,'(A)',advance='no') 'format="ascii"> '
330  write(vtu,'(E14.7)',advance='no') time
331  write(vtu,'(A)',advance='yes') ' </DataArray> '
332  write(vtu,'(A)',advance='no') ' <DataArray ' ! cycle
333  write(vtu,'(A)',advance='no') 'type="Int32" '
334  write(vtu,'(A)',advance='no') 'Name="CYCLE" '
335  write(vtu,'(A)',advance='no') 'NumberOfTuples="1" '
336  write(vtu,'(A)',advance='no') 'format="ascii"> '
337  write(vtu,'(I0)',advance='no') istep
338  write(vtu,'(A)',advance='yes') ' </DataArray> '
339 
340  write(vtu,'(A)',advance='yes') ' </FieldData>'
341  write(vtu,'(A)',advance='no') ' <Piece '
342  write(vtu,'(A)',advance='no') 'NumberOfPoints="'
343  write(vtu,'(I0)',advance='no') npt_total
344  write(vtu,'(A)',advance='yes') '" NumberOfCells="0"> '
345 
346 ! -----------
347 ! COORDINATES
348 ! -----------
349  iint = 0
350  write(vtu,'(A)',advance='yes') ' <Points>'
351  call lpm_io_vtu_data(vtu,"Position",3 ,iint)
352  iint = iint + 3 *isize*npt_total + isize
353  write(vtu,'(A)',advance='yes') ' </Points>'
354 
355 ! ----
356 ! DATA
357 ! ----
358  write(vtu,'(A)',advance='yes') ' <PointData>'
359 
360  call lpm_io_vtu_data(vtu,'lpm-y',lpm_lrs,iint)
361  iint = iint + lpm_lrs*isize*npt_total + isize
362 
363  call lpm_io_vtu_data(vtu,'lpm-rprop',lpm_lrp,iint)
364  iint = iint + lpm_lrp*isize*npt_total + isize
365 
366  call lpm_io_vtu_data(vtu,'lpm-iprop',3,iint)
367  iint = iint + 3*isize*npt_total + isize
368 
369  write(vtu,'(A)',advance='yes') ' </PointData> '
370 
371 ! ----------
372 ! END MATTER
373 ! ----------
374  write(vtu,'(A)',advance='yes') ' <Cells> '
375  write(vtu,'(A)',advance='no') ' <DataArray '
376  write(vtu,'(A)',advance='no') 'type="Int32" '
377  write(vtu,'(A)',advance='no') 'Name="connectivity" '
378  write(vtu,'(A)',advance='yes') 'format="ascii"/> '
379  write(vtu,'(A)',advance='no') ' <DataArray '
380  write(vtu,'(A)',advance='no') 'type="Int32" '
381  write(vtu,'(A)',advance='no') 'Name="offsets" '
382  write(vtu,'(A)',advance='yes') 'format="ascii"/> '
383  write(vtu,'(A)',advance='no') ' <DataArray '
384  write(vtu,'(A)',advance='no') 'type="Int32" '
385  write(vtu,'(A)',advance='no') 'Name="types" '
386  write(vtu,'(A)',advance='yes') 'format="ascii"/> '
387  write(vtu,'(A)',advance='yes') ' </Cells> '
388  write(vtu,'(A)',advance='yes') ' </Piece> '
389  write(vtu,'(A)',advance='yes') ' </UnstructuredGrid> '
390 
391 ! -----------
392 ! APPEND DATA
393 ! -----------
394  write(vtu,'(A)',advance='no') ' <AppendedData encoding="raw">'
395  close(vtu)
396 
397  open(unit=vtu,file=vtufile,access='stream',form="unformatted"
398  > ,position='append')
399  write(vtu) '_'
400  close(vtu)
401 
402  inquire(file=vtufile,size=ivtu_size)
403  endif
404 
405  call bcast(ivtu_size, isize)
406 
407  ! byte-displacements
408  idisp_pos = ivtu_size + isize*(3 *stride_len + 1)
409  idisp_sln = ivtu_size + isize*(3 *npt_total + lpm_lrs*stride_len
410  > + 2)
411  idisp_lrp = ivtu_size + isize*(3 *npt_total + lpm_lrs*npt_total
412  > + lpm_lrp*stride_len + 3)
413  idisp_lip = ivtu_size + isize*(3 *npt_total + lpm_lrs*npt_total
414  > + lpm_lrp*npt_total + 3*stride_len + 4 )
415 
416  ! how much to write
417  icount_pos = 3 *nxx
418  icount_sln = lpm_lrs*nxx
419  icount_lrp = lpm_lrp*nxx
420  icount_lip = 3 *nxx
421 
422  iorank = -1
423 
424  ! integer write
425  if (nid .eq. 0) then
426  open(unit=vtu,file=vtufile,access='stream',form="unformatted"
427  > ,position='append')
428  write(vtu) if_pos
429  close(vtu)
430  endif
431 
432  rdum = dnekclock_sync()
433 
434  ! write
435  call byte_open_mpi(vtufile,pth,.false.,ierr)
436  call byte_set_view(idisp_pos,pth)
437  call byte_write_mpi(rout_pos,icount_pos,iorank,pth,ierr)
438  call byte_close_mpi(pth,ierr)
439 
440  rdum = dnekclock_sync()
441 
442  ! integer write
443  if (nid .eq. 0) then
444  open(unit=vtu,file=vtufile,access='stream',form="unformatted"
445  > ,position='append')
446  write(vtu) if_sln
447  close(vtu)
448  endif
449 
450  rdum = dnekclock_sync()
451 
452  ! write
453  call byte_open_mpi(vtufile,pth,.false.,ierr)
454  call byte_set_view(idisp_sln,pth)
455  call byte_write_mpi(rout_sln,icount_sln,iorank,pth,ierr)
456  call byte_close_mpi(pth,ierr)
457 
458  ! integer write
459  if (nid .eq. 0) then
460  open(unit=vtu,file=vtufile,access='stream',form="unformatted"
461  > ,position='append')
462  write(vtu) if_lrp
463  close(vtu)
464  endif
465 
466  rdum = dnekclock_sync()
467 
468  ! write
469  call byte_open_mpi(vtufile,pth,.false.,ierr)
470  call byte_set_view(idisp_lrp,pth)
471  call byte_write_mpi(rout_lrp,icount_lrp,iorank,pth,ierr)
472  call byte_close_mpi(pth,ierr)
473 
474  ! integer write
475  if (nid .eq. 0) then
476  open(unit=vtu,file=vtufile,access='stream',form="unformatted"
477  > ,position='append')
478  write(vtu) if_lip
479  close(vtu)
480  endif
481 
482  ! write
483  call byte_open_mpi(vtufile,pth,.false.,ierr)
484  call byte_set_view(idisp_lip,pth)
485  call byte_write_mpi(rout_lip,icount_lip,iorank,pth,ierr)
486  call byte_close_mpi(pth,ierr)
487 
488  if (nid .eq. 0) then
489  vtu=867+nid
490  open(unit=vtu,file=vtufile,status='old',position='append')
491 
492  write(vtu,'(A)',advance='yes') '</AppendedData>'
493  write(vtu,'(A)',advance='yes') '</VTKFile>'
494 
495  close(vtu)
496  endif
497 
498  return
499  end
500 !-----------------------------------------------------------------------
501  subroutine lpm_io_vtu_data(vtu,dataname,ncomp,idist)
502 
503  integer vtu,ncomp
504  integer*4 idist
505  character (len = *) dataname
506  character*50 dumstr
507 
508  write(vtu,'(A)',advance='no') ' <DataArray '
509  write(vtu,'(A)',advance='no') 'type="Float32" '
510  write(vtu,'(A)',advance='no') 'Name="'
511  write(vtu,'(A)',advance='no') dataname
512  write(vtu,'(A)',advance='no') '" NumberOfComponents="'
513  write(vtu,'(I0)',advance='no') ncomp
514  write(vtu,'(A)',advance='no') '" format="append" '
515  write(vtu,'(A)',advance='no') 'offset="'
516  write(vtu,'(I0)',advance='no') idist
517  write(vtu,'(A)',advance='yes') '"/>'
518 
519  return
520  end
521 !-----------------------------------------------------------------------
522  subroutine lpm_io_read_restart(filename)
523  include 'SIZE'
524  include 'SOLN'
525  include 'INPUT'
526  include 'MASS'
527  include 'GEOM'
528  include 'TSTEP'
529  include 'PARALLEL'
530 # include "LPM"
531 
532  parameter(lpm_lrf=3+lpm_lrp+lpm_lrs)
533  COMMON /lpm_rfpts_c/ lpm_rfpts
534  real*4 lpm_rfpts(lpm_lrf,lpm_lpart)
535 
536  common /nekmpi/ mid,mp,nekcomm,nekgroup,nekreal
537 
538  character*12 filename
539 
540  logical partl ! This is a dummy placeholder, used in cr()
541 
542  integer vtu,vtu1,prevs(2,mp)
543  integer*4 iint
544  integer*8 stride_len
545 
546  integer*8 disp
547  integer*4 nptot
548  integer pth, nread
549  real*4 rnptot
550 
551 ! ------------------------------
552 ! LOAD TOTAL NUMBER OF PARTICLES
553 ! ------------------------------
554  pth = 185
555 
556  write(6,*) filename
557 
558  disp = 0
559  icount = 0
560  if (nid .eq. 0) icount = 1
561  iorank = -1
562  call byte_open_mpi(trim(filename),pth,.true.,ierr)
563  call byte_set_view(disp,pth)
564  call byte_read_mpi(rnptot,icount,iorank,pth,ierr)
565  call byte_close_mpi(pth,ierr)
566 
567  nptot = int(rnptot)
568  call bcast(nptot, isize)
569 
570 ! --------------------------------------------------
571 ! FIRST GET HOW MANY PARTICLES WERE BEFORE THIS RANK
572 ! --------------------------------------------------
573  npmax = min(nptot/lpm_lpart+1,mp)
574  stride_len = 0
575  if (nid .le. npmax-1 .and. nid. ne. 0) stride_len = nid*lpm_lpart
576 
577  nread = lpm_lpart
578  if (nid .gt. npmax-1) nread = 0
579 
580  ndiff = nptot - (npmax-1)*lpm_lpart
581  if (nid .eq. npmax-1) nread = ndiff
582 
583 ! -------------------------
584 ! Parallel MPI file read in
585 ! -------------------------
586  disp = isize + stride_len*lpm_lrf*isize ! is there real*4 var?
587  icount = nread*lpm_lrf
588  iorank = -1
589  call byte_open_mpi(filename,pth,.true.,ierr)
590  call byte_set_view(disp,pth)
591  call byte_read_mpi(lpm_rfpts,icount,iorank,pth,ierr)
592  call byte_close_mpi(pth,ierr)
593 
594 ! ----------------------------------------
595 ! Assign values read in to rpart and ipart
596 ! ----------------------------------------
597  i = lpm_npart
598  do ii = 1,nread
599  i = lpm_npart + ii
600 
601  ic = 1
602  lpm_iprop(5,i) = int(lpm_rfpts(ic,ii))
603  ic = ic + 1
604  lpm_iprop(6,i) = int(lpm_rfpts(ic,ii))
605  ic = ic + 1
606  lpm_iprop(7,i) = int(lpm_rfpts(ic,ii))
607 
608  do j=1,lpm_lrs
609  ic = ic + 1
610  lpm_y(j,i) = lpm_rfpts(ic,ii)
611  enddo
612  do j=1,lpm_lrp
613  ic = ic + 1
614  lpm_rprop(j,i) = lpm_rfpts(ic,ii)
615  enddo
616  enddo
617  lpm_npart = i
618 
619  call lpm_comm_findpts
620  call lpm_comm_crystal
621 
622  return
623  end
624 c----------------------------------------------------------------------
subroutine byte_open_mpi(fnamei, mpi_fh, ifro, ierr)
Definition: byte_mpi.f:13
subroutine byte_write_mpi(buf, icount, iorank, mpi_fh, ierr)
Definition: byte_mpi.f:65
subroutine byte_read_mpi(buf, icount, iorank, mpi_fh, ierr)
Definition: byte_mpi.f:48
subroutine byte_close_mpi(mpi_fh, ierr)
Definition: byte_mpi.f:84
subroutine byte_set_view(ioff_in, mpi_fh)
Definition: byte_mpi.f:97
subroutine bcast(buf, len)
Definition: comm_mpi.f:431
real *8 function dnekclock_sync()
Definition: comm_mpi.f:401
subroutine lpm_comm_crystal
Definition: lpm_comm.f:136
subroutine lpm_comm_findpts
Definition: lpm_comm.f:19
subroutine lpm_io_write(filein1, iobig)
Definition: lpm_io.f:175
subroutine lpm_io_read_restart(filename)
Definition: lpm_io.f:523
subroutine lpm_io_vtu_data(vtu, dataname, ncomp, idist)
Definition: lpm_io.f:502
subroutine lpm_io_read(filein1, npart)
Definition: lpm_io.f:3
subroutine iadd(i1, iscal, n)
Definition: math.f:337
function iglsum(a, n)
Definition: math.f:926