fileio_vtk.f90
Go to the documentation of this file.
1!#############################################################################
2!# #
3!# fosite - 3D hydrodynamical simulation program #
4!# module: fileio_vtk.f90 #
5!# #
6!# Copyright (C) 2010-2024 #
7!# Björn Sperling <sperling@astrophysik.uni-kiel.de> #
8!# Tobias Illenseer <tillense@astrophysik.uni-kiel.de> #
9!# Jannes Klee <jklee@astrophysik.uni-kiel.de> #
10!# #
11!# This program is free software; you can redistribute it and/or modify #
12!# it under the terms of the GNU General Public License as published by #
13!# the Free Software Foundation; either version 2 of the License, or (at #
14!# your option) any later version. #
15!# #
16!# This program is distributed in the hope that it will be useful, but #
17!# WITHOUT ANY WARRANTY; without even the implied warranty of #
18!# MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE, GOOD TITLE or #
19!# NON INFRINGEMENT. See the GNU General Public License for more #
20!# details. #
21!# #
22!# You should have received a copy of the GNU General Public License #
23!# along with this program; if not, write to the Free Software #
24!# Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. #
25!# #
26!#############################################################################
27!----------------------------------------------------------------------------!
51!----------------------------------------------------------------------------!
52MODULE fileio_vtk_mod
53 USE fileio_base_mod
54 USE fileio_gnuplot_mod
55 USE geometry_base_mod
56 USE mesh_base_mod
57 USE physics_base_mod
58 USE fluxes_base_mod
59 USE timedisc_base_mod
60 USE sources_generic_mod
61 USE common_dict
62#ifdef PARALLEL
63#ifdef HAVE_MPI_MOD
64 USE mpi
65#endif
66#endif
67 IMPLICIT NONE
68#ifdef PARALLEL
69#ifdef HAVE_MPIF_H
70 include 'mpif.h'
71#endif
72#endif
73 !--------------------------------------------------------------------------!
74 PRIVATE
75 ! Private Attributes section starts here:
78 INTEGER, PARAMETER :: maxcomp = 9
80 INTEGER, PARAMETER :: maxcols = 64
81 CHARACTER, PARAMETER :: lf = achar(10)
84#ifdef PARALLEL
85 INTEGER, PARAMETER :: rank_str_len = 5
86#endif
87
88! CHARACTER(LEN=16),DIMENSION(6),PARAMETER :: fluxkey = (/'bflux_WEST ', &
89! 'bflux_EAST ', &
90! 'bflux_SOUTH ', &
91! 'bflux_NORTH ', &
92! 'bflux_BOTTOM', &
93! 'bflux_TOP ' /)
94 !--------------------------------------------------------------------------!
95#ifdef PARALLEL
96
97 TYPE, EXTENDS(filehandle_fortran) :: filehandle_vts
99 CONTAINS
101 PROCEDURE :: initfilehandle
102 PROCEDURE :: getbasename
103 PROCEDURE :: getrankstring
104! FINAL :: Finalize_vts
105 END TYPE filehandle_vts
106#endif
107
108 TYPE, EXTENDS(fileio_gnuplot) :: fileio_vtk
110 TYPE(filehandle_fortran) :: pvdfile
111#ifdef PARALLEL
112 TYPE(filehandle_fortran) :: pvtsfile
113#endif
114 CHARACTER(LEN=32) :: buf
115 CHARACTER(LEN=32) :: endianness
116 CHARACTER(LEN=12) :: realfmt
117 INTEGER :: realsize
118 INTEGER :: ip1,jp1,kp1
119 REAL, DIMENSION(:,:,:,:), POINTER :: &
120 binout
121 REAL, DIMENSION(:,:) , POINTER :: &
122 bflux
123 REAL, DIMENSION(:,:,:,:), POINTER :: &
124 vtkcoords
125#ifdef PARALLEL
126 CHARACTER(LEN=64), DIMENSION(:), POINTER :: &
127 extent
128#endif
129 CONTAINS
131 PROCEDURE :: initfileio_deferred => initfileio_vtk
132 PROCEDURE :: writeheader
133! PROCEDURE :: ReadHeader
134 PROCEDURE :: writedataset_deferred => writedataset_vtk
135 PROCEDURE :: writeparaviewfile
136 final :: finalize
137 END TYPE
138 !--------------------------------------------------------------------------!
139 PUBLIC :: &
140 fileio_vtk
141
142
143CONTAINS
144
149 SUBROUTINE initfileio_vtk(this,Mesh,Physics,Timedisc,Sources,config,IO)
150 IMPLICIT NONE
151 !------------------------------------------------------------------------!
152 CLASS(fileio_vtk), INTENT(INOUT) :: this
153 CLASS(mesh_base), INTENT(IN) :: Mesh
154 CLASS(physics_base), INTENT(IN) :: Physics
155 CLASS(timedisc_base),INTENT(IN) :: Timedisc
156 CLASS(sources_list), ALLOCATABLE, INTENT(IN) :: Sources
157 TYPE(dict_typ), INTENT(IN), POINTER :: config
158 TYPE(dict_typ), INTENT(IN), POINTER :: IO
159 !------------------------------------------------------------------------!
160 TYPE(output_typ),DIMENSION(:),POINTER :: poutput
161 TYPE(real_t) :: dummy1
162 TYPE(dict_typ), POINTER :: node
163 CHARACTER(LEN=MAX_CHAR_LEN),DIMENSION(2) :: skip
164 INTEGER :: k,i,j
165#ifdef PARALLEL
166 INTEGER, DIMENSION(:), POINTER :: sendbuf,recvbuf
167 INTEGER :: n
168#endif
169 REAL, DIMENSION(:,:,:,:,:), POINTER :: corners
170 !------------------------------------------------------------------------!
171#ifdef PARALLEL
172 ALLOCATE(filehandle_vts::this%datafile)
173#endif
174 ! start init form base class in beginning
175 CALL this%InitFileIO(mesh,physics,timedisc,sources,config,io,"VTK","vts",textfile=.false.)
176
177 ! some sanity checks
178 IF (this%datafile%cycles.NE.this%count+1) &
179 CALL this%Error('fileio_vtk::InitFileIO','filecycles = count+1 required for VTK output')
180
181 ! determine dimensions for the coordinate array
182 IF (abs(mesh%xmax-mesh%xmin).LT.2*epsilon(mesh%xmin) &
183 .OR.(mesh%INUM.EQ.1.AND.mesh%ROTSYM.EQ.1)) THEN
184 ! suppress x-direction in output
185 this%IP1 = 0
186 ELSE
187 this%IP1 = 1
188 END IF
189 IF (abs(mesh%ymax-mesh%ymin).LT.2*epsilon(mesh%ymin) &
190 .OR.(mesh%JNUM.EQ.1.AND.mesh%ROTSYM.EQ.2)) THEN
191 ! suppress y-direction in output
192 this%JP1 = 0
193 ELSE
194 this%JP1 = 1
195 END IF
196 IF (abs(mesh%zmax-mesh%zmin).LT.2*epsilon(mesh%zmin) &
197 .OR.(mesh%KNUM.EQ.1.AND.mesh%ROTSYM.EQ.3)) THEN
198 ! suppress z-direction in output
199 this%KP1 = 0
200 ELSE
201 this%KP1 = 1
202 END IF
203
204 this%MAXCOLS = maxcols
205
206 ! allocate memory for auxilliary data fields
207 ALLOCATE(this%binout(maxcomp, &
208 mesh%IMIN:mesh%IMAX, &
209 mesh%JMIN:mesh%JMAX, &
210 mesh%KMIN:mesh%KMAX), &
211 this%vtkcoords(3, &
212 mesh%IMIN:mesh%IMAX+this%IP1, &
213 mesh%JMIN:mesh%JMAX+this%JP1, &
214 mesh%KMIN:mesh%KMAX+this%KP1), &
215 this%output(this%MAXCOLS), &
216 this%tsoutput(this%MAXCOLS), &
217 this%bflux(physics%VNUM,6), &
218 stat = this%err)
219 IF (this%err.NE.0) &
220 CALL this%Error("fileio_vtk::InitFileio_vtk", "Unable to allocate memory.")
221
222 ! determine size in bits of default real numbers
223 ! Fortran 2008 standard function!
224 this%realsize = storage_size(this%binout)
225 SELECT CASE(this%realsize)
226 CASE(32,64,128)
227 ! generate string with default real number format (size in bits)
228 WRITE(this%linebuf,fmt='(I4)',iostat=this%err) this%realsize
229 WRITE(this%realfmt,fmt='(A)',iostat=this%err) '"Float' // trim(adjustl(this%linebuf)) // '"'
230 ! size of default real numbers in byte
231 this%realsize = this%realsize/8
232 CASE DEFAULT
233 CALL this%Error("fileio_vtk::InitFileio_vtk", &
234 "only single, double or quadruple precision real numbers (4,8,16 bytes) supported")
235 END SELECT
236
237 ! determine the system endianness
238 CALL this%GetEndianness(this%endianness,'"LittleEndian"','"BigEndian"')
239
240#ifdef PARALLEL
241 ! allocate memory for this%extent, send & receive buffers;
242 ! recvbuf and this%extent are only needed on the receiver node (rank 0)
243 ! but parallel profiling using scalasca with mpich3 crashes
244 ! if recvbuf is not allocated on each node hence we allocate a minimal
245 ! array of size 1 on all nodes except for the receiver node;
246 ! remark: according to the MPI-2 standard recvbuf must only be allocated
247 ! on the receiver node
248 IF (this%GetRank() .EQ. 0 ) THEN
249 n = this%GetNumProcs()
250 ALLOCATE(this%extent(0:(n-1)),sendbuf(6),recvbuf(6*n),stat = this%err)
251 ELSE
252 ! see comment above
253 ALLOCATE(sendbuf(6),recvbuf(1),stat = this%err)
254 END IF
255 IF (this%err.NE.0) &
256 CALL this%Error( "fileio_vtk::InitFileio_vtk", "Unable to allocate memory.")
257
258 ! store information about grid extent on each node
260 sendbuf(1) = mesh%IMIN
261 sendbuf(2) = mesh%IMAX+this%IP1
262 sendbuf(3) = mesh%JMIN
263 sendbuf(4) = mesh%JMAX+this%JP1
264 sendbuf(5) = mesh%KMIN
265 sendbuf(6) = mesh%KMAX+this%KP1
266
267 ! gather information about grid extent on each node at the rank 0 node
268 CALL mpi_gather(sendbuf, 6, mpi_integer, recvbuf, 6, mpi_integer, &
269 0, mpi_comm_world, this%err)
270
271 ! store information about grid extent at the rank 0 node
272 IF (this%GetRank() .EQ. 0 ) THEN
273 DO i=0,n-1
274 WRITE(this%extent(i),fmt='(6(I7))',iostat=this%err)&
275 recvbuf(i*6+1),recvbuf(i*6+2),recvbuf(i*6+3),recvbuf(i*6+4),recvbuf(i*6+5),recvbuf(i*6+6)
276 END DO
277 END IF
278 ! free buffer memory
279 DEALLOCATE(sendbuf,recvbuf,stat=this%err)
280 IF (this%err.NE.0) &
281 CALL this%Error( "fileio_vtk::InitFileio_vtk", "Deallocation of sendbuf/recvbuf failed.")
282#endif
283
284 ! check whether mesh coordinates should be cartesian or curvilinear
285 IF (this%cartcoords) THEN
286 corners => mesh%RemapBounds(mesh%cart%corners)
287 ELSE
288 corners => mesh%RemapBounds(mesh%curv%corners)
289 END IF
290
291 ! store mesh corners in appropriate array for VTK coordinate output
292 DO k=mesh%KMIN,mesh%KMAX+mesh%KP1
293 DO j=mesh%JMIN,mesh%JMAX+mesh%JP1
294 DO i=mesh%IMIN,mesh%IMAX+mesh%IP1
295 this%vtkcoords(1:3,i,j,k) = corners(i,j,k,1,1:3)
296 END DO
297 END DO
298 END DO
299
300 ! add coordinates in 1D/2D simulations with non-vanishing extent
301 IF (mesh%INUM.EQ.1.AND.this%IP1.EQ.1) THEN
302 ! no transport in x-direction, but x-extent > 0
303 i = mesh%IMIN
304 DO k=mesh%KMIN,mesh%KMAX+mesh%KP1
305 DO j=mesh%JMIN,mesh%JMAX+mesh%JP1
306 this%vtkcoords(1:3,i+1,j,k) = corners(i,j,k,2,1:3)
307 END DO
308 END DO
309 END IF
310
311 IF (mesh%JNUM.EQ.1.AND.this%JP1.EQ.1) THEN
312 ! no transport in y-direction, but y-extent > 0
313 j = mesh%JMIN
314 DO k=mesh%KMIN,mesh%KMAX+mesh%KP1
315 DO i=mesh%IMIN,mesh%IMAX+mesh%IP1
316 this%vtkcoords(1:3,i,j+1,k) = corners(i,j,k,3,1:3)
317 END DO
318 END DO
319 END IF
320
321 IF (mesh%KNUM.EQ.1.AND.this%KP1.EQ.1) THEN
322 ! no transport in z-direction, but z-extent > 0
323 k = mesh%KMIN
324 DO j=mesh%JMIN,mesh%JMAX+mesh%JP1
325 DO i=mesh%IMIN,mesh%IMAX+mesh%IP1
326 this%vtkcoords(1:3,i,j,k+1) = corners(i,j,k,5,1:3)
327 END DO
328 END DO
329 END IF
330
331 ! get pointer to simulation time to make it the first entry
332 ! in the time step data set
333 CALL getattr(io,"/timedisc/time",dummy1)
334 IF (ASSOCIATED(dummy1%p)) THEN
335 this%tsoutput(1)%val => dummy1%p
336 this%tsoutput(1)%key = "TIME"
337 this%TSCOLS = 1
338 ELSE
339 this%TSCOLS = 0
340 END IF
341
342 ! create list of output data arrays
343 node => io
344 k = 0
345 skip(1:2) = [CHARACTER(LEN=MAX_CHAR_LEN) :: "corners", "time"]
346 CALL this%GetOutputlist(mesh,node,k,this%TSCOLS,skip)
347
348 ! shrink this%output
349 poutput => this%output
350 ALLOCATE(this%output,source=poutput(1:k),stat=this%err)
351 IF (this%err.NE.0) &
352 CALL this%Error("fileio_gnuplot::InitFileIO","memory allocation failed for this%output")
353 DEALLOCATE(poutput)
354
355 ! set size of output arrays in bytes
356 DO k = 1, SIZE(this%output)
357 this%output(k)%numbytes = 0
358 DO i = 1, SIZE(this%output(k)%p)
359 this%output(k)%numbytes = this%output(k)%numbytes + SIZE(this%output(k)%p(i)%val)*this%realsize
360 END DO
361 END DO
362
363 ! create file handle for the pvd file; only on rank 0 in parallel mode
364#ifdef PARALLEL
365 IF (this%GetRank().EQ.0) THEN
366 CALL this%pvtsfile%InitFilehandle(this%datafile%filename,this%datafile%path,"pvts",&
367 textfile=.true.,onefile=.false.,cycles=this%count+1)
368#endif
369 CALL this%pvdfile%InitFilehandle(this%datafile%filename,this%datafile%path,"pvd",&
370 textfile=.true.,onefile=.true.,cycles=1)
371#ifdef PARALLEL
372 END IF
373#endif
374 END SUBROUTINE initfileio_vtk
375
378 SUBROUTINE writeparaviewfile(this)
379 IMPLICIT NONE
380 !------------------------------------------------------------------------!
381 CLASS(fileio_vtk), INTENT(INOUT) :: this
382 !------------------------------------------------------------------------!
383 INTEGER :: k
384 REAL :: ftime
385#ifdef PARALLEL
386 INTEGER :: i
387#endif
388 !------------------------------------------------------------------------!
389#ifdef PARALLEL
390 IF (this%GetRank().EQ.0) THEN
391#endif
392 CALL this%pvdfile%OpenFile(replace,this%step)
393
394 ! write vtk header
395 WRITE(unit=this%pvdfile%GetUnitNumber(),fmt='(A)',iostat=this%err) &
396 '<?xml version="1.0"?>'//lf &
397 //'<VTKFile type="Collection" version="0.1" byte_order=' &
398 //trim(this%endianness)//'>'//lf &
399 //repeat(' ',2)//'<Collection>'
400
401 ! write entries for each time step
402 DO k=0,this%step
403 ftime = this%stoptime*k/this%count
404 IF (this%err.EQ. 0) &
405 WRITE(unit=this%pvdfile%GetUnitNumber(),fmt='(A,ES12.5,A)',iostat=this%err) &
406 repeat(' ',4)//'<DataSet timestep="',ftime,'" part="0" file="' // &
407#ifdef PARALLEL
408 trim(this%pvtsfile%GetBasename(k))//&
409 "." // trim(this%pvtsfile%extension)//'"/>'
410#else
411 trim(this%datafile%GetBasename(k))//&
412 "." // trim(this%datafile%extension)//'"/>'
413#endif
414#ifdef PARALLEL
415 ! in parallel mode each node generates its own file, but the rank 0 node
416 ! creates the pvd and pvts files
417! DO i=0,this%GetNumProcs()-1
418! IF (this%err.EQ. 0) WRITE(this%pvtsfile%GetUnitNumber(),FMT='(A,E11.5,A,I4.4,A)',IOSTAT=this%err) &
419! REPEAT(' ',4)//'<DataSet timestep="',&
420! ftime,'" part="', i ,'" file="'//TRIM(this%pvtsfile%GetBasename(i))//'"/>'
421! END DO
422#else
423#endif
424 END DO
425
426 IF (this%err.EQ. 0) &
427 WRITE(unit=this%pvdfile%GetUnitNumber(),fmt='(A)',iostat=this%err) &
428 repeat(' ',2)//'</Collection>'//lf//'</VTKFile>'
429
430 CALL this%pvdfile%CloseFile(this%step)
431
432 IF (this%err.NE.0) CALL this%Error("fileio_vtk::WriteParaviewFile","cannot write pvd-file")
433
434#ifdef PARALLEL
435 END IF
436#endif
437 END SUBROUTINE writeparaviewfile
438
439
440! !> \public Specific routine to open a file for vtk I/O
441! !!
442! SUBROUTINE OpenFile(this,action,ftype)
443! IMPLICIT NONE
444! !------------------------------------------------------------------------!
445! CLASS(fileio_vtk), INTENT(INOUT) :: this !< \param [in,out] this fileio type
446! INTEGER, INTENT(IN) :: action !< \param [in] action mode of file access
447! CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: ftype !< \param [in] file type
448! !------------------------------------------------------------------------!
449! CHARACTER(LEN=32) :: sta,act,pos,fext
450! !------------------------------------------------------------------------!
451! this%err=1
452! SELECT CASE(action)
453! CASE(READONLY)
454! sta = 'OLD'
455! act = 'READ'
456! pos = 'REWIND'
457! CASE(READEND)
458! sta = 'OLD'
459! act = 'READ'
460! pos = 'APPEND'
461! CASE(REPLACE)
462! sta = 'REPLACE'
463! act = 'WRITE'
464! pos = 'REWIND'
465! CASE(APPEND)
466! sta = 'OLD'
467! act = 'READWRITE'
468! pos = 'APPEND'
469! CASE DEFAULT
470! CALL this%Error("OpenFile","Unknown access mode.")
471! END SELECT
472!
473! ! check file type to open
474! IF (PRESENT(ftype)) THEN
475! fext=TRIM(ftype)
476! ELSE
477! ! default
478! fext='vts'
479! END IF
480!
481! SELECT CASE(TRIM(fext))
482! CASE('vts')
483! ! open the VTS file
484! OPEN(this%unit, FILE=this%GetFilename(),STATUS=TRIM(sta), &
485! ACCESS = 'STREAM' , &
486! ACTION=TRIM(act),POSITION=TRIM(pos),IOSTAT=this%err)
487! #ifdef PARALLEL
488! CASE('pvts')
489! ! open PVTS file (only parallel mode)
490! IF (this%GetRank().EQ.0) THEN
491! this%extension='pvts' ! change file name extension
492! OPEN(this%unit+100, FILE=this%GetFilename(-1),STATUS=TRIM(sta), &
493! FORM='FORMATTED',ACTION=TRIM(act),POSITION=TRIM(pos),IOSTAT=this%err)
494! this%extension='vts' ! reset default extension
495! END IF
496! #endif
497! CASE('pvd')
498! OPEN(this%unit+200, FILE=TRIM(this%filename)//'.'//TRIM(fext),STATUS=TRIM(sta), &
499! FORM='FORMATTED',ACTION=TRIM(act),POSITION=TRIM(pos),IOSTAT=this%err)
500! CASE DEFAULT
501! CALL this%Error("OpenFile","Unknown file type")
502! END SELECT
503!
504! IF (this%err.GT.0) &
505! CALL this%Error("OpenFile","cannot open " // TRIM(this%extension) // "-file")
506! END SUBROUTINE OpenFile
507
508! !> \public Specific routine to close a file for vtk I/O
509! !!
510! SUBROUTINE CloseFile(this,ftype)
511! IMPLICIT NONE
512! !------------------------------------------------------------------------!
513! CLASS(fileio_vtk), INTENT(INOUT) :: this !< \param [in,out] this fileio type
514! CHARACTER(LEN=*), OPTIONAL :: ftype !< \param [in] file type
515! !------------------------------------------------------------------------!
516! CHARACTER(LEN=4) :: fext
517! !------------------------------------------------------------------------!
518! INTENT(IN) :: ftype
519! !------------------------------------------------------------------------!
520! this%err=1
521! ! check file type to open
522! IF (PRESENT(ftype)) THEN
523! fext=TRIM(ftype)
524! ELSE
525! ! default
526! fext='vts'
527! END IF
528!
529! ! close file depending on the file type
530! SELECT CASE(TRIM(fext))
531! CASE('vts')
532! CLOSE(this%unit,IOSTAT=this%err)
533! #ifdef PARALLEL
534! CASE('pvts')
535! IF (this%GetRank() .EQ. 0 ) THEN
536! CLOSE(this%unit+100,IOSTAT=this%err)
537! END IF
538! #endif
539! CASE('pvd')
540! CLOSE(this%unit+200,IOSTAT=this%err)
541! CASE DEFAULT
542! CALL this%Error("OpenFile","Unknown file type")
543! END SELECT
544! ! set default extension
545!
546! IF(this%err .NE. 0) &
547! CALL this%Error( "CloseFileIO", "cannot close "//TRIM(fext)//"-file")
548! END SUBROUTINE CloseFile
549
550! !> Extract a subkey of a string (obsolete)
551! !!
552! !! Extract a subkey from a string (key) of type "abcd/SUBKEY/efgh"
553! FUNCTION GetSubKey(key) RESULT(subkey)
554! IMPLICIT NONE
555! !------------------------------------------------------------------------!
556! CHARACTER(LEN=*) :: key
557! !------------------------------------------------------------------------!
558! CHARACTER(LEN=MAXKEY) :: subkey
559! !------------------------------------------------------------------------!
560! INTENT(IN) :: key
561! !------------------------------------------------------------------------!
562! ! format of key like: "abcd/SUBKEY/efgh"
563! subkey = key(SCAN(key,"/",.TRUE.)+1:)
564! END FUNCTION GetSubKey
565
568 SUBROUTINE writeheader(this,Mesh,Physics,Header,IO)
569 IMPLICIT NONE
570 !------------------------------------------------------------------------!
571 CLASS(fileio_vtk), INTENT(INOUT) :: this
572 CLASS(mesh_base), INTENT(IN) :: Mesh
573 CLASS(physics_base), INTENT(IN) :: Physics
574 TYPE(dict_typ), POINTER :: IO
575 TYPE(dict_typ), POINTER :: Header
576 !------------------------------------------------------------------------!
577 INTEGER :: k
578 !------------------------------------------------------------------------!
579 ! write header in vts files
580 WRITE(this%linebuf,fmt='(A,6(I7),A)') '<?xml version="1.0"?>' //lf// &
581 '<VTKFile type="StructuredGrid" version="0.1" byte_order='&
582 //trim(this%endianness)//'>'//lf//&
583 ' <StructuredGrid WholeExtent="',&
584 mesh%IMIN,mesh%IMAX+this%IP1,mesh%JMIN,mesh%JMAX+this%JP1, &
585 mesh%KMIN,mesh%KMAX+this%KP1,'">'//lf// &
586 ' <FieldData>'//lf// &
587 ' <InformationKey type="String" Name="fosite version" format="ascii">'//lf// &
588 ' '//trim(version)//lf// &
589 ' </InformationKey>'//lf
590 WRITE(unit=this%datafile%GetUnitNumber(),iostat=this%err) trim(this%linebuf)
591
592#ifdef PARALLEL
593 ! write header in pvts file (only rank 0)
594 IF (this%GetRank() .EQ. 0 ) THEN
595 CALL this%pvtsfile%OpenFile(replace,this%step)
596 WRITE(unit=this%pvtsfile%GetUnitNumber(),fmt='(A,6(I7),A)',iostat=this%err) '<?xml version="1.0"?>' //lf// &
597 '<VTKFile type="PStructuredGrid" version="0.1" byte_order='&
598 //trim(this%endianness)//'>'//lf//&
599 ' <PStructuredGrid WholeExtent="',&
600 1,mesh%INUM+this%IP1,1,mesh%JNUM+this%JP1,1,mesh%KNUM+this%KP1,'">'//lf//&
601 ' <PFieldData>'//lf// &
602 ' <InformationKey type="String" Name="fosite version" format="ascii">'//lf// &
603 ' '//trim(version)//lf//&
604 ' </InformationKey>'
605 IF (this%err.NE.0) &
606 CALL this%Error("fileio_vtk::WriteHeader","Writing pvts file header failed")
607 END IF
608#endif
609
610 ! loop over all time step data registred in GetOutputlist
611 DO k=1,this%TSCOLS
612 WRITE(this%linebuf,fmt='(A,ES14.7,A)') &
613 ' <DataArray type='//trim(this%realfmt)//' Name="'//trim(this%tsoutput(k)%key)// &
614 '" NumberOfTuples="1" format="ascii">'//lf// &
615 repeat(' ',8),this%tsoutput(k)%val,lf//&
616 ' </DataArray>'
617 IF (this%err.EQ.0) THEN
618#ifndef PARALLEL
619 ! write time step data to VTS file in serial mode
620 WRITE(unit=this%datafile%GetUnitNumber(),iostat=this%err) trim(this%linebuf) // lf
621#else
622 ! write time step data to PVTS file in parallel mode
623 IF (this%GetRank() .EQ. 0 ) THEN
624 WRITE(unit=this%pvtsfile%GetUnitNumber(),fmt='(A)',advance='NO',iostat=this%err) trim(this%linebuf)
625 CALL this%pvtsfile%CloseFile(this%step)
626 END IF
627#endif
628 END IF
629 END DO
630
631 IF (this%err.NE.0) CALL this%Error("fileio_vtk::WriteHeader","cannot write (P)VTS file header")
632 END SUBROUTINE writeheader
633
636! SUBROUTINE ReadHeader(this,success)
637! IMPLICIT NONE
638! !------------------------------------------------------------------------!
639! CLASS(fileio_vtk), INTENT(INOUT) :: this
640! LOGICAL :: success
641! !------------------------------------------------------------------------!
642! INTENT(OUT) :: success
643! !------------------------------------------------------------------------!
644! CALL this%OpenFile(READONLY,this%step)
645! success = .FALSE.
646! CALL this%CloseFile(this%step)
647! END SUBROUTINE ReadHeader
648
651 SUBROUTINE writedataset_vtk(this,Mesh,Physics,Fluxes,Timedisc,Header,IO)
652 IMPLICIT NONE
653 !------------------------------------------------------------------------!
654 CLASS(fileio_vtk), INTENT(INOUT) :: this
655 CLASS(mesh_base), INTENT(IN) :: Mesh
656 CLASS(physics_base), INTENT(INOUT) :: Physics
657 CLASS(fluxes_base), INTENT(IN) :: Fluxes
658 CLASS(timedisc_base), INTENT(IN) :: Timedisc
659 TYPE(dict_typ), POINTER :: Header
660 TYPE(dict_typ), POINTER :: IO
661 !------------------------------------------------------------------------!
662 INTEGER :: i,j,k,m,l
663 INTEGER :: n, offset
664#ifdef PARALLEL
665 CHARACTER(LEN=256) :: filename
666#endif
667 !------------------------------------------------------------------------!
668 ! this part is formerly from fileio_generic.f90
669 ! calculate boundary fluxes, if they were requested for write
670! IF(ASSOCIATED(Timedisc%bflux)) THEN
671! print *, "test"
672! DO k=1,4
673! Timedisc%bflux(:,k) = Fluxes%GetBoundaryFlux(Mesh,Physics,k)
674! END DO
675!#ifdef PARALLEL
676! CALL MPI_BCAST(Timedisc%bflux, Physics%VNUM*4, DEFAULT_MPI_REAL, 0, &
677! Mesh%comm_cart, ierror)
678!#endif
679! END IF
680
681 CALL this%WriteParaviewFile()
682
683 !------------------------------------------------------------------------!
684 ! REMARK: VTS data files are opened in unformated or stream mode,
685 ! hence we always write the formated data into a line buffer
686 ! and then write the line buffer to the data file.
687 ! PVTS and PVD files are opened in formated mode.
688 ! -----------------------------------------------------------------------!
689 ! 1. META data
690 ! -----------------------------------------------------------------------!
691 ! a) global information; currently only scalar real numbers supported
692#ifdef PARALLEL
693 IF (this%GetRank() .EQ. 0 ) THEN
694 CALL this%pvtsfile%OpenFile(append,this%step)
695 IF(this%err .NE. 0) CALL this%Error( "fileio_vtk::WriteDataset_vtk", "Can't open pvts file")
696 END IF
697#endif
698
700
701! simple hack to write positions of the binary star system
702! IF (HasKey(IO, "/sources/grav/binary/binpos/value"))THEN
703! CALL GetAttr(IO, "/sources/grav/binary/binpos", dummy2)
704!
705! WRITE(this%tslinebuf,fmt='(ES14.7,A,ES14.7)',IOSTAT=this%err)&
706! dummy2(1,1), repeat(' ',4),dummy2(2,1)
707!
708! WRITE(this%unit, IOSTAT=this%err)&
709! LF//repeat(' ',6)//'<DataArray type="Float64" Name="position primary"' &
710! //' NumberOfTuples="2" format="ascii">' &
711! //LF//repeat(' ',8)//TRIM(this%tslinebuf) &
712! //LF//repeat(' ',6)//'</DataArray>'
713!
714! #ifdef PARALLEL
715! IF (GetRank(this) .EQ. 0 ) &
716! WRITE(this%unit+100, IOSTAT=this%err)&
717! LF//repeat(' ',6)//'<DataArray type="Float64" Name="position primary"' &
718! //' NumberOfTuples="2" format="ascii">' &
719! //LF//repeat(' ',8)//TRIM(this%tslinebuf) &
720! //LF//repeat(' ',6)//'</DataArray>'
721!
722! #endif
723!
724! WRITE(this%tslinebuf,fmt='(ES14.7,A,ES14.7)',IOSTAT=this%err)&
725! dummy2(1,2), repeat(' ',4),dummy2(2,2)
726!
727! WRITE(this%unit, IOSTAT=this%err)&
728! LF//repeat(' ',6)//'<DataArray type="Float64" Name="position secondary"' &
729! //' NumberOfTuples="2" format="ascii">' &
730! //LF//repeat(' ',8)//TRIM(this%tslinebuf) &
731! //LF//repeat(' ',6)//'</DataArray>'
732!
733! #ifdef PARALLEL
734! IF (GetRank(this) .EQ. 0 ) &
735! WRITE(this%unit+100, IOSTAT=this%err)&
736! LF//repeat(' ',6)//'<DataArray type="Float64" Name="position secondary"' &
737! //' NumberOfTuples="2" format="ascii">' &
738! //LF//repeat(' ',8)//TRIM(this%tslinebuf) &
739! //LF//repeat(' ',6)//'</DataArray>'
740! #endif
741!
742! END IF
743
744 ! -----------------------------------------------------------------------!
745 ! b) coordinates
746 offset = 0
747 WRITE(this%linebuf,fmt='(A,(6(I7)),A,I8,A)')&
748 ' </FieldData>'//lf// &
749 ' <Piece Extent="', &
750 mesh%IMIN,mesh%IMAX+this%IP1,mesh%JMIN,mesh%JMAX+this%JP1,&
751 mesh%KMIN,mesh%KMAX+this%KP1,'">'//lf// &
752 ' <Points>'//lf//&
753 ' <DataArray type='//trim(this%realfmt)//' NumberOfComponents="3"'//&
754 ' Name="Point" format="appended" offset="',offset,'">'//lf// &
755 ' </DataArray>'//lf// &
756 ' </Points>'//lf// &
757 ' <CellData>'//lf
758 IF (this%err.EQ.0) &
759 WRITE(unit=this%datafile%GetUnitNumber(),iostat=this%err) trim(this%linebuf)
760 ! add size of data set (in bytes) + 4 (4 byte integer for size information)
761 ! -> "jump address" for next data set in the file in bytes
762 offset = offset + SIZE(this%vtkcoords)*this%realsize + 4
763
764#ifdef PARALLEL
765 IF (this%GetRank() .EQ. 0 ) THEN
766 WRITE(unit=this%pvtsfile%GetUnitNumber(),fmt='(A)',iostat=this%err)&
767 repeat(' ',4)//'</PFieldData>' &
768 //lf//repeat(' ',4)//'<PPoints>' &
769 //lf//repeat(' ',6)//'<DataArray type='//trim(this%realfmt) &
770 //' NumberOfComponents="3" Name="Point"/>' &
771 //lf//repeat(' ',4)//'</PPoints>' &
772 //lf//repeat(' ',4)//'<PCellData>'
773 END IF
774#endif
775
776 ! -----------------------------------------------------------------------!
777 ! c) output data fields
778 DO k = 1, SIZE(this%output)
779 WRITE(this%linebuf,fmt='(A,I3,A,I8,A)')&
780 repeat(' ',8)//'<DataArray type='//this%realfmt//&
781 ' NumberOfComponents="', SIZE(this%output(k)%p), &
782 '" Name="'//trim(this%output(k)%key)//&
783 '" format="appended" offset="',offset,'"/>'//lf
784 IF (this%err.EQ.0) &
785 WRITE(unit=this%datafile%GetUnitNumber(),iostat=this%err) trim(this%linebuf)
786
787#ifdef PARALLEL
788 IF (this%GetRank() .EQ. 0 ) THEN
789 WRITE(unit=this%pvtsfile%GetUnitNumber(),fmt='(A,I3,A)',iostat=this%err)&
790 repeat(' ',6)//'<DataArray type='//this%realfmt//&
791 ' NumberOfComponents="', SIZE(this%output(k)%p),&
792 '" Name="'//trim(this%output(k)%key)//'"/>'
793 END IF
794#endif
795 ! add size of data set (in bytes) + 4 (4 byte integer for size information)
796 ! -> "jump address" for next data set in the file in bytes
797 offset = offset + this%output(k)%numbytes + 4
798 END DO
799
800 IF (this%err.EQ.0) &
801 WRITE(unit=this%datafile%GetUnitNumber(),iostat=this%err) &
802 repeat(' ',6)//'</CellData>' &
803 //lf//repeat(' ',4)//'</Piece>' &
804 //lf//repeat(' ',2)//'</StructuredGrid>' &
805 //lf//repeat(' ',2)//'<GlobalData>'
806
807#ifdef PARALLEL
808 IF (this%GetRank() .EQ. 0 ) THEN
809 WRITE(unit=this%pvtsfile%GetUnitNumber(),fmt='(A)',iostat=this%err) repeat(' ',4)//'</PCellData>'
810
811 SELECT TYPE(df=>this%datafile)
812 CLASS IS(filehandle_vts)
813 DO i=0,this%GetNumProcs()-1
814 filename = trim(df%path) // &
815 trim(df%filename) // &
816 trim(df%GetRankString(i)) // &
817 trim(df%GetStepString(this%step)) // &
818 "." // trim(this%datafile%extension)
819 WRITE(unit=this%pvtsfile%GetUnitNumber(),fmt='(A)',iostat=this%err) &
820 repeat(' ',4)//'<Piece Extent="'//trim(this%extent(i))&
821 //'" Source="'//trim(filename)//'"/>'
822 END DO
823 CLASS DEFAULT
824 CALL this%Error("fileio_vtk::WriteDataset","datafile must be of type filehandle_vts in parallel mode")
825 END SELECT
826
827 WRITE(unit=this%pvtsfile%GetUnitNumber(),fmt='(A)',iostat=this%err)&
828 repeat(' ',2)//'</PStructuredGrid>' &
829 //lf//'</VTKFile>'
830
831 END IF
832#endif
833
834 ! -----------------------------------------------------------------------!
835 ! d) boundary fluxes
836! DO i=1,4
837! ! in PARALLEL mode the result is sent to process with rank 0
838! this%bflux(:,i) = GetBoundaryFlux(Fluxes,Mesh,Physics,i)
839!
840! WRITE(this%unit,IOSTAT=this%err)&
841! LF//repeat(' ',4)//'<DataArray type='//this%realfmt//&
842! 'Name="'//TRIM(fluxkey(i))//'" format="ascii" >' &
843! //LF//repeat(' ',6)
844! DO k=1,Physics%VNUM
845! WRITE(this%linebuf((k-1)*20+1:k*20),fmt='(E16.9,A)',IOSTAT=this%err)&
846! this%bflux(k,i),repeat(' ',4)
847! END DO
848!
849! WRITE(this%unit,IOSTAT=this%err)TRIM(this%linebuf(1:Physics%VNUM*20)) &
850! //LF//repeat(' ',4)//'</DataArray>'
851! END DO
852
853 ! end META data
854 ! -----------------------------------------------------------------------------------!
855
856 IF (this%err.EQ.0) &
857 WRITE(unit=this%datafile%GetUnitNumber(),iostat=this%err)&
858 lf//repeat(' ',2)//'</GlobalData>' &
859 //lf//' <AppendedData encoding="raw">' &
860 //lf//'_'
861 IF (this%err.NE.0) &
862 CALL this%Error("fileio_vtk::WriteDataset", "writing data array specs failed")
863
864#ifdef PARALLEL
865 IF (this%GetRank() .EQ. 0 ) THEN
866 CALL this%pvtsfile%CloseFile(this%step)
867 IF(this%err .NE. 0) CALL this%Error( "fileio_vtk::WriteDataset_vtk", "Can't close pvts file")
868 END IF
869#endif
870
871 ! -----------------------------------------------------------------------!
872 ! 2. REAL data
873 ! -----------------------------------------------------------------------!
874 ! a) coordinates:
875 ! (i) size of data array in bytes (must be a 4 byte integer)
876 ! (ii) coordinate data array (generated in InitFileio)
877 IF (this%err.EQ.0) &
878 WRITE(unit=this%datafile%GetUnitNumber(),iostat=this%err) &
879 int(SIZE(this%vtkcoords)*this%realsize + 4,kind=4),this%vtkcoords(:,:,:,:)
880 IF (this%err.NE.0) &
881 CALL this%Error("fileio_vtk::WriteDataset", "writing coordinates failed")
882
883 ! -----------------------------------------------------------------------!
884 ! b) data fields
885 ! (i) size of data array in bytes (must be a 4 byte integer);
886 ! this is determined in GetOutputlist and stored in
887 ! this%output(:)%numbytes
888 ! (ii) the data given by a list of pointers which is also
889 ! generated in GetOutputlist
890 ! loop over all output data fields
891 DO l = 1, SIZE(this%output)
892 ! reorganize the data to comply with VTK style ordering
893 n = SIZE(this%output(l)%p)
894 DO m=1,n
895 DO k=1,mesh%KMAX-mesh%KMIN+1
896 DO j=1,mesh%JMAX-mesh%JMIN+1
897 DO i=1,mesh%IMAX-mesh%IMIN+1
898 this%binout(m,mesh%IMIN+i-1,mesh%JMIN+j-1,mesh%KMIN+k-1) &
899 = this%output(l)%p(m)%val(i,j,k)
900 END DO
901 END DO
902 END DO
903 END DO
904 IF (this%err.EQ.0) &
905 WRITE(unit=this%datafile%GetUnitNumber(),iostat=this%err) &
906 int(this%output(l)%numbytes,kind=4), &
907! this%binout(1:n,:,:,:)
908 ((((this%binout(m,i,j,k),m=1,n),i=mesh%IMIN,mesh%IMAX), &
909 j=mesh%JMIN,mesh%JMAX),k=mesh%KMIN,mesh%KMAX)
910
911 IF (this%err.NE.0) &
912 CALL this%Error("fileio_vtk::WriteDataset", "writing data")
913 END DO
914 WRITE(unit=this%datafile%GetUnitNumber(),iostat=this%err) lf//repeat(' ',2)//&
915 '</AppendedData>'//lf//'</VTKFile>'//lf
916 END SUBROUTINE writedataset_vtk
917
920 SUBROUTINE finalize(this)
921 IMPLICIT NONE
922 !------------------------------------------------------------------------!
923 TYPE(fileio_vtk), INTENT(INOUT) :: this
924 !------------------------------------------------------------------------!
925 ! REMARK: this%output and this%tsoutput are deallocated in the gnuplot finalizer
926 ! which is called automatically
927 IF (ASSOCIATED(this%binout)) DEALLOCATE(this%binout)
928 IF (ASSOCIATED(this%vtkcoords)) DEALLOCATE(this%vtkcoords)
929 IF (ASSOCIATED(this%bflux)) DEALLOCATE(this%bflux)
930 NULLIFY(this%binout,this%vtkcoords,this%bflux)
931 END SUBROUTINE finalize
932
933#ifdef PARALLEL
934
935 SUBROUTINE initfilehandle(this,filename,path,extension,textfile,onefile,cycles,unit)
936 IMPLICIT NONE
937 !-------------------------------------------------------------------!
938 CLASS(filehandle_vts), INTENT(INOUT) :: this
939 CHARACTER(LEN=*), INTENT(IN) :: filename
940 CHARACTER(LEN=*), INTENT(IN) :: path
941 CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: extension
942 LOGICAL, OPTIONAL, INTENT(IN) :: textfile
943 LOGICAL, OPTIONAL, INTENT(IN) :: onefile
944 INTEGER, OPTIONAL, INTENT(IN) :: cycles
945 INTEGER, OPTIONAL, INTENT(IN) :: unit
946 !-------------------------------------------------------------------!
947 CALL this%filehandle_fortran%InitFilehandle(filename,path,extension,textfile,onefile,cycles,unit)
948 ! check number of parallel processes
949 IF (rank_str_len.LT.3) &
950 CALL this%Error("filehandle_vts::InitFilehandle", &
951 "we need at least 3 bytes for rank string: 2 for '-r' + 1 for process number (0..9)")
952 IF (this%GetNumProcs().GT.10**(rank_str_len-3)) &
953 CALL this%Error("filehandle_vts::InitFilehandle", &
954 "number of processes for multiple file output exceeds limits")
955 END SUBROUTINE initfilehandle
956
958 FUNCTION getbasename(this,step) RESULT (fname)
959 IMPLICIT NONE
960 !------------------------------------------------------------------------!
961 CLASS(filehandle_vts), INTENT(IN) :: this
962 INTEGER, INTENT(IN) :: step
963 CHARACTER(LEN=256) :: fname
964 !------------------------------------------------------------------------!
965 fname = trim(this%path) // trim(this%filename) // trim(this%GetRankString()) &
966 // trim(this%GetStepString(step))
967 END FUNCTION getbasename
968
970 FUNCTION getrankstring(this,rank)
971 IMPLICIT NONE
972 !------------------------------------------------------------------------!
973 CLASS(filehandle_vts), INTENT(IN) :: this
974 INTEGER, OPTIONAL, INTENT(IN) :: rank
975 CHARACTER(LEN=RANK_STR_LEN) :: getrankstring
976 !------------------------------------------------------------------------!
977 CHARACTER(LEN=16) :: fmtstr
978 !-------------------------------------------------------------------!
979 ! process rank as string with "-r" + leading zeros
980 WRITE (fmtstr ,'(A,I1,A)') "(A2,I0.",rank_str_len-2,")"
981 IF (PRESENT(rank)) THEN
982 WRITE (getrankstring,fmt=trim(fmtstr)) "-r",rank
983 ELSE
984 WRITE (getrankstring,fmt=trim(fmtstr)) "-r",this%GetRank()
985 END IF
986 END FUNCTION getrankstring
987
988! !> \public Destructor of VTS file handle
989! !!
990! SUBROUTINE Finalize_vts(this)
991! IMPLICIT NONE
992! !------------------------------------------------------------------------!
993! TYPE(filehandle_vts), INTENT(INOUT) :: this !< \param [in,out] this filehandle
994! !------------------------------------------------------------------------!
995! END SUBROUTINE Finalize_vts
996#endif
997END MODULE fileio_vtk_mod
common_dict
Dictionary for generic data types.
Definition: common_dict.f90:61
common_dict::this
type(logging_base), save this
Definition: common_dict.f90:149
fileio_base_mod
Generic file I/O module.
Definition: fileio_base.f90:54
fileio_base_mod::append
integer, parameter append
read/write access at end
Definition: fileio_base.f90:277
fileio_base_mod::initfilehandle
subroutine initfilehandle(this, filename, path, extension, textfile, onefile, cycles, unit)
basic initialization of Fortran file handle
Definition: fileio_base.f90:331
fileio_base_mod::replace
integer, parameter replace
read/write access replacing file
Definition: fileio_base.f90:276
fileio_base_mod::getbasename
character(len=256) function getbasename(this, step)
get the file name without path and extension but with step string appended; if thisonefile suppress t...
Definition: fileio_base.f90:787
fileio_gnuplot_mod
I/O for GNUPLOT readable tabular files.
Definition: fileio_gnuplot.f90:44
fileio_gnuplot_mod::finalize
subroutine finalize(this)
Closes the file I/O.
Definition: fileio_gnuplot.f90:814
fileio_gnuplot_mod::writeheader
subroutine writeheader(this, Mesh, Physics, Header, IO)
Writes the configuration as a header to the file.
Definition: fileio_gnuplot.f90:545
fileio_gnuplot_mod::lf
character, parameter lf
line feed
Definition: fileio_gnuplot.f90:73
fileio_vtk_mod
I/O for VTK files in XML format (vtkStructuredGrid)
Definition: fileio_vtk.f90:52
fileio_vtk_mod::writedataset_vtk
subroutine writedataset_vtk(this, Mesh, Physics, Fluxes, Timedisc, Header, IO)
Reads the header (not yet implemented)
Definition: fileio_vtk.f90:652
fileio_vtk_mod::initfileio_vtk
subroutine initfileio_vtk(this, Mesh, Physics, Timedisc, Sources, config, IO)
Constructor for the VTK file I/O.
Definition: fileio_vtk.f90:150
fileio_vtk_mod::maxcols
integer, parameter maxcols
max of different output arrays
Definition: fileio_vtk.f90:80
fileio_vtk_mod::rank_str_len
integer, parameter rank_str_len
length of multi process rank string
Definition: fileio_vtk.f90:85
fileio_vtk_mod::writeparaviewfile
subroutine writeparaviewfile(this)
Write the Paraview global description file.
Definition: fileio_vtk.f90:379
fileio_vtk_mod::maxcomp
integer, parameter maxcomp
max. of allowed components 9 is a tensor (rank 2, dim 3)
Definition: fileio_vtk.f90:78
fileio_vtk_mod::getrankstring
character(len=rank_str_len) function getrankstring(this, rank)
convert process rank to string
Definition: fileio_vtk.f90:971
fluxes_base_mod
base module for numerical flux functions
Definition: fluxes_base.f90:39
geometry_base_mod
base class for geometrical properties
Definition: geometry_base.f90:36
mesh_base_mod
basic mesh module
Definition: mesh_base.f90:72
physics_base_mod
Basic physics module.
Definition: physics_base.f90:55
sources_generic_mod
module to manage list of source terms
Definition: sources_generic.f90:37
timedisc_base_mod
Definition: timedisc_base.f90:57
common_dict::dict_typ
Definition: common_dict.f90:123
common_dict::real_t
Definition: common_dict.f90:131
fileio_base_mod::filehandle_fortran
class for Fortran file handle
Definition: fileio_base.f90:99
fileio_gnuplot_mod::fileio_gnuplot
FileIO gnuplot class.
Definition: fileio_gnuplot.f90:101
fileio_gnuplot_mod::output_typ
Definition: fileio_gnuplot.f90:87
fileio_vtk_mod::filehandle_vts
names of fluxes
Definition: fileio_vtk.f90:97
fileio_vtk_mod::fileio_vtk
FileIO class for VTK output.
Definition: fileio_vtk.f90:108
fluxes_base_mod::fluxes_base
Definition: fluxes_base.f90:61
mesh_base_mod::mesh_base
mesh data structure
Definition: mesh_base.f90:122
physics_base_mod::physics_base
Definition: physics_base.f90:69
sources_generic_mod::sources_list
container class to manage the list of source terms
Definition: sources_generic.f90:58
timedisc_base_mod::timedisc_base
Definition: timedisc_base.f90:84
Generated by   doxygen 1.9.4