m_write_silo.f90

 
module m_write_silo
 
  implicit none
  private
 
  include 'silo_f9x.inc'
 
  integer, parameter :: dp       = kind(0.0d0)
  integer, parameter :: line_len = 200
  integer, parameter :: db_type  = db_pdb
 
  public :: silo_mkdir
  public :: silo_create_file
  public :: silo_open_file
  public :: silo_close_file
  public :: silo_set_time_varying
  public :: silo_add_curve
  public :: silo_add_grid
  public :: silo_add_var
  public :: silo_set_mmesh_grid
  public :: silo_set_mmesh_var
 
contains
 
  subroutine silo_create_file(filename, dbix)
    character(len=*), intent(in) :: filename
    integer, intent(out)         :: dbix
    integer                      :: ierr
    character(len=line_len)      :: fileinfo
 
    fileinfo = "A silo file"
    ierr = dbcreate(trim(filename), len_trim(filename), db_clobber, &
         db_local, fileinfo, len_trim(fileinfo), db_type, dbix)
    if (ierr /= 0) then
       print *, "Error creating file ", trim(filename)
       error stop
    end if
  end subroutine silo_create_file
 
  subroutine silo_open_file(filename, dbix)
    character(len=*), intent(in) :: filename
    integer :: dbix, ierr
 
    ierr = dbopen(trim(filename), len_trim(filename), db_type, &
         db_append, dbix)
    if (ierr /= 0) then
       print *, "Error opening file ", trim(filename)
       error stop
    end if
  end subroutine silo_open_file
 
  subroutine silo_close_file(dbix)
    integer, intent(in) :: dbix
    integer :: ierr
 
    ierr = dbclose(dbix)
    if (ierr /= 0) then
       print *, "Error closing file with index", dbix
       error stop
    end if
  end subroutine silo_close_file
 
  subroutine silo_mkdir(dbix, dirname)
    character(len=*), intent(in) :: dirname
    integer, intent(in) :: dbix
    integer :: ierr, iostat
 
    ierr = dbmkdir(dbix, trim(dirname), len_trim(dirname), iostat)
    if (ierr /= 0) then
       print *, "Error creating directory ", dirname
       error stop
    end if
  end subroutine silo_mkdir
 
  subroutine silo_set_time_varying(dbix)
    integer, intent(in)         :: dbix
    integer                     :: ierr
    integer, parameter          :: int_bool(1) = 1
    integer, parameter          :: dims(1) = 1
    character(len=*), parameter :: name1       = "/ConnectivityIsTimeVarying"
    character(len=*), parameter :: name2       = "/MetadataIsTimeVarying"
 
    interface
       integer function dbwrite(dbid, varname, lvarname, var, dims, &
            ndims, datatype)
         use, intrinsic :: iso_c_binding
         integer(c_int) :: dbid, var(*), lvarname, dims(*), ndims, datatype
         character(kind=c_char) :: varname(*)
       end function dbwrite
    end interface
 
    ierr = dbwrite(dbix, name1, len(name1), int_bool, dims, 1, db_int);
    if (ierr /= 0) print *, "Error writing ", trim(name1)
    ierr = dbwrite(dbix, name2, len(name2), int_bool, dims, 1, db_int);
    if (ierr /= 0) print *, "Error writing ", trim(name2)
  end subroutine silo_set_time_varying
 
  subroutine silo_add_curve(dbix, curvename, xvals, yvals)
    character(len=*), intent(in) :: curvename
    integer, intent(in)          :: dbix
    real(dp), intent(in)         :: xvals(:), yvals(:)
    integer                      :: iostat, ierr, dboptix
 
    interface
      integer (c_int) function dbputcurve(dbid, curvename, lcurvename, &
        xvals, yvals, datatype, npoints, optlist_id, status)
        use, intrinsic :: iso_c_binding
        integer(c_int) :: dbid, lcurvename, datatype, npoints, status, optlist_id
        real(c_double) ::  xvals(*), yvals(*)
        character(kind=c_char) :: curvename(*)
      end function dbputcurve
 
      integer (c_int) function dbaddiopt(optlist_id, option, ivalue)
        use, intrinsic :: iso_c_binding
        integer(c_int), intent(in) :: optlist_id, option, ivalue(*)
      end function dbaddiopt
    end interface
 
    if (size(shape(xvals)) /= 1 .or. size(shape(yvals)) /= 1) then
      print *, "Input is not a one-dimensional array. "
    end if
 
    if (size(xvals) /= size(yvals)) then
      print *, "x and y arrays not of the same dimensions"
      print *, size(xvals), " and ", size(yvals)
    end if
 
    ierr = dbmkoptlist(20, dboptix)
    if (ierr /= 0) print *, &
         "Error creating options list in SILO_add_curve ", dboptix
 
    ierr = dbputcurve(dbix, trim(curvename), len_trim(curvename), &
              xvals, yvals, db_double, size(xvals), dboptix, iostat)
 
    if (ierr /= 0) print *, &
      "ERROR: curve object not added to SILO file"
 
    ierr = dbfreeoptlist(dboptix)
    if (ierr /= 0) print *, &
    "Error dbfreeoptlist in SILO_add_curve", ierr
  end subroutine silo_add_curve
 
  subroutine silo_add_grid(dbix, gridname, n_dim, N_r, r_min, dr, &
       lo_offset, hi_offset, n_cycle)
    character(len=*), intent(in) :: gridname
    integer, intent(in)          :: dbix, n_dim, n_r(:)
    integer, intent(in)          :: lo_offset(n_dim), hi_offset(n_dim)
    real(dp), intent(in)         :: r_min(:), dr(:)
    real(dp), allocatable        :: x_coords(:), y_coords(:), z_coords(:)
    integer                      :: i, ierr, iostat, dboptix
    integer, optional, intent(in) :: n_cycle
 
    interface
       function dbputqm(dbid, name, lname, xname, lxname, yname, &
            lyname, zname, lzname, x, y, z, dims, ndims, &
            datatype, coordtype, optlist_id, status)
         use, intrinsic :: iso_c_binding
         integer(c_int) :: dbid, lname, lxname, lyname, lzname, dims(*), ndims
         integer(c_int) :: datatype, coordtype, optlist_id, status, dbputqm
         real(c_double) :: x(*), y(*), z(*)
         character(kind=c_char) :: name(*), xname(*), yname(*), zname(*)
       end function dbputqm
 
       integer (c_int) function dbaddiopt(optlist_id, option, ivalue)
         use, intrinsic :: iso_c_binding
         integer(c_int), intent(in) :: optlist_id, option, ivalue(*)
       end function dbaddiopt
    end interface
 
    if (n_dim < 1 .or. n_dim > 3) then
       print *, "Cannot add grid for which n_dim < 1 or n_dim > 3"
       return
    end if
 
    allocate(x_coords(n_r(1)))
    do i = 1, n_r(1)
       x_coords(i) = r_min(1) + (i-1) * dr(1)
    end do
 
    if (n_dim > 1) then
       allocate(y_coords(n_r(2)))
       do i = 1, n_r(2)
          y_coords(i) = r_min(2) + (i-1) * dr(2)
       end do
    else
       allocate(y_coords(0))
    end if
 
    if (n_dim > 2) then
       allocate(z_coords(n_r(3)))
       do i = 1, n_r(3)
          z_coords(i) = r_min(3) + (i-1) * dr(3)
       end do
    else
       allocate(z_coords(0))
    end if
 
    ! Make option list
    ierr = dbmkoptlist(20, dboptix)
    if (ierr /= 0) print *, &
         "Error creating options list in SILO_add_grid ", dboptix
 
    ! Make iteration number explicit
    ierr = dbaddiopt(dboptix, dbopt_cycle, [n_cycle])
    if (ierr /= 0) print *, &
         "Error passing iteration number SILO_add_grid ", dboptix
 
    ! Set integer options
    ierr = dbaddiopt(dboptix, dbopt_nspace, [n_dim])
    if (ierr /= 0) print *, &
            "Error dbaddiopt in SILO_add_grid: DBOPT_NSPACE", ierr
 
    ierr = dbaddiopt(dboptix, dbopt_lo_offset, lo_offset)
    if (ierr /= 0) print *, &
         "Error dbaddiopt in SILO_add_grid: DBOPT_LO_OFFSET", ierr
 
    ierr = dbaddiopt(dboptix, dbopt_hi_offset, hi_offset)
    if (ierr /= 0) print *, &
         "Error dbaddiopt in SILO_add_grid: DBOPT_HI_OFFSET", ierr
 
    ! Write the grid structure
    ierr = dbputqm(dbix, trim(gridname), len_trim(gridname), &
         'x', 1, 'y', 1, 'z', 1, x_coords, y_coords, z_coords, &
         n_r, n_dim, db_double, db_collinear, dboptix, iostat)
    if (ierr /= 0) print *, &
            "Error dbputqm is SILO_add_grid", ierr
 
    ierr = dbfreeoptlist(dboptix)
    if (ierr /= 0) print *, &
            "Error dbfreeoptlist is SILO_add_grid", ierr
  end subroutine silo_add_grid
 
  subroutine silo_add_var(dbix, dataname, gridname, &
       d_packed, d_shape, n_cycle)
    character(len=*), intent(in) :: gridname, dataname
    real(dp), intent(in)         :: d_packed(:)
    integer, intent(in)          :: dbix, d_shape(:)
 
    integer                      :: dboptix, ierr, iostat
    real(dp)                     :: dummy(1)
    integer, optional, intent(in) :: n_cycle
 
    interface
       function dbputqv1(dbid, name, lname, meshname, lmeshname, &
            var, dims, ndims, mixvar, mixlen, datatype, &
            centering, optlist_id, status)
         use, intrinsic :: iso_c_binding
         integer(c_int) :: dbid, lname, lmeshname, dims(*), ndims, mixlen
         integer(c_int) :: centering, optlist_id, status, datatype, dbputqv1
         real(c_double) :: var(*), mixvar(*)
         character(kind=c_char) :: name(*), meshname(*)
       end function dbputqv1
 
       integer (c_int) function dbaddiopt(optlist_id, option, ivalue)
         use, intrinsic :: iso_c_binding
         integer(c_int), intent(in) :: optlist_id, option, ivalue(*)
       end function dbaddiopt
    end interface
 
    if (size(d_packed) /= product(d_shape)) then
       error stop "Error: d_packed does not correspond to d_shape"
    end if
 
    if (size(d_shape) < 1 .or. size(d_shape) > 3) then
       error stop "Error: size(d_shape) < 1 or size(d_shape) > 3"
    end if
 
    ierr = dbmkoptlist(10, dboptix)
    if (ierr /= 0) then
       error stop "Error creating options list in SILO_add_var"
    end if
 
    ! Make iteration number explicit
    if (present(n_cycle)) then
      ierr = dbaddiopt(dboptix, dbopt_cycle, [n_cycle])
      if (ierr /= 0) print *, &
           "Error passing iteration number SILO_add_var ", dboptix
    end if
 
    ierr = dbaddiopt(dboptix, dbopt_hide_from_gui, [1])
 
    ! Write the data to the grid
    ierr = dbputqv1(dbix, trim(dataname), len_trim(dataname), &
         trim(gridname), len_trim(gridname), d_packed, d_shape, &
         size(d_shape), dummy, 0, db_double, db_zonecent, dboptix, iostat)
    if (ierr /= 0) then
       print *, "Error dbputqv1 in SILO_add_var", ierr
       error stop
    end if
 
    ierr = dbfreeoptlist(dboptix)
  end subroutine silo_add_var
 
  subroutine silo_set_mmesh_grid(dbix, mmname, gridnames, n_cycle, time)
    character(len=*), intent(in)   :: mmname, gridnames(:)
    integer, intent(in)            :: dbix
    integer, intent(in), optional  :: n_cycle
    real(dp), intent(in), optional :: time
 
    integer                        :: i, ierr,length
    integer                        :: dboptix, iostat, old_str_len
    integer                        :: n_grids, name_len, total_len
    integer, allocatable           :: m_types(:), name_lengths(:)
    character(len=:), allocatable  :: mnames
 
    interface
       function dbputmmesh(dbid, name, lname, nmesh, meshnames, lmeshnames, &
            meshtypes, optlist_id, status)
         use, intrinsic :: iso_c_binding
         integer(c_int) :: dbputmmesh, lmeshnames(*)
         integer(c_int) :: dbid, lname, nmesh, meshtypes(*), optlist_id, status
         character(kind=c_char) :: name(*), meshnames(*)
       end function dbputmmesh
    end interface
 
    n_grids = size(gridnames)
    if (n_grids < 1) then
       error stop "SILO_set_mmesh_grid: error too few grids (<1)"
    end if
 
    name_len  = len(gridnames(1))
    total_len = name_len * n_grids
    allocate(character(total_len) :: mnames)
    allocate(name_lengths(n_grids))
    allocate(m_types(n_grids))
 
    do i = 1, n_grids
       mnames((i-1)*name_len+1:i*name_len) = trim(gridnames(i)) // char(0)
    end do
 
    old_str_len  = dbset2dstrlen(name_len)
    m_types      = db_quad_rect
    name_lengths = name_len
 
    ierr = dbmkoptlist(10, dboptix)
 
    if (present(n_cycle)) then
       ierr = dbaddiopt(dboptix, dbopt_cycle, n_cycle)
    end if
 
    if (present(time)) then
       ierr = dbadddopt(dboptix, dbopt_dtime, time)
    end if
 
    ierr = dbputmmesh(dbix, trim(mmname), len_trim(mmname), n_grids, &
         mnames(1:total_len), name_lengths, m_types, dboptix, iostat)
    if (ierr /= 0) then
       error stop "Error calling dbputmmesh"
    end if
 
    ierr = dbfreeoptlist(dboptix)
    length = dbset2dstrlen(old_str_len)
  end subroutine silo_set_mmesh_grid
 
  subroutine silo_set_mmesh_var(dbix, mvname, mmname, datanames, n_cycle, time)
    character(len=*), intent(in)   :: mvname, mmname, datanames(:)
    integer, intent(in)            :: dbix
    integer, intent(in), optional  :: n_cycle
    real(dp), intent(in), optional :: time
 
 
    integer                        :: i, ierr, dboptix, iostat,length
    integer                        :: old_str_len, n_grids, name_len, total_len
    integer, allocatable           :: m_types(:), name_lengths(:)
    character(:), allocatable      :: dnames
 
    interface
       function dbputmvar(dbid, name, lname, nlevels, meshnames, &
            lmnames, meshtypes, optlist_id, status)
         use, intrinsic :: iso_c_binding
         integer(c_int) :: dbputmvar, lmnames(*)
         integer(c_int) :: dbid, lname, nlevels, meshtypes(*)
         integer(c_int) :: optlist_id, status
         character(kind=c_char) :: name(*), meshnames(*)
       end function dbputmvar
    end interface
 
    n_grids = size(datanames)
    if (n_grids < 1) then
       error stop "SILO_set_mmesh_var: error too few grids (<1)"
    end if
 
    name_len = len(datanames(1))
    total_len = name_len * n_grids
    allocate(character(total_len) :: dnames)
    allocate(name_lengths(n_grids))
    allocate(m_types(n_grids))
 
    do i = 1, n_grids
       dnames((i-1)*name_len+1:i*name_len) = trim(datanames(i)) // char(0)
    end do
    old_str_len  = dbset2dstrlen(name_len)
    m_types      = db_quadvar
    name_lengths = name_len
 
    ierr = dbmkoptlist(10, dboptix)
    if (ierr /= 0) then
       error stop "Error creating options list in SILO_set_mmesh_var"
    end if
 
    if (present(n_cycle)) then
       ierr = dbaddiopt(dboptix, dbopt_cycle, n_cycle)
    end if
 
    if (present(time)) then
       ierr = dbadddopt(dboptix, dbopt_dtime, time)
    end if
 
    ! Jannis: valgrind complains about a memory leak due to this call, but if
    ! that leak is there, it seems to be Silo's fault.
    ierr = dbaddcopt(dboptix, dbopt_mmesh_name, &
         trim(mmname), len_trim(mmname))
    if (ierr /= 0) print *, &
            "Error dbaddiopt in SILO_set_mmesh_var: DBOPT_MMESH_NAME", ierr
 
    ierr = dbputmvar(dbix, trim(mvname), len_trim(mvname), n_grids, &
         dnames(1:total_len), name_lengths, m_types, dboptix, iostat)
    if (ierr /= 0) then
       error stop "Error calling dbputmvar"
    end if
 
    ierr = dbfreeoptlist(dboptix)
    length = dbset2dstrlen(old_str_len)
  end subroutine silo_set_mmesh_var
 
  ! !> Add a curve-object (pairs of x-y values) to the Silo file
  ! subroutine SILO_add_curve(dbix, curvename, xvals, yvals, xname, yname)
  !   character(len=*), intent(in) :: curvename
  !   integer, intent(in)          :: dbix
  !   real(dp), intent(in)         :: xvals(:), yvals(:)
  !   character(len=*), intent(in) :: xname, yname
  !   integer                      :: iostat, ierr, dboptix
  !
  !   interface
  !      integer (c_int) function dbputcurve(dbid, curvename, lcurvename, &
  !           xvals, yvals, datatype, npoints, optlist_id, status)
  !        use, intrinsic :: iso_c_binding
  !        integer(c_int) :: dbid, lcurvename, datatype, npoints, status, optlist_id
  !        real(c_double) ::  xvals(*), yvals(*)
  !        character(kind=c_char) :: curvename(*)
  !      end function dbputcurve
  !
  !      integer (c_int) function dbaddiopt(optlist_id, option, ivalue)
  !        use, intrinsic :: iso_c_binding
  !        integer(c_int), intent(in) :: optlist_id, option, ivalue(*)
  !      end function dbaddiopt
  !   end interface
  !
  !   if (size(xvals) /= size(yvals)) &
  !      error stop "SILO_add_curve: x and y arrays unequal size"
  !
  !   ierr = dbmkoptlist(20, dboptix)
  !   if (ierr /= 0) error stop "creating options list in SILO_add_curve "
  !
  !   ierr = dbaddcopt(dboptix, DBOPT_XLABEL, trim(xname), len_trim(xname))
  !   if (ierr /= 0) error stop "adding option in SILO_add_curve"
  !   ierr = dbaddcopt(dboptix, DBOPT_YLABEL, trim(yname), len_trim(yname))
  !   if (ierr /= 0) error stop "adding option in SILO_add_curve"
  !
  !   ierr = dbputcurve(dbix, trim(curvename), len_trim(curvename), &
  !        xvals, yvals, DB_DOUBLE, size(xvals), dboptix, iostat)
  !   if (ierr /= 0) error stop "curve object not added to SILO file"
  !
  !   ierr = dbfreeoptlist(dboptix)
  !   if (ierr /= 0) error stop "dbfreeoptlist in SILO_add_curve"
  ! end subroutine SILO_add_curve
 
end module m_write_silo