fileio_xdmf.f90

Go to the documentation of this file.

!#############################################################################
!#                                                                           #
!# fosite - 3D hydrodynamical simulation program                             #
!# module: fileio_xdmf.f90                                                   #
!#                                                                           #
!# Copyright (C) 2013-2015                                                   #
!# Manuel Jung <mjung@astrophysik.uni-kiel.de>                               #
!# Jannes Klee <jklee@astrophysik.uni-kiel.de>                               #
!#                                                                           #
!# This program is free software; you can redistribute it and/or modify      #
!# it under the terms of the GNU General Public License as published by      #
!# the Free Software Foundation; either version 2 of the License, or (at     #
!# your option) any later version.                                           #
!#                                                                           #
!# This program is distributed in the hope that it will be useful, but       #
!# WITHOUT ANY WARRANTY; without even the implied warranty of                #
!# MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE, GOOD TITLE or        #
!# NON INFRINGEMENT.  See the GNU General Public License for more            #
!# details.                                                                  #
!#                                                                           #
!# You should have received a copy of the GNU General Public License         #
!# along with this program; if not, write to the Free Software               #
!# Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.                 #
!#                                                                           #
!#############################################################################
!----------------------------------------------------------------------------!
!----------------------------------------------------------------------------!
MODULE fileio_xdmf_mod
  USE fileio_base_mod
  USE fileio_binary_mod
  USE geometry_base_mod
  USE mesh_base_mod
  USE physics_base_mod
  USE timedisc_base_mod
  USE fluxes_base_mod
  USE sources_base_mod
  USE common_dict
  !--------------------------------------------------------------------------!
  PRIVATE
  CHARACTER, PARAMETER         :: lf = achar(10)        

  TYPE, EXTENDS(fileio_binary) :: fileio_xdmf
    CHARACTER(LEN=14)          :: endian_xdmf
    INTEGER                    :: unit_xdmf
  CONTAINS
    PROCEDURE :: initfileio_xdmf
    PROCEDURE :: iteratedict
    PROCEDURE :: writexmf
    PROCEDURE :: writenode_xdmf
    PROCEDURE :: writekey_xdmf
    PROCEDURE :: writeattribute
    PROCEDURE :: writemeshxml
    PROCEDURE :: writevector
    PROCEDURE :: writedataitem
    PROCEDURE :: finalize
  END TYPE
  !--------------------------------------------------------------------------!
  PUBLIC :: &
       ! types
       fileio_xdmf
       ! constants
       ! methods
  !--------------------------------------------------------------------------!

CONTAINS
  SUBROUTINE initfileio_xdmf(this,Mesh,Physics,Timedisc,Sources,config,IO)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(fileio_xdmf),   INTENT(INOUT) :: this
    CLASS(mesh_base),     INTENT(IN)    :: Mesh
    CLASS(physics_base),  INTENT(IN)    :: Physics
    CLASS(timedisc_base), INTENT(IN)    :: Timedisc
    CLASS(sources_base),  INTENT(IN), POINTER &
                                        :: Sources
    TYPE(Dict_TYP),       INTENT(IN), POINTER &
                                        :: config
    TYPE(Dict_TYP),       INTENT(IN), POINTER &
                                        :: IO
    !------------------------------------------------------------------------!
    INTEGER               :: err
    !------------------------------------------------------------------------!
    !\attention interface changed, however "xdmf" is not passed anymore
    CALL this%InitFileIO_binary(mesh,physics,timedisc,sources,config,io)

    CALL this%GetEndianness(this%endian_xdmf, '"Little"', '"Big"')
    IF(this%realsize.GT.8) &
      CALL this%Error("WriteXMF","Only single and double precision are allowed")

    WRITE(this%realfmt,'(A,I1,A)',iostat=err) '"',this%realsize,'"'

    CALL this%WriteXMF(mesh,io)

  END SUBROUTINE initfileio_xdmf


  RECURSIVE SUBROUTINE iteratedict(this,Mesh,config,offset,filename,path)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(fileio_xdmf), INTENT(INOUT)  :: this
    CLASS(mesh_base),   INTENT(IN)     :: mesh
    TYPE(dict_typ),     INTENT(IN), POINTER  &
                                       :: config
    CHARACTER(LEN=*),   INTENT(IN)     :: filename
    CHARACTER(LEN=*),   INTENT(IN), OPTIONAL &
                                       :: path
    INTEGER, INTENT(INOUT)             :: offset
    !------------------------------------------------------------------------!
    CHARACTER(LEN=MAX_CHAR_LEN)        :: str, key
    TYPE(dict_typ),POINTER             :: node
    !------------------------------------------------------------------------!
    IF(PRESENT(path)) THEN
        str = path
    ELSE
        str = ""
    ENDIF
    node => config
    DO WHILE(ASSOCIATED(node))
      key = trim(str)//"/"//trim(getkey(node))

      CALL this%WriteNode_xdmf(mesh,key,node,offset,filename)

      IF(haschild(node)) THEN
        CALL this%IterateDict(mesh, getchild(node), offset, filename, key)
      END IF
      node=>getnext(node)
    END DO
  END SUBROUTINE iteratedict

  SUBROUTINE writekey_xdmf(this,offset,key,type,bytes,dims)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(fileio_xdmf), INTENT(INOUT) :: this
    CHARACTER(LEN=*),   INTENT(IN)    :: key
    INTEGER,            INTENT(INOUT) :: offset
    INTEGER,            INTENT(IN)    :: type
    INTEGER,            INTENT(IN)    :: bytes
    INTEGER,DIMENSION(5), OPTIONAL, INTENT(IN) &
                                      :: dims
    !------------------------------------------------------------------------!
    INTEGER :: keylen,l
    !------------------------------------------------------------------------!
    keylen = len_trim(key)
    IF(PRESENT(dims)) THEN
      l = 3
      IF(dims(4).GT.1) l = 4
      IF(dims(5).GT.1) l = 5
    ELSE
      l = 0
    END IF
    offset = offset + (3+l) * this%intsize + keylen
  END SUBROUTINE writekey_xdmf


  SUBROUTINE writenode_xdmf(this,Mesh,key,node,offset,filename)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(fileio_xdmf), INTENT(INOUT) :: this
    CLASS(mesh_base),   INTENT(IN)    :: Mesh
    CHARACTER(LEN=MAX_CHAR_LEN), INTENT(IN) &
                                      :: key
    TYPE(Dict_TYP),     INTENT(IN), POINTER &
                                      :: node
    INTEGER ,           INTENT(INOUT) :: offset
    CHARACTER(LEN=*),   INTENT(IN)    :: filename
    !------------------------------------------------------------------------!
    REAL,DIMENSION(:,:,:),    POINTER :: ptr3 => null()
    REAL,DIMENSION(:,:,:,:),  POINTER :: ptr4 => null()
    REAL,DIMENSION(:,:,:,:,:),POINTER :: ptr5 => null()
    INTEGER,DIMENSION(5)              :: dims
    INTEGER                           :: bytes
    CHARACTER(LEN=1),DIMENSION(:),POINTER  &
                                      :: val
    LOGICAL                           :: ref
    INTEGER                           :: type
    !------------------------------------------------------------------------!
    ref = .false.
    dims(:) = 1

    type = getdatatype(node)
    SELECT CASE(type)
    CASE(dict_real_p)
      bytes = this%realsize
      type = dict_real
    CASE(dict_real_threed)
      CALL getattr(node,key,ptr3)
      dims(1:3) = shape(ptr3)
    CASE(dict_real_fourd)
      CALL getattr(node,key,ptr4)
      dims(1:4) = shape(ptr4)
    CASE(dict_real_fived)
      CALL getattr(node,key,ptr5)
      dims(1:5) = shape(ptr5)
    CASE DEFAULT
      val => getdata(node)
      IF(ASSOCIATED(val)) THEN
        bytes = SIZE(val)
      ELSE
        bytes = 0
      END IF
    END SELECT

    IF(product(dims).GT.1) THEN
      CALL this%SetMeshDims(mesh,dims)
      bytes = product(dims) * this%realsize
      CALL this%WriteKey_xdmf(offset,key,type,bytes,dims)
    ELSE
      CALL this%WriteKey_xdmf(offset,key,type,bytes)
    END IF

    SELECT CASE(getdatatype(node))
    CASE(dict_real_threed)
      CALL this%WriteAttribute(mesh,trim(key),dims(1:3),filename,offset,ref)
! It is not clear, how other datatypes than DICT_REAL_TWOD should be mapped in
! XDMF.
!    CASE(DICT_REAL_THREED)
!      CALL this%WriteAttribute(Mesh,TRIM(key),dims(1:3),filename,offset,ref)
!    CASE(DICT_REAL_FOURD)
!      CALL this%WriteAttribute(Mesh,TRIM(key),dims(1:4),filename,offset,ref)
    END SELECT

    offset = offset + bytes

  END SUBROUTINE writenode_xdmf


  SUBROUTINE writedataitem(this,dims,filename,offset)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(fileio_xdmf), INTENT(INOUT) :: this
    CHARACTER(LEN=*),   INTENT(IN)    :: dims
    CHARACTER(LEN=*),   INTENT(IN)    :: filename
    INTEGER,            INTENT(IN)    :: offset
    !------------------------------------------------------------------------!
    CHARACTER(LEN=16)                 :: seek
    INTEGER                           :: err
    !------------------------------------------------------------------------!
    WRITE(this%unit, iostat=err)&
           '<DataItem Dimensions=' // trim(dims) // ' ' &
        // 'NumberType="Float" Precision=' // trim(this%realfmt) // lf &
        // 'Format="Binary" Endian=' // trim(this%endian_xdmf)
    WRITE(seek,'(I16)',iostat=err) offset
    WRITE(this%unit, iostat=err)&
       lf // 'Seek="' // trim(adjustl(seek)) // '"'
    WRITE(this%unit, iostat=err)&
           '>' // lf &
        // trim(filename) // lf &
        // '</DataItem>' // lf

  END SUBROUTINE writedataitem

  SUBROUTINE writeattribute(this,Mesh,name,dims,filename,offset,ref)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(fileio_xdmf), INTENT(INOUT) :: this
    CLASS(mesh_base),   INTENT(IN)    :: Mesh
    CHARACTER(LEN=*),   INTENT(IN)    :: name
    CHARACTER(LEN=*),   INTENT(IN)    :: filename
    INTEGER,            INTENT(IN), DIMENSION(:) &
                                      :: dims
    INTEGER,            INTENT(IN)    :: offset
    LOGICAL,            INTENT(IN)    :: ref
    !------------------------------------------------------------------------!
    CHARACTER(LEN=32)                 :: type,center
    INTEGER                           :: dsize, err
    !------------------------------------------------------------------------!
    !data size
    dsize = SIZE(dims)

    ! TODO CHECK THIS - 3D HACK
    ! If the data has mesh sizes, it is cell data
    IF((dims(1).EQ.mesh%INUM).AND.&
       (dims(2).EQ.mesh%JNUM).AND.&
       (dims(3).EQ.mesh%KNUM)) THEN
      center = "Cell"
      dsize = dsize - 3
    ELSE IF((dims(1).EQ.mesh%INUM+1).AND.&
            (dims(2).EQ.mesh%JNUM+1).AND.&
            (dims(3).EQ.mesh%KNUM+1)) THEN
      center = "Node"
      dsize = dsize - 3
    ELSE ! otherwise data for the whole grid
      center = "Grid"
    END IF

    ! The remaining dimensions define the data type
    SELECT CASE(dsize)
    CASE(0)
        type = "Scalar"
    CASE(1)
        type = "Vector"
    CASE DEFAULT
        type = "Matrix"
    END SELECT

    WRITE(this%unit, iostat=err)&
         '<Attribute Name="' // trim(name) // '" '&
      // 'AttributeType="' // trim(type) // '" ' &
      // 'Center="' // trim(center) // '">' // lf

    CALL this%WriteDataItem(getdimsstr(mesh,dims), trim(filename), offset)

    WRITE(this%unit, iostat=err)&
         '</Attribute>' // lf
  END SUBROUTINE writeattribute

  SUBROUTINE writevector(this,Mesh,name,dims,ref1,ref2,ref3,step)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(fileio_xdmf), INTENT(INOUT) :: this
    CLASS(mesh_base),   INTENT(IN)    :: Mesh
    CHARACTER(LEN=*),   INTENT(IN)    :: name
    INTEGER,            INTENT(IN), DIMENSION(3) &
                                      :: dims
    CHARACTER(LEN=*),   INTENT(IN)    :: ref1,ref2,ref3
    CHARACTER(LEN=*),   INTENT(IN)    :: step
    !------------------------------------------------------------------------!
    CHARACTER(LEN=128)                :: dstr
    INTEGER                           :: err
    !------------------------------------------------------------------------!
    dstr = getdimsstr(mesh,dims)
    WRITE(this%unit, iostat=err)&
         '<Attribute Name="' // trim(name) // '" ' &
      // 'AttributeType="Vector" Center="Cell">' // lf&
      // '<DataItem ItemType="Function" Function="JOIN($0, $1, $2)" '&
      // 'Dimensions=' // trim(dstr) // '>' // lf &
      ! first component
      // '<DataItem Reference="/Xdmf/Domain/Grid'&
      // "/Grid[@Name='step" // step // "']" &
      // '/Attribute[@Name=' &
      // "'" // trim(ref1) // "'" // ']/DataItem[1]"/>' // lf &
      ! second component
      // '<DataItem Reference="/Xdmf/Domain/Grid'&
      // "/Grid[@Name='step" // step // "']" &
      // '/Attribute[@Name=' &
      // "'" // trim(ref2) // "'" // ']/DataItem[1]"/>' // lf &
      ! third component
      // '<DataItem Reference="/Xdmf/Domain/Grid'&
      // "/Grid[@Name='step" // step // "']" &
      // '/Attribute[@Name=' &
      // "'" // trim(ref3) // "'" // ']/DataItem[1]"/>' // lf &
      ! end
      // '</DataItem>' // lf &
      // '</Attribute>' // lf
  END SUBROUTINE writevector

  SUBROUTINE writemeshxml(this,Mesh,filename,offset)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(fileio_xdmf), INTENT(INOUT) :: this
    CLASS(mesh_base),   INTENT(IN)    :: Mesh
    CHARACTER(LEN=*),   INTENT(IN)    :: filename
    INTEGER,            INTENT(INOUT) :: offset
    !------------------------------------------------------------------------!
    CHARACTER(LEN=64) :: dims
    CHARACTER(LEN=8)  :: inum, jnum, knum
    INTEGER           :: err
    !------------------------------------------------------------------------!
    WRITE(inum,'(I8)') mesh%INUM+1
    WRITE(jnum,'(I8)') mesh%JNUM+1
    WRITE(knum,'(I8)') mesh%KNUM+1
    WRITE(dims,'(A,A,A,A,A)') trim(adjustl(knum)), ' ',trim(adjustl(jnum)), &
                              ' ',trim(adjustl(inum))

! TODO 3D HACK should not be needed anymore, since now everything is 3D now
! \attention here was a bianglespherical hack
    WRITE(this%unit, iostat=err) &
            '<Topology TopologyType="3DSMesh" ' &
         // 'NumberOfElements="' // trim(dims) // '"/>' // lf &
         // '<Geometry GeometryType="X_Y_Z">' // lf
    WRITE(this%unit, iostat=err) &
         '<DataItem Reference="/Xdmf/Domain/Grid/Grid/Attribute[@Name='//"'/mesh/grid_x'"//']/DataItem[1]"/>' // lf &
      // '<DataItem Reference="/Xdmf/Domain/Grid/Grid/Attribute[@Name='//"'/mesh/grid_y'"//']/DataItem[1]"/>' // lf &
      // '<DataItem Reference="/Xdmf/Domain/Grid/Grid/Attribute[@Name='//"'/mesh/grid_z'"//']/DataItem[1]"/>' // lf

    WRITE(this%unit, iostat=err) &
         '</Geometry>' // lf

  END SUBROUTINE writemeshxml

  SUBROUTINE writexmf(this,Mesh,IO)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(fileio_xdmf), INTENT(INOUT) :: this
    CLASS(mesh_base),   INTENT(IN)    :: Mesh
    TYPE(Dict_TYP),     INTENT(IN), POINTER &
                                      :: IO
    !------------------------------------------------------------------------!
    INTEGER                 :: i, offset, err
    REAL                    :: ftime
    CHARACTER(LEN=4)        :: step
    CHARACTER(LEN=32)       :: time
    CHARACTER(LEN=256)      :: filename
    TYPE(Dict_TYP), POINTER :: meshIO => null()
    !------------------------------------------------------------------------!

    IF (this%GetRank().EQ.0) THEN
      ! write a xmf-file
      OPEN(this%unit, file=trim(this%filename)//'.xmf', &
           status     = 'REPLACE',      &
           action     = 'WRITE',        &
           access     = 'STREAM',       &
           position   = 'REWIND',       &
           iostat     = err)
       IF (err.NE. 0) CALL this%Error("WriteXMF","Can't open xmf-file")

      WRITE(this%unit, iostat=err)&
           '<?xml version="1.0" ?>' // lf &
        // '<!DOCTYPE Xdmf SYSTEM "Xdmf.dtd" []>' // lf &
        // '<Xdmf Version="2.0">' // lf &
        // '<Domain>' // lf &
        // '<Grid Name="mesh" GridType="Collection" CollectionType="Temporal">' // lf

! this works in paraview, but not in visit :(. So we have to use the simpler method for now
!        // '<Time TimeType="List">' // LF &
!        // '<DataItem Format="XML" NumberType="Float" Dimensions="11">' // LF &
!        // '0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0' // LF &
!        // '</DataItem>' // LF &
!        // '</Time>' // LF
!      CALL WriteDataItem_xdmf(this, '"1"', TRIM(filename)//'_'//step//'.bin')
!      WRITE(this%unit, IOSTAT=err)&

      CALL getattr(io,"mesh",meshio)

      DO i=0,this%count
        offset = 13
        ftime = this%stoptime/(this%count)*i
        WRITE(step,'(I4.4)') i
        WRITE(time,'(ES16.9)') ftime
        WRITE(filename, '(A)') trim(this%filename)//"_"//step//".bin"

        WRITE(this%unit, iostat=err)&
             '<Grid Name="step' // step // '" GridType="Uniform">' // lf &
          // '<Time Value="' // trim(time) // '" />' // lf

        CALL this%WriteMeshXML(mesh,filename,offset)
        CALL this%WriteVector(mesh,"/timedisc/velocity", &
                              (/mesh%INUM,mesh%JNUM,mesh%KNUM/), &
                              "/timedisc/xvelocity","/timedisc/yvelocity", &
                              "/timedisc/zvelocity", step)
        CALL this%IterateDict(mesh, io, offset, filename)

        WRITE(this%unit, iostat=err) &
             '</Grid>' // lf
      END DO

      WRITE(this%unit, iostat=err) &
           '</Grid>' // lf


      WRITE(this%unit, iostat=err) &
           '</Domain>' // lf &
        // '</Xdmf>' // lf

      IF(err.NE.0) CALL this%Error("WriteXMF", "Can't write xmf-file")
      CLOSE(this%unit,iostat=err)
      IF(err.NE.0) CALL this%Error("WriteXMF", "Can't close xmf-file")
   END IF

  END SUBROUTINE writexmf


  FUNCTION getdimsstr(Mesh,dims) RESULT(res)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(mesh_base),     INTENT(IN) :: mesh
    INTEGER,DIMENSION(:), INTENT(IN) :: dims
    CHARACTER(LEN=128)               :: res
    !------------------------------------------------------------------------!
    INTEGER                          :: i, l
    CHARACTER(LEN=128), DIMENSION(:), ALLOCATABLE &
                                     :: buf
    !------------------------------------------------------------------------!
    l = SIZE(dims)
    ALLOCATE(buf(l))
    IF(l.GE.3) THEN
      WRITE(buf(1),'(I8)') dims(3)
      WRITE(buf(2),'(I8)') dims(2)
      WRITE(buf(3),'(I8)') dims(1)
      DO i=4,l
        WRITE(buf(i),'(I8)') dims(i)
      END DO
      SELECT CASE(l)
      CASE(2)
        WRITE(res,'(A,A,A,A,A)') '"',trim(adjustl(buf(1))),&
          ' ',trim(adjustl(buf(2))),'"'
      CASE(3)
        WRITE(res,'(A,A,A,A,A,A,A)') '"',trim(adjustl(buf(1))),&
          ' ',trim(adjustl(buf(2))),' ',trim(adjustl(buf(3))),'"'
      CASE(4)
        WRITE(res,'(A,A,A,A,A,A,A,A,A)') '"',trim(adjustl(buf(1))),&
          ' ',trim(adjustl(buf(2))),' ',trim(adjustl(buf(3))), &
          ' ',trim(adjustl(buf(4))),'"'
      CASE(5)
        WRITE(res,'(A,A,A,A,A,A,A,A,A,A,A)') '"',trim(adjustl(buf(1))),&
          ' ',trim(adjustl(buf(2))),' ',trim(adjustl(buf(3))), &
          ' ',trim(adjustl(buf(4))),' ',trim(adjustl(buf(5))),'"'
      END SELECT
    ELSE
      WRITE(buf(1),'(I8)') dims(1)
      WRITE(res,'(A,A,A)') '"',trim(adjustl(buf(1))),'"'
    END IF
    DEALLOCATE(buf)
  END FUNCTION getdimsstr

  SUBROUTINE finalize(this)
    IMPLICIT NONE
    !-----------------------------------------------------------------------!
    CLASS(fileio_xdmf), INTENT(INOUT) :: this
    !-----------------------------------------------------------------------!
    CALL this%fileio_binary%Finalize()
  END SUBROUTINE

END MODULE fileio_xdmf_mod