IO: Implement write_snapshot_xdmf().

This subroutine stores an XDMF file which is a wrapper to HDF5 files for
easy reading simulation data to visualization tools like Paraview or
Visit.

This subroutine is based on the one written by Pierre Kestener.

Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
This commit is contained in:
Grzegorz Kowal 2016-12-21 12:15:32 -02:00
parent 723a20054b
commit 5a9d8f0c46

View File

@ -458,6 +458,7 @@ module io
! import external variables
!
use evolution , only : time
use mpitools , only : master
! local variables are not implicit by default
!
@ -479,6 +480,9 @@ module io
! store variable snapshot file
!
call write_snapshot_h5()
if (with_xdmf) then
call write_snapshot_xdmf()
end if
#endif /* HDF5 */
! increase the snapshot counter and calculate the next snapshot time
@ -7264,6 +7268,235 @@ module io
!-------------------------------------------------------------------------------
!
end subroutine read_5d_array_double_h5
!
!===============================================================================
!
! subroutine WRITE_SNAPSHOT_XDMF:
! ------------------------------
!
! Subroutine writes an XDMF file per snapshot per MPI process.
! XDMF file is just a wrapper which helps to load data in a visualization
! tools like Paraview or Visit.
!
! Based on the subroutine by Pierre Kestener
! (see https://bitbucket.org/pkestene/amun-code).
!
!
!===============================================================================
!
subroutine write_snapshot_xdmf()
! import procedures and variables from other modules
!
use blocks , only : block_data, list_data
use blocks , only : get_dblocks
use equations , only : nv, pvars
use mpitools , only : nproc
use coordinates , only : in, jn, kn
use coordinates , only : adx, ady, adz
use coordinates , only : ng
use evolution , only : time
! local variables are not implicit by default
!
implicit none
! local pointers
!
type(block_data), pointer :: pdata
! local variables
!
character(len=64) :: fname, hname
character(len=128) :: stmp, ttmp, sdim, bdim, pdim
integer(kind=4) :: l, p
integer(kind=4) :: ip, jp, kp
! local arrays
!
integer, dimension(12) :: slab
! local parameters
!
integer, parameter :: xdmf = 101
!
!-------------------------------------------------------------------------------
!
! prepare the XDMF and HDF5 file names
write(fname, "(a1,i6.6,'_',i5.5,'.xdmf')") ftype, isnap, nproc
write(hname, "(a1,i6.6,'_',i5.5,'.h5' )") ftype, isnap, nproc
! open the XDMF file
!
open (unit = xdmf, file = fname, status = 'replace')
! write the header
!
write(xdmf, "(a)") '<?xml version="1.0" encoding="UTF-8"?>'
write(xdmf, "(a)") '<Xdmf Version="2.2"' &
// ' xmlns:xi="http://www.w3.org/2003/XInclude">'
write(xdmf, "(a)") ' <Domain>'
write(stmp, "(a,i5.5)") 'region_', nproc
write(xdmf, "(a)") ' <Grid Name="' // trim(adjustl(stmp)) &
// '" GridType="Collection" CollectionType="Spatial">'
write(stmp, "(1g15.8)") time
write(xdmf, "(a)") ' <Time TimeType="Single" Value="' &
// trim(adjustl(stmp)) // '"/>'
! do not proceed if there are not data block on this processor
!
if (get_dblocks() > 0) then
! prepare dimensions
!
ip = in + 1
jp = jn + 1
#if NDIMS == 3
kp = kn + 1
#endif /* NDIMS == 3 */
#if NDIMS == 3
write(stmp, "(1i8)") kn
write(ttmp, "(1i8)") jn
stmp = trim(adjustl(stmp)) // ' ' // trim(adjustl(ttmp))
write(ttmp, "(1i8)") in
#else /* NDIMS == 3 */
write(stmp, "(1i8)") jn
write(ttmp, "(1i8)") in
#endif /* NDIMS == 3 */
bdim = trim(adjustl(stmp)) // ' ' // trim(adjustl(ttmp))
#if NDIMS == 3
write(stmp, "(1i8)") 1
#else /* NDIMS == 3 */
write(stmp, "(1i8)") kn
#endif /* NDIMS == 3 */
write(ttmp, "(1i8)") jn
stmp = trim(adjustl(stmp)) // ' ' // trim(adjustl(ttmp))
write(ttmp, "(1i8)") in
stmp = trim(adjustl(stmp)) // ' ' // trim(adjustl(ttmp))
write(ttmp, "(1i8)") get_dblocks()
sdim = trim(adjustl(stmp)) // ' ' // trim(adjustl(ttmp))
! prepare slab indices
!
#if NDIMS == 3
if (with_ghosts) then
slab(:) = (/ ng, ng, ng, -1, 1, 1, 1, 1, kn, jn, in, 1 /)
else
slab(:) = (/ 0, 0, 0, -1, 1, 1, 1, 1, kn, jn, in, 1 /)
end if
#else /* NDIMS == 3 */
if (with_ghosts) then
slab(:) = (/ 0, ng, ng, -1, 1, 1, 1, 1, 1, jn, in, 1 /)
else
slab(:) = (/ 0, 0, 0, -1, 1, 1, 1, 1, 1, jn, in, 1 /)
end if
#endif /* NDIMS == 3 */
! iterate over all data blocks
!
l = 0
pdata => list_data
do while(associated(pdata))
! store block geometry information
!
write(stmp, "(1i9.9)") pdata%meta%id
write(xdmf, "(a)") ' <Grid Name="block_' &
// trim(adjustl(stmp)) // '">'
#if NDIMS == 3
write(stmp, "(1i8)") kp
write(ttmp, "(1i8)") jp
stmp = trim(adjustl(stmp)) // ' ' // trim(adjustl(ttmp))
write(ttmp, "(1i8)") ip
stmp = trim(adjustl(stmp)) // ' ' // trim(adjustl(ttmp))
write(xdmf, "(a)") ' <Topology TopologyType="3DCoRectMesh"' &
// ' Dimensions="' // trim(adjustl(stmp)) // '"/>'
write(xdmf, "(a)") ' <Geometry GeometryType="ORIGIN_DXDYDZ">'
write(stmp, "(3e16.8)") pdata%meta%zmin, pdata%meta%ymin &
, pdata%meta%xmin
write(xdmf, "(a)") ' <DataItem Name="Origin" NumberType="Float"' &
// ' Precision="4" Dimensions="3" Format="XML">' &
// trim(adjustl(stmp)) // '</DataItem>'
write(stmp, "(3e16.8)") adz(pdata%meta%level), ady(pdata%meta%level) &
, adx(pdata%meta%level)
write(xdmf, "(a)") ' <DataItem Name="Spacing" NumberType="Float"' &
// ' Precision="4" Dimensions="3" Format="XML">' &
// trim(adjustl(stmp)) // '</DataItem>'
#else /* NDIMS == 3 */
write(stmp, "(1i8)") jp
write(ttmp, "(1i8)") ip
stmp = trim(adjustl(stmp)) // ' ' // trim(adjustl(ttmp))
write(xdmf, "(a)") ' <Topology TopologyType="2DCoRectMesh"' &
// ' Dimensions="' // trim(adjustl(stmp)) // '"/>'
write(xdmf, "(a)") ' <Geometry GeometryType="ORIGIN_DXDY">'
write(stmp, "(2e16.8)") pdata%meta%ymin, pdata%meta%xmin
write(xdmf, "(a)") ' <DataItem Name="Origin" NumberType="Float"' &
// ' Precision="4" Dimensions="2" Format="XML">' &
// trim(adjustl(stmp)) // '</DataItem>'
write(stmp, "(2e16.8)") ady(pdata%meta%level), adx(pdata%meta%level)
write(xdmf, "(a)") ' <DataItem Name="Spacing" NumberType="Float"' &
// ' Precision="4" Dimensions="2" Format="XML">' &
// trim(adjustl(stmp)) // '</DataItem>'
#endif /* NDIMS == 3 */
write(xdmf, "(a)") ' </Geometry>'
! convert slab dimensions to string
!
slab(4) = l
write(pdim, "(1i8)") slab(1)
do p = 2, 12
write(ttmp, "(1i8)") slab(p)
pdim = trim(adjustl(pdim)) // ' ' // trim(adjustl(ttmp))
end do
! loop over all variables and store their information
!
do p = 1, nv
write(xdmf, "(a)") ' <Attribute Name="' &
// trim(adjustl(pvars(p))) &
// '" AttributeType="Scalar" Center="Cell">'
write(xdmf, "(a)") ' <DataItem ItemType="Hyperslab"' &
// ' Dimensions="' // trim(adjustl(bdim)) &
// '" Type="HyperSlab">'
write(xdmf, "(a)") ' <DataItem Dimensions="3 4" Format="XML">' &
// trim(adjustl(pdim)) // '</DataItem>'
write(ttmp, "(a,':/variables/',a)") trim(hname), trim(pvars(p))
write(xdmf, "(a)") ' <DataItem NumberType="Float"' &
// ' Precision="8" Dimensions="' &
// trim(adjustl(sdim)) // '" Format="HDF">' &
// trim(adjustl(ttmp)) // '</DataItem>'
write(xdmf, "(a)") ' </DataItem>'
write(xdmf, "(a)") ' </Attribute>'
end do
! close grid structure for the current block
!
write(xdmf,"(a)") ' </Grid>'
l = l + 1
pdata => pdata%next
end do
end if ! get_dblocks() > 0
! close the XDMF structures
!
write(xdmf, "(a)") ' </Grid>'
write(xdmf, "(a)") ' </Domain>'
write(xdmf, "(a)") '</Xdmf>'
! close the XDMF file
!
close(xdmf)
!-------------------------------------------------------------------------------
!
end subroutine write_snapshot_xdmf
#endif /* HDF5 */
!===============================================================================