vtkhdf_file.F90

Go to the documentation of this file.

! Copyright (c) 2026, The Neko Authors
! All rights reserved.
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions
! are met:
!
!   * Redistributions of source code must retain the above copyright
!     notice, this list of conditions and the following disclaimer.
!
!   * Redistributions in binary form must reproduce the above
!     copyright notice, this list of conditions and the following
!     disclaimer in the documentation and/or other materials provided
!     with the distribution.
!
!   * Neither the name of the authors nor the names of its
!     contributors may be used to endorse or promote products derived
!     from this software without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
! FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
! COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
! INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
! BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
! ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
! POSSIBILITY OF SUCH DAMAGE.
!
module vtkhdf_file
  use num_types, only : rp, sp, dp, qp, i8
  use generic_file, only : generic_file_t
  use checkpoint, only : chkp_t
  use utils, only : neko_error, neko_warning, filename_split, &
       nonlinear_index, linear_index, mkdir
  use mesh, only : mesh_t
  use field, only : field_t
  use field_list, only : field_list_t
  use field_series, only : field_series_t, field_series_ptr_t
  use dofmap, only : dofmap_t
  use logger, only : neko_log
  use comm, only : pe_rank, pe_size, neko_comm
  use device, only : device_to_host, host_to_device
  use mpi_f08, only : mpi_info_null, mpi_allreduce, mpi_allgather, &
       mpi_in_place, mpi_integer, mpi_sum, mpi_max, mpi_comm_size, mpi_exscan, &
       mpi_barrier, mpi_bcast, mpi_logical
  use vtk, only : vtk_ordering
#ifdef HAVE_HDF5
  use hdf5, only : &
       hid_t, hsize_t, size_t, &
       h5open_f, h5close_f, &
       h5fcreate_f, h5fopen_f, h5fclose_f, h5fflush_f, h5fget_obj_count_f, &
       h5f_obj_all_f, h5f_scope_global_f, &
       h5gcreate_f, h5gopen_f, h5gclose_f, &
       h5acreate_f, h5aopen_f, h5awrite_f, h5aread_f, h5aclose_f, h5aexists_f, &
       h5aopen_by_name_f, &
       h5dcreate_f, h5dopen_f, h5dwrite_f, h5dread_f, h5dclose_f, &
       h5dget_create_plist_f, &
       h5tcopy_f, h5tclose_f, h5tset_strpad_f, h5tset_size_f, &
       h5screate_f, h5screate_simple_f, h5sclose_f, &
       h5sselect_hyperslab_f, h5sselect_all_f, h5sget_simple_extent_dims_f, &
       h5dget_space_f, h5dset_extent_f, &
       h5pcreate_f, h5pclose_f, h5pset_fapl_mpio_f, h5pset_dxpl_mpio_f, &
       h5pget_virtual_count_f, h5pget_virtual_filename_f, &
       h5lexists_f, h5ldelete_f, &
       h5p_file_access_f, h5p_dataset_xfer_f, h5p_dataset_create_f, &
       h5f_acc_trunc_f, h5f_acc_rdwr_f, h5f_acc_rdonly_f, h5pset_chunk_f, &
       h5t_std_u8le, h5t_native_integer, h5t_fortran_s1, h5t_str_nullterm_f, &
       h5kind_to_type, h5_real_kind, h5_integer_kind, h5pset_virtual_f, &
       h5s_scalar_f, h5s_select_set_f, h5fd_mpio_collective_f, h5p_default_f, &
       h5s_unlimited_f
#endif
  implicit none
  private
 
  type, public, extends(generic_file_t) :: vtkhdf_file_t
     logical :: amr_enabled = .false.
     logical :: subdivide = .false.
     logical :: enable_vds = .false.
     integer :: precision = -1
   contains
     procedure :: get_vtkhdf_fname => vtkhdf_file_get_fname
     procedure :: get_next_output_fname => vtkhdf_file_get_next_output_fname
     procedure :: read => vtkhdf_file_read
     procedure :: write => vtkhdf_file_write
     procedure :: set_overwrite => vtkhdf_file_set_overwrite
     procedure :: enable_amr => vtkhdf_file_enable_amr
     procedure :: set_precision => vtkhdf_file_set_precision
     procedure :: set_subdivide => vtkhdf_file_set_subdivide
  end type vtkhdf_file_t
 
  integer, dimension(2), parameter :: vtkhdf_version = [2, 6]
 
contains
 
  ! -------------------------------------------------------------------------- !
  ! Well defined subroutines
 
  subroutine vtkhdf_file_set_overwrite(this, overwrite)
    class(vtkhdf_file_t), intent(inout) :: this
    logical, intent(in) :: overwrite
    this%overwrite = overwrite
  end subroutine vtkhdf_file_set_overwrite
 
  subroutine vtkhdf_file_enable_amr(this)
    class(vtkhdf_file_t), intent(inout) :: this
    this%amr_enabled = .false.
  end subroutine vtkhdf_file_enable_amr
 
  subroutine vtkhdf_file_set_precision(this, precision)
    class(vtkhdf_file_t), intent(inout) :: this
    integer, intent(in) :: precision
    this%precision = precision
  end subroutine vtkhdf_file_set_precision
 
  function vtkhdf_file_get_fname(this) result(base_fname)
    class(vtkhdf_file_t), intent(in) :: this
    character(len=1024) :: base_fname
    character(len=1024) :: fname
    character(len=1024) :: path, name, suffix
 
    fname = trim(this%get_base_fname())
    call filename_split(fname, path, name, suffix)
 
    write(base_fname, '(A,A,"_",I0,A)') &
         trim(path), trim(name), this%get_start_counter(), trim(suffix)
 
  end function vtkhdf_file_get_fname
 
  function vtkhdf_file_get_next_output_fname(this) result(fname)
    class(vtkhdf_file_t), intent(in) :: this
    character(len=1024) :: fname
 
    fname = this%get_vtkhdf_fname()
 
  end function vtkhdf_file_get_next_output_fname
 
  subroutine vtkhdf_file_set_subdivide(this, subdivide)
    class(vtkhdf_file_t), intent(inout) :: this
    logical, intent(in) :: subdivide
    this%subdivide = subdivide
  end subroutine vtkhdf_file_set_subdivide
 
#ifdef HAVE_HDF5
  ! -------------------------------------------------------------------------- !
  ! HDF5 Required subroutines
 
  subroutine vtkhdf_file_write(this, data, t)
    class(vtkhdf_file_t), intent(inout) :: this
    class(*), target, intent(in) :: data
    real(kind=rp), intent(in), optional :: t
    type(mesh_t), pointer :: msh
    type(dofmap_t), pointer :: dof
    type(field_list_t) :: fields
    integer :: ierr, mpi_info, mpi_comm, i, n_fields
    integer(hid_t) :: plist_id, file_id, attr_id, vtkhdf_grp
    integer(hid_t) :: filespace, H5T_NEKO_STRING
    integer(hsize_t), dimension(1) :: vdims
    integer(size_t) :: type_len
    integer :: lx, ly, lz
    integer :: local_points, local_cells, local_conn
    integer :: total_points, total_cells, total_conn
    integer :: point_offset
    integer :: max_local_points
    integer, allocatable :: part_points(:), part_cells(:), part_conns(:)
    character(len=1024) :: fname
    character(len=16) :: type_str
    logical :: exists
    integer :: counter
 
    ! Determine mesh and field data
    select type(data)
    type is (field_t)
       msh => data%msh
       dof => data%dof
       n_fields = 1
       call fields%init(1)
       ! Call the specific bindings; resolving the assign generic here
       ! requires matching a TARGET actual to a POINTER, INTENT(IN) dummy
       ! (F2008), which the Fujitsu compiler rejects
       call fields%assign_to_field(1, data)
    type is (field_list_t)
       msh => data%msh(1)
       dof => data%dof(1)
       call fields%assign_to_list(data)
    class default
       call neko_error('Invalid data type for vtkhdf_file_write')
    end select
 
    ! Check conditions to ensure the input data is supported.
    if (.not. associated(msh)) then
       call neko_error('Mesh must be associated for vtkhdf_file_write')
    end if
    if (dof%Xh%lx .lt. 2 .or. dof%Xh%ly .lt. 2) then
       call neko_error('VTKHDF linear output requires lx, ly >= 2')
    end if
    if (msh%gdim .eq. 3 .and. dof%Xh%lz .lt. 2) then
       call neko_error('VTKHDF linear output requires lz >= 2 in 3D')
    end if
    if (msh%gdim .lt. 2 .or. msh%gdim .gt. 3) then
       call neko_error('VTKHDF output only supports 2D and 3D meshes')
    end if
 
    ! Ensure precision is set and are valid.
    if (this%precision .gt. rp) then
       this%precision = rp
       call neko_warning('Requested precision is higher than working precision')
    else if (this%precision .eq. -1) then
       this%precision = rp
    end if
 
    call this%increment_counter()
    fname = trim(this%get_vtkhdf_fname())
    counter = this%get_counter() - this%get_start_counter()
 
    mpi_info = mpi_info_null%mpi_val
    mpi_comm = neko_comm%mpi_val
 
    call h5open_f(ierr)
    call h5pcreate_f(h5p_file_access_f, plist_id, ierr)
    call h5pset_fapl_mpio_f(plist_id, mpi_comm, mpi_info, ierr)
 
    if (counter .eq. 0) then
       ! First write: always create a fresh file to avoid stale data
       call h5fcreate_f(fname, h5f_acc_trunc_f, &
            file_id, ierr, access_prp = plist_id)
    else
       call h5fopen_f(fname, h5f_acc_rdwr_f, file_id, ierr, &
            access_prp = plist_id)
    end if
 
    ! Create/open VTKHDF root group with vtkhdf_version and type attributes
    call h5lexists_f(file_id, "VTKHDF", exists, ierr)
    if (exists) then
       call h5gopen_f(file_id, "VTKHDF", vtkhdf_grp, ierr)
    else
       call h5gcreate_f(file_id, "VTKHDF", vtkhdf_grp, ierr)
 
       ! Write Version attribute
       vdims = 2_hsize_t
       call h5screate_simple_f(1, vdims, filespace, ierr)
       call h5acreate_f(vtkhdf_grp, "Version", h5t_native_integer, filespace, &
            attr_id, ierr)
       call h5awrite_f(attr_id, h5t_native_integer, vtkhdf_version, vdims, ierr)
       call h5aclose_f(attr_id, ierr)
       call h5sclose_f(filespace, ierr)
 
       ! Write Type attribute "UnstructuredGrid" as a fixed-length string
       type_str = "UnstructuredGrid"
       type_len = int(len_trim(type_str), kind=size_t)
       vdims = 1_hsize_t
       call h5screate_f(h5s_scalar_f, filespace, ierr)
 
       call h5tcopy_f(h5t_fortran_s1, h5t_neko_string, ierr)
       call h5tset_size_f(h5t_neko_string, type_len, ierr)
       call h5tset_strpad_f(h5t_neko_string, h5t_str_nullterm_f, ierr)
 
       call h5acreate_f(vtkhdf_grp, "Type", h5t_neko_string, filespace, &
            attr_id, ierr)
       call h5awrite_f(attr_id, h5t_neko_string, [type_str], vdims, ierr)
       call h5aclose_f(attr_id, ierr)
 
       call h5tclose_f(h5t_neko_string, ierr)
       call h5sclose_f(filespace, ierr)
    end if
 
    if (present(t)) then
       call vtkhdf_write_steps(vtkhdf_grp, counter, t)
    end if
 
    if (associated(msh)) then
       call vtkhdf_write_mesh(vtkhdf_grp, dof, msh, &
            this%amr_enabled, counter, this%subdivide, t)
    end if
 
    ! Write field data in PointData group
    if (fields%size() .gt. 0) then
       call vtkhdf_write_pointdata(vtkhdf_grp, fields, this%precision, &
            counter, fname, t)
    end if
 
    call h5gclose_f(vtkhdf_grp, ierr)
    call h5pclose_f(plist_id, ierr)
 
    ! A leaked id defers the real (weak) close indefinitely, so the
    ! metadata cache would never reach disk; flush explicitly and warn
    ! if anything besides the file itself is still open.
    block
      integer(size_t) :: obj_count
      character(len=80) :: wrn_buf
      call h5fget_obj_count_f(file_id, h5f_obj_all_f, obj_count, ierr)
      if (obj_count .gt. 1_size_t .and. pe_rank .eq. 0) then
         write(wrn_buf, '(A,I0,A)') 'VTKHDF: ', obj_count - 1, &
              ' HDF5 id(s) still open at file close'
         call neko_warning(trim(wrn_buf))
      end if
    end block
    call h5fflush_f(file_id, h5f_scope_global_f, ierr)
    call h5fclose_f(file_id, ierr)
    call h5close_f(ierr)
 
    call fields%free()
 
  end subroutine vtkhdf_file_write
 
  ! -------------------------------------------------------------------------- !
  ! Internal helper subroutines for VTKHDF writing
 
  subroutine vtkhdf_write_mesh(vtkhdf_grp, dof, msh, amr, counter, subdivide, t)
    type(dofmap_t), intent(in) :: dof
    type(mesh_t), intent(in) :: msh
    integer(hid_t), intent(in) :: vtkhdf_grp
    logical, intent(in) :: amr
    integer, intent(in) :: counter
    logical, intent(in) :: subdivide
    real(kind=rp), intent(in), optional :: t
 
    integer(kind=1) :: VTK_cell_type
    integer :: ierr, i, ii, jj, kk, el, local_idx
    integer :: lx, ly, lz, npts_per_cell, nodes_per_cell, cells_per_element
    integer :: local_points, local_cells, local_conn
    integer :: total_points, total_cells, total_conn
    integer :: point_offset, max_local_points
    integer :: total_offsets, cell_offset, conn_offset, offsets_offset
    integer :: max_local_cells, max_local_conn
    integer(hid_t) :: xf_id, dset_id, dcpl_id, grp_id, attr_id
    integer(hid_t) :: filespace, memspace, H5T_NEKO_DOUBLE
    integer(hsize_t), dimension(1) :: dcount, vdims, maxdims, doffset, chunkdims
    integer(hsize_t), dimension(2) :: dcount2, vdims2, maxdims2, doffset2
    integer(kind=i8) :: i8_value
    logical :: exists
    integer, dimension(3) :: component_sizes
    integer, dimension(3) :: component_offsets
    integer, dimension(3) :: component_max_sizes
 
    lx = dof%Xh%lx
    ly = dof%Xh%ly
    lz = dof%Xh%lz
 
    if (subdivide .and. msh%gdim .eq. 3) then
       vtk_cell_type = int(12, kind=1) ! VTK_HEXAHEDRON
       cells_per_element = (lx - 1) * (ly - 1) * (lz - 1)
       nodes_per_cell = 8
    else if (subdivide .and. msh%gdim .eq. 2) then
       vtk_cell_type = int(9, kind=1) ! VTK_QUAD
       cells_per_element = (lx - 1) * (ly - 1)
       nodes_per_cell = 4
    else if (msh%gdim .eq. 3) then
       vtk_cell_type = int(72, kind=1) ! VTK_LAGRANGE_HEXAHEDRON
       cells_per_element = 1
       nodes_per_cell = lx * ly * lz
    else if (msh%gdim .eq. 2) then
       vtk_cell_type = int(70, kind=1) ! VTK_LAGRANGE_QUADRILATERAL
       cells_per_element = 1
       nodes_per_cell = lx * ly
    end if
 
    ! --- Build the number of cells and the connectivity
    local_points = dof%size()
    local_cells = msh%nelv * cells_per_element
    local_conn = local_cells * nodes_per_cell
 
    total_points = dof%global_size()
    total_cells = msh%glb_nelv * cells_per_element
    total_conn = total_cells * nodes_per_cell
 
    component_sizes = [local_points, local_cells, local_conn]
    component_offsets = 0
    component_max_sizes = 0
 
    call mpi_exscan(component_sizes, component_offsets, 3, mpi_integer, &
         mpi_sum, neko_comm, ierr)
    call mpi_allreduce(component_sizes, component_max_sizes, 3, mpi_integer, &
         mpi_max, neko_comm, ierr)
 
    point_offset = component_offsets(1)
    cell_offset = component_offsets(2)
    conn_offset = component_offsets(3)
    max_local_points = component_max_sizes(1)
    max_local_cells = component_max_sizes(2)
    max_local_conn = component_max_sizes(3)
 
    offsets_offset = cell_offset + pe_rank
    total_offsets = total_cells + pe_size
 
    ! Create collective transfer property list
    call h5pcreate_f(h5p_dataset_xfer_f, xf_id, ierr)
    call h5pset_dxpl_mpio_f(xf_id, h5fd_mpio_collective_f, ierr)
 
    ! --- NumberOfPoints, NumberOfCells, NumberOfConnectivityIds ---
    ! These datasets must accumulate nPieces entries per timestep,
    ! giving a total size of nSteps * nPieces. VTK's reader computes
    ! numberOfPieces = dims[0] / nSteps, so missing entries cause
    ! garbage reads.
    block
      integer(hsize_t), dimension(1) :: nof_dims, nof_maxdims
      integer(hsize_t), dimension(1) :: nof_count, nof_offset, nof_chunk
      integer(hid_t) :: nof_filespace, nof_memspace, nof_dcpl
      integer(kind=i8) :: nof_values(3)
 
      nof_count(1) = 1_hsize_t
      nof_offset(1) = int(counter, hsize_t) * int(pe_size, hsize_t) &
           + int(pe_rank, hsize_t)
      nof_chunk(1) = max(1_hsize_t, int(pe_size, hsize_t))
      nof_values = [int(local_points, kind=i8), int(local_cells, kind=i8), &
           int(local_conn, kind=i8)]
 
      call h5pcreate_f(h5p_dataset_create_f, nof_dcpl, ierr)
      call h5pset_chunk_f(nof_dcpl, 1, nof_chunk, ierr)
 
      call vtkhdf_write_numberof(vtkhdf_grp, "NumberOfPoints", &
           nof_values(1), nof_offset, nof_count, nof_dcpl, &
           counter, xf_id, ierr)
      call vtkhdf_write_numberof(vtkhdf_grp, "NumberOfCells", &
           nof_values(2), nof_offset, nof_count, nof_dcpl, &
           counter, xf_id, ierr)
      call vtkhdf_write_numberof(vtkhdf_grp, "NumberOfConnectivityIds", &
           nof_values(3), nof_offset, nof_count, nof_dcpl, &
           counter, xf_id, ierr)
 
      call h5pclose_f(nof_dcpl, ierr)
    end block
 
    ! --- Points dataset (global coordinates) ---
    call h5lexists_f(vtkhdf_grp, "Points", exists, ierr)
    if (.not. exists) then
 
       vdims2 = [3_hsize_t, int(total_points, hsize_t)]
       maxdims2 = [3_hsize_t, h5s_unlimited_f]
       chunkdims(1) = int(max(1, min(max_local_points, total_points)), hsize_t)
       dcount2 = [3_hsize_t, int(local_points, hsize_t)]
       doffset2 = [0_hsize_t, int(point_offset, hsize_t)]
       h5t_neko_double = h5kind_to_type(dp, h5_real_kind)
 
       call h5pcreate_f(h5p_dataset_create_f, dcpl_id, ierr)
       call h5screate_simple_f(2, dcount2, memspace, ierr)
       call h5screate_simple_f(2, vdims2, filespace, ierr, maxdims2)
 
       call h5pset_chunk_f(dcpl_id, 2, [3_hsize_t, chunkdims(1)], ierr)
       call h5dcreate_f(vtkhdf_grp, "Points", h5t_neko_double, &
            filespace, dset_id, ierr, dcpl_id = dcpl_id)
       call h5sselect_hyperslab_f(filespace, h5s_select_set_f, &
            doffset2, dcount2, ierr)
 
       block
         real(kind=dp), allocatable :: coords(:,:)
 
         allocate(coords(3, local_points))
         do concurrent(local_idx = 1:local_points)
            block
              integer :: idx(4)
              real(kind=dp) :: x, y, z
              idx = nonlinear_index(local_idx, lx, ly, lz)
 
              x = real(dof%x(idx(1), idx(2), idx(3), idx(4)), dp)
              y = real(dof%y(idx(1), idx(2), idx(3), idx(4)), dp)
              z = real(dof%z(idx(1), idx(2), idx(3), idx(4)), dp)
 
              coords(1, local_idx) = x
              coords(2, local_idx) = y
              coords(3, local_idx) = z
            end block
         end do
         call h5dwrite_f(dset_id, h5t_neko_double, coords, dcount2, ierr, &
              file_space_id = filespace, mem_space_id = memspace, &
              xfer_prp = xf_id)
         deallocate(coords)
       end block
 
       call h5dclose_f(dset_id, ierr)
       call h5sclose_f(filespace, ierr)
       call h5sclose_f(memspace, ierr)
       call h5pclose_f(dcpl_id, ierr)
    end if
 
    ! --- Connectivity dataset ---
    call h5lexists_f(vtkhdf_grp, "Connectivity", exists, ierr)
    if (exists) call h5ldelete_f(vtkhdf_grp, "Connectivity", ierr)
 
    vdims = int(total_conn, hsize_t)
    maxdims = h5s_unlimited_f
    chunkdims = int(max(1, min(max_local_conn, total_conn)), hsize_t)
    dcount = int(local_conn, hsize_t)
    doffset = int(conn_offset, hsize_t)
 
    call h5pcreate_f(h5p_dataset_create_f, dcpl_id, ierr)
    call h5screate_simple_f(1, dcount, memspace, ierr)
    call h5screate_simple_f(1, vdims, filespace, ierr, maxdims)
 
    call h5pset_chunk_f(dcpl_id, 1, chunkdims, ierr)
    call h5dcreate_f(vtkhdf_grp, "Connectivity", h5t_native_integer, &
         filespace, dset_id, ierr, dcpl_id = dcpl_id)
    call h5sselect_hyperslab_f(filespace, h5s_select_set_f, &
         doffset, dcount, ierr)
 
    block
      integer, allocatable :: connectivity(:)
 
      allocate(connectivity(local_conn))
      call vtkhdf_build_connectivity(connectivity, vtk_cell_type, msh, dof, &
           subdivide)
      call h5dwrite_f(dset_id, h5t_native_integer, connectivity, dcount, &
           ierr, file_space_id = filespace, mem_space_id = memspace, &
           xfer_prp = xf_id)
      deallocate(connectivity)
    end block
 
    call h5dclose_f(dset_id, ierr)
    call h5sclose_f(filespace, ierr)
    call h5sclose_f(memspace, ierr)
    call h5pclose_f(dcpl_id, ierr)
 
    ! --- Offsets dataset ---
    call h5lexists_f(vtkhdf_grp, "Offsets", exists, ierr)
    if (exists) call h5ldelete_f(vtkhdf_grp, "Offsets", ierr)
 
    vdims = int(total_offsets, hsize_t)
    maxdims = h5s_unlimited_f
    chunkdims = int(max(1, min(max_local_cells + 1, total_offsets)), hsize_t)
    dcount = int(local_cells + 1, hsize_t)
    doffset = int(offsets_offset, hsize_t)
 
    call h5pcreate_f(h5p_dataset_create_f, dcpl_id, ierr)
    call h5screate_simple_f(1, dcount, memspace, ierr)
    call h5screate_simple_f(1, vdims, filespace, ierr, maxdims)
 
    call h5pset_chunk_f(dcpl_id, 1, chunkdims, ierr)
    call h5dcreate_f(vtkhdf_grp, "Offsets", h5t_native_integer, &
         filespace, dset_id, ierr, dcpl_id = dcpl_id)
    call h5sselect_hyperslab_f(filespace, h5s_select_set_f, &
         doffset, dcount, ierr)
 
    block
      integer, allocatable :: offsets(:)
 
      allocate(offsets(local_cells + 1))
      do concurrent(i = 1:local_cells)
         offsets(i) = (i - 1) * nodes_per_cell
      end do
      offsets(local_cells + 1) = local_conn
      call h5dwrite_f(dset_id, h5t_native_integer, offsets, dcount, ierr, &
           file_space_id = filespace, mem_space_id = memspace, xfer_prp = xf_id)
      deallocate(offsets)
    end block
 
    call h5dclose_f(dset_id, ierr)
    call h5sclose_f(filespace, ierr)
    call h5sclose_f(memspace, ierr)
    call h5pclose_f(dcpl_id, ierr)
 
    ! --- Types dataset (VTK cell types) ---
    call h5lexists_f(vtkhdf_grp, "Types", exists, ierr)
    if (exists) call h5ldelete_f(vtkhdf_grp, "Types", ierr)
 
    vdims = int(total_cells, hsize_t)
    maxdims = h5s_unlimited_f
    chunkdims = int(max(1, min(max_local_cells, total_cells)), hsize_t)
    dcount = int(local_cells, hsize_t)
    doffset = int(cell_offset, hsize_t)
 
    call h5pcreate_f(h5p_dataset_create_f, dcpl_id, ierr)
    call h5screate_simple_f(1, dcount, memspace, ierr)
    call h5screate_simple_f(1, vdims, filespace, ierr, maxdims)
 
    call h5pset_chunk_f(dcpl_id, 1, chunkdims, ierr)
    call h5dcreate_f(vtkhdf_grp, "Types", h5t_std_u8le, &
         filespace, dset_id, ierr, dcpl_id = dcpl_id)
    call h5sselect_hyperslab_f(filespace, h5s_select_set_f, &
         doffset, dcount, ierr)
 
    block
      integer(kind=1), allocatable :: cell_types(:)
      allocate(cell_types(local_cells), source=vtk_cell_type)
      call h5dwrite_f(dset_id, h5t_std_u8le, cell_types, dcount, ierr, &
           file_space_id = filespace, mem_space_id = memspace, xfer_prp = xf_id)
      deallocate(cell_types)
    end block
 
    call h5dclose_f(dset_id, ierr)
    call h5sclose_f(filespace, ierr)
    call h5sclose_f(memspace, ierr)
    call h5pclose_f(dcpl_id, ierr)
 
    if (present(t)) then
       ! Open Steps group
       call h5gopen_f(vtkhdf_grp, "Steps", grp_id, ierr)
 
       ! --- NumberOfParts ---
       call vtkhdf_write_i8_at(grp_id, "NumberOfParts", int(pe_size, kind=i8), &
            counter)
 
       ! --- PartOffsets ---
       i8_value = 0_i8
       if (amr) then
          i8_value = int(counter - 1, kind=i8) * int(pe_size, kind=i8)
          call vtkhdf_write_i8_at(grp_id, "PartOffsets", i8_value, counter)
 
          i8_value = int(counter - 1, kind=i8) * int(total_points, kind=i8)
          call vtkhdf_write_i8_at(grp_id, "PointOffsets", i8_value, counter)
 
          i8_value = int(counter - 1, kind=i8) * int(total_cells, kind=i8)
          call vtkhdf_write_i8_at(grp_id, "CellOffsets", i8_value, counter)
 
          i8_value = int(counter - 1, kind=i8) * int(total_conn, kind=i8)
          call vtkhdf_write_i8_at(grp_id, "ConnectivityIdOffsets", i8_value, &
               counter)
 
       else
          i8_value = 0_i8
          call vtkhdf_write_i8_at(grp_id, "PartOffsets", i8_value, counter)
          call vtkhdf_write_i8_at(grp_id, "PointOffsets", i8_value, counter)
          call vtkhdf_write_i8_at(grp_id, "CellOffsets", i8_value, counter)
          call vtkhdf_write_i8_at(grp_id, "ConnectivityIdOffsets", i8_value, &
               counter)
       end if
 
       call h5gclose_f(grp_id, ierr)
    end if
 
    call h5pclose_f(xf_id, ierr)
 
  end subroutine vtkhdf_write_mesh
 
  subroutine vtkhdf_write_steps(vtkhdf_grp, counter, t)
    integer(hid_t), intent(in) :: vtkhdf_grp
    integer, intent(in) :: counter
    real(kind=rp), intent(in) :: t
 
    integer(hid_t) :: xf_id, H5T_NEKO_DOUBLE
    integer :: ierr
    integer(hid_t) :: grp_id, dset_id, dcpl_id, filespace, memspace, attr_id
    integer(hsize_t), dimension(1) :: step_dims, step_maxdims
    integer(hsize_t), dimension(1) :: step_count, step_offset, chunkdims, ddim
    real(kind=dp), dimension(1) :: time_value
    integer(kind=i8) :: i8_value
    logical :: exists, attr_exists
 
    ! Create collective transfer property list
    call h5pcreate_f(h5p_dataset_xfer_f, xf_id, ierr)
    call h5pset_dxpl_mpio_f(xf_id, h5fd_mpio_collective_f, ierr)
    h5t_neko_double = h5kind_to_type(dp, h5_real_kind)
 
    ! Create or open Steps group
    call h5lexists_f(vtkhdf_grp, "Steps", exists, ierr)
    if (exists) then
       call h5gopen_f(vtkhdf_grp, "Steps", grp_id, ierr)
    else
       call h5gcreate_f(vtkhdf_grp, "Steps", grp_id, ierr)
    end if
 
    ! --- Values dataset (time values, real type) ---
    call h5lexists_f(grp_id, "Values", exists, ierr)
    if (exists) then
       call h5dopen_f(grp_id, "Values", dset_id, ierr)
       call h5dget_space_f(dset_id, filespace, ierr)
       call h5sget_simple_extent_dims_f(filespace, step_dims, step_maxdims, &
            ierr)
       call h5sclose_f(filespace, ierr)
 
       ! We have not written this timestep yet, expand the array
       if (step_dims(1) .eq. int(counter, hsize_t)) then
          step_dims(1) = int(counter + 1, hsize_t)
          call h5dset_extent_f(dset_id, step_dims, ierr)
       else if (step_dims(1) .lt. int(counter, hsize_t)) then
          call neko_error("VTKHDF: Time steps written out of order.")
       end if
    else
       step_dims(1) = 1_hsize_t
       step_maxdims(1) = h5s_unlimited_f
       chunkdims(1) = 1_hsize_t
 
       call h5screate_simple_f(1, step_dims, filespace, ierr, step_maxdims)
       call h5pcreate_f(h5p_dataset_create_f, dcpl_id, ierr)
       call h5pset_chunk_f(dcpl_id, 1, chunkdims, ierr)
       call h5dcreate_f(grp_id, "Values", h5t_neko_double, &
            filespace, dset_id, ierr, dcpl_id = dcpl_id)
       call h5sclose_f(filespace, ierr)
       call h5pclose_f(dcpl_id, ierr)
    end if
 
    step_count(1) = 1_hsize_t
    step_offset(1) = int(counter, hsize_t)
 
    call h5dget_space_f(dset_id, filespace, ierr)
    call h5screate_simple_f(1, step_count, memspace, ierr)
    call h5sselect_hyperslab_f(filespace, h5s_select_set_f, &
         step_offset, step_count, ierr)
 
    time_value(1) = real(t, kind=dp)
    call h5dwrite_f(dset_id, h5t_neko_double, time_value, step_count, ierr, &
         file_space_id = filespace, mem_space_id = memspace, xfer_prp = xf_id)
 
    call h5sclose_f(memspace, ierr)
    call h5sclose_f(filespace, ierr)
    call h5dclose_f(dset_id, ierr)
 
    ! --- NSteps attribute ---
    ddim(1) = 1_hsize_t
    call h5aexists_f(grp_id, "NSteps", attr_exists, ierr)
    if (attr_exists) then
       call h5aopen_f(grp_id, "NSteps", attr_id, ierr)
    else
       call h5screate_f(h5s_scalar_f, filespace, ierr)
       call h5acreate_f(grp_id, "NSteps", h5t_native_integer, filespace, &
            attr_id, ierr, h5p_default_f, h5p_default_f)
       call h5sclose_f(filespace, ierr)
    end if
 
    call h5awrite_f(attr_id, h5t_native_integer, counter + 1, ddim, ierr)
 
    call h5aclose_f(attr_id, ierr)
    call h5gclose_f(grp_id, ierr)
    call h5pclose_f(xf_id, ierr)
 
  end subroutine vtkhdf_write_steps
 
  subroutine vtkhdf_write_pointdata(vtkhdf_grp, fields, precision, counter, &
       fname, t)
    integer(hid_t), intent(in) :: vtkhdf_grp
    type(field_list_t), intent(inout) :: fields
    integer, intent(in) :: precision
    integer, intent(in) :: counter
    character(len=*), intent(in) :: fname
    real(kind=rp), intent(in), optional :: t
 
    integer(kind=i8) :: time_offset
    integer :: local_points, point_offset, total_points
    integer(hid_t) :: precision_hdf
    integer :: ierr, i, j
    integer :: n_fields
    integer(hid_t) :: pointdata_grp, grp_id, step_grp_id
    integer(hid_t) :: dset_id, dcpl_id, filespace
    integer(hsize_t), dimension(1) :: pd_dims1, pd_maxdims1
    integer(hsize_t), dimension(2) :: pd_dims2, pd_maxdims2
    type(field_t), pointer :: u, v, w
    character(len=128) :: field_name
    logical :: exists, is_vector
 
    ! VDS and per-timestep external file variables
    character(len=1024) :: ext_fname, ext_path, src_pattern
    character(len=1024) :: main_path, main_name, main_suffix
    integer(hid_t) :: ext_file_id, ext_plist_id, vds_src_space
    integer(hid_t) :: write_target, attr_id, H5T_NEKO_STRING
    integer :: mpi_info, mpi_comm
 
    ! Collected field info for VDS phase
    integer :: fields_written
    character(len=128), allocatable :: name_list(:)
    logical, allocatable :: vector_list(:)
 
    mpi_info = mpi_info_null%mpi_val
    mpi_comm = neko_comm%mpi_val
 
    n_fields = fields%size()
 
    ! Compute local/global point counts and MPI offsets
    local_points = fields%item_size(1)
    total_points = fields%items(1)%ptr%dof%global_size()
    point_offset = 0
    call mpi_exscan(local_points, point_offset, 1, mpi_integer, &
         mpi_sum, neko_comm, ierr)
 
    ! Sync all the fields
    do i = 1, n_fields
       if (associated(fields%items(i)%ptr)) then
          call fields%items(i)%ptr%copy_from(device_to_host, sync = i .eq. n_fields)
       end if
    end do
 
    fields_written = 0
    allocate(name_list(n_fields))
    allocate(vector_list(n_fields))
 
    ! Create PointData group if missing
    call h5lexists_f(vtkhdf_grp, "PointData", exists, ierr)
    if (.not. exists) then
       call h5gcreate_f(vtkhdf_grp, "PointData", pointdata_grp, ierr)
       call h5gclose_f(pointdata_grp, ierr)
    end if
 
    ! ------------------------------------------------------------------------ !
    ! Construct the target where data is written
 
    if (present(t)) then
 
       ! Derive base path from main filename for external files
       call filename_split(fname, main_path, main_name, main_suffix)
       write(ext_path, '(A,A,".data/")') trim(main_path), trim(main_name)
       write(ext_fname, '(A,I0,".h5")') trim(ext_path), counter
       write(src_pattern, '(A,".data/%b.h5")') trim(main_name)
 
       if (pe_rank .eq. 0) then
          call mkdir(trim(ext_path))
       end if
       call mpi_barrier(neko_comm, ierr)
 
       call h5pcreate_f(h5p_file_access_f, ext_plist_id, ierr)
       call h5pset_fapl_mpio_f(ext_plist_id, mpi_comm, mpi_info, ierr)
       call h5fcreate_f(trim(ext_fname), h5f_acc_trunc_f, write_target, ierr, &
            access_prp = ext_plist_id)
       call h5pclose_f(ext_plist_id, ierr)
 
    else
       ! Non-temporal: write directly into the main file's PointData group
       call h5gopen_f(vtkhdf_grp, "PointData", write_target, ierr)
    end if
 
    ! ------------------------------------------------------------------------ !
    ! Write field data
 
    do i = 1, n_fields
       field_name = fields%name(i)
       if (field_name .eq. 'p') field_name = 'Pressure'
 
       ! Determine if this is a velocity component to group as a vector
       is_vector = .false.
       if (field_name .eq. 'u' .or. field_name .eq. 'v' .or. &
            field_name .eq. 'w') then
          u => null()
          v => null()
          w => null()
          do j = 1, n_fields
             select case (trim(fields%name(j)))
             case ('u')
                u => fields%get(j)
             case ('v')
                v => fields%get(j)
             case ('w')
                w => fields%get(j)
             end select
          end do
 
          if (associated(u) .and. associated(v) .and. associated(w)) then
             is_vector = .true.
             field_name = 'Velocity'
          end if
       end if
 
       ! Skip duplicate fields (e.g. fluid_rho added by both fluid output
       ! and field_writer when sharing the same output file)
       exists = .false.
       do j = 1, fields_written
          if (trim(name_list(j)) .eq. trim(field_name)) then
             exists = .true.
             exit
          end if
       end do
 
       ! Skip duplicate fields
       if (exists) cycle
 
       ! Track unique field names for VDS phase
       fields_written = fields_written + 1
       name_list(fields_written) = field_name
       vector_list(fields_written) = is_vector
 
       ! Write field data to the target
       if (is_vector) then
          call write_vector_field(write_target, field_name, u%x, v%x, w%x, &
               local_points, precision, total_points, point_offset)
       else
          call write_scalar_field(write_target, field_name, fields%x(i), &
               local_points, precision, total_points, point_offset)
       end if
    end do
 
    ! Close write target
    if (present(t)) then
       call h5fclose_f(write_target, ierr)
    else
       call h5gclose_f(write_target, ierr)
    end if
 
    ! ------------------------------------------------------------------------ !
    ! Manage temporal datasets through VDS
 
    if (present(t)) then
 
       ! Temporal: set up per-timestep external file as write target
       call h5gopen_f(vtkhdf_grp, "Steps", step_grp_id, ierr)
       time_offset = int(counter, kind=i8) * int(total_points, kind=i8)
 
       ! Write PointDataOffsets under Steps for each unique field
       call h5lexists_f(step_grp_id, "PointDataOffsets", exists, ierr)
       if (exists) then
          call h5gopen_f(step_grp_id, "PointDataOffsets", grp_id, ierr)
       else
          call h5gcreate_f(step_grp_id, "PointDataOffsets", grp_id, ierr)
       end if
       do i = 1, fields_written
          call vtkhdf_write_i8_at(grp_id, trim(name_list(i)), &
               time_offset, counter)
       end do
       call h5gclose_f(grp_id, ierr)
       call h5gclose_f(step_grp_id, ierr)
 
       ! ===== Create or extend VDS in main file =====
       call h5gopen_f(vtkhdf_grp, "PointData", pointdata_grp, ierr)
 
       do i = 1, fields_written
          field_name = name_list(i)
          is_vector = vector_list(i)
 
          if (counter .eq. 0) then
             ! First write: create VDS with pattern-based mapping
             call h5pcreate_f(h5p_dataset_create_f, dcpl_id, ierr)
             precision_hdf = h5kind_to_type(precision, h5_real_kind)
 
             if (is_vector) then
                pd_dims2 = [3_hsize_t, int(total_points, hsize_t)]
                call h5screate_simple_f(2, pd_dims2, vds_src_space, ierr)
                call h5sselect_all_f(vds_src_space, ierr)
 
                pd_maxdims2 = [3_hsize_t, h5s_unlimited_f]
                call h5screate_simple_f(2, pd_dims2, filespace, ierr, &
                     pd_maxdims2)
 
                call h5sselect_hyperslab_f(filespace, h5s_select_set_f, &
                     [0_hsize_t, 0_hsize_t], &
                     [1_hsize_t, h5s_unlimited_f], &
                     ierr, &
                     stride = [3_hsize_t, int(total_points, hsize_t)], &
                     block = [3_hsize_t, int(total_points, hsize_t)])
 
                call h5pset_virtual_f(dcpl_id, filespace, trim(src_pattern), &
                     trim(field_name), vds_src_space, ierr)
                call h5sclose_f(vds_src_space, ierr)
 
                call h5dcreate_f(pointdata_grp, trim(field_name), &
                     precision_hdf, filespace, dset_id, ierr, &
                     dcpl_id = dcpl_id)
                call h5sclose_f(filespace, ierr)
             else
                pd_dims1 = int(total_points, hsize_t)
                call h5screate_simple_f(1, pd_dims1, vds_src_space, ierr)
                call h5sselect_all_f(vds_src_space, ierr)
 
                pd_maxdims1(1) = h5s_unlimited_f
                call h5screate_simple_f(1, pd_dims1, filespace, ierr, &
                     pd_maxdims1)
 
                call h5sselect_hyperslab_f(filespace, h5s_select_set_f, &
                     [0_hsize_t], &
                     [h5s_unlimited_f], &
                     ierr, &
                     stride = [int(total_points, hsize_t)], &
                     block = [int(total_points, hsize_t)])
 
                call h5pset_virtual_f(dcpl_id, filespace, trim(src_pattern), &
                     trim(field_name), vds_src_space, ierr)
                call h5sclose_f(vds_src_space, ierr)
 
                call h5dcreate_f(pointdata_grp, trim(field_name), &
                     precision_hdf, filespace, dset_id, ierr, &
                     dcpl_id = dcpl_id)
                call h5sclose_f(filespace, ierr)
             end if
 
             call h5pclose_f(dcpl_id, ierr)
             call h5dclose_f(dset_id, ierr)
 
          else
             ! Subsequent write: extend VDS to include new timestep
             call h5dopen_f(pointdata_grp, trim(field_name), dset_id, ierr)
 
             if (is_vector) then
                pd_dims2 = [3_hsize_t, &
                     int(counter + 1, hsize_t) * int(total_points, hsize_t)]
                call h5dset_extent_f(dset_id, pd_dims2, ierr)
             else
                pd_dims1 = int(counter + 1, hsize_t) * int(total_points, hsize_t)
                call h5dset_extent_f(dset_id, pd_dims1, ierr)
             end if
 
             call h5dclose_f(dset_id, ierr)
          end if
       end do
 
       call h5gclose_f(pointdata_grp, ierr)
    end if
 
    ! ------------------------------------------------------------------------ !
    ! Add attributes to the PointData
 
    do i = 1, fields_written
       field_name = name_list(i)
       is_vector = vector_list(i)
       pd_dims1 = 1_hsize_t
 
       call h5gopen_f(vtkhdf_grp, "PointData", pointdata_grp, ierr)
       call h5dopen_f(pointdata_grp, trim(field_name), dset_id, ierr)
 
       call h5aexists_f(dset_id, "Attribute", exists, ierr)
       if (exists) then
          call h5dclose_f(dset_id, ierr)
          call h5gclose_f(pointdata_grp, ierr)
          cycle
       end if
 
       ! Write the attribute what this PointData represents
       call h5screate_f(h5s_scalar_f, filespace, ierr)
 
       call h5tcopy_f(h5t_fortran_s1, h5t_neko_string, ierr)
       call h5tset_size_f(h5t_neko_string, int(6, size_t), ierr)
       call h5tset_strpad_f(h5t_neko_string, h5t_str_nullterm_f, ierr)
 
       call h5acreate_f(dset_id, "Attribute", h5t_neko_string, filespace, &
            attr_id, ierr)
       if (is_vector) then
          call h5awrite_f(attr_id, h5t_neko_string, ["Vector"], pd_dims1, ierr)
       else
          call h5awrite_f(attr_id, h5t_neko_string, ["Scalar"], pd_dims1, ierr)
       end if
 
       call h5aclose_f(attr_id, ierr)
       call h5tclose_f(h5t_neko_string, ierr)
       call h5sclose_f(filespace, ierr)
       call h5dclose_f(dset_id, ierr)
       call h5gclose_f(pointdata_grp, ierr)
 
    end do
 
    ! ------------------------------------------------------------------------ !
    ! Cleanup before returning
 
    deallocate(name_list)
    deallocate(vector_list)
 
  end subroutine vtkhdf_write_pointdata
 
  ! -------------------------------------------------------------------------- !
  ! Helper functions and routines
 
  subroutine vtkhdf_build_connectivity(conn, vtk_type, msh, dof, subdivide)
    integer, intent(inout) :: conn(:)
    integer(kind=1), intent(in) :: vtk_type
    type(mesh_t), intent(in) :: msh
    type(dofmap_t), intent(in) :: dof
    logical, intent(in) :: subdivide
    integer :: lx, ly, lz, nelv
    integer :: ie, ii, n_pts_per_elem, n_conn_per_elem
    integer, allocatable :: node_order(:)
 
    nelv = msh%nelv
    lx = dof%Xh%lx
    ly = dof%Xh%ly
    lz = dof%Xh%lz
    n_pts_per_elem = lx * ly * lz
 
    if (subdivide .and. vtk_type .eq. int(12, kind=1)) then
       node_order = subdivide_to_hex_ordering(lx, ly, lz)
    else if (subdivide .and. vtk_type .eq. int(9, kind=1)) then
       node_order = subdivide_to_quad_ordering(lx, ly)
    else
       node_order = vtk_ordering(vtk_type, lx, ly, lz)
    end if
 
    n_conn_per_elem = size(node_order)
 
    do concurrent(ie = 1:nelv, ii = 1:n_conn_per_elem)
       block
         integer :: idx, base
         idx = (ie - 1) * n_conn_per_elem
         base = (ie - 1) * n_pts_per_elem
         conn(idx + ii) = base + node_order(ii)
       end block
    end do
 
    deallocate(node_order)
 
  end subroutine vtkhdf_build_connectivity
 
  subroutine vtkhdf_write_numberof(grp, dset_name, value, offset, cnt, &
       dcpl, index, xf_id, ierr)
    integer(hid_t), intent(in) :: grp, dcpl, xf_id
    character(len=*), intent(in) :: dset_name
    integer(kind=i8), intent(in) :: value
    integer, intent(in):: index
    integer(hsize_t), dimension(1), intent(in) :: offset, cnt
    integer, intent(out) :: ierr
 
    integer(hid_t) :: dset_id, fspace, mspace
    integer(hsize_t), dimension(1) :: dims, maxdims
    integer(hid_t) :: H5T_NEKO_INTEGER
    logical :: exists
 
    h5t_neko_integer = h5kind_to_type(i8, h5_integer_kind)
 
    call h5lexists_f(grp, dset_name, exists, ierr)
    if (exists) then
       call h5dopen_f(grp, dset_name, dset_id, ierr)
       dims(1) = int(index + 1, hsize_t) * int(pe_size, hsize_t)
       call h5dset_extent_f(dset_id, dims, ierr)
    else
       dims(1) = int(pe_size, hsize_t)
       maxdims(1) = h5s_unlimited_f
       call h5screate_simple_f(1, dims, fspace, ierr, maxdims)
       call h5dcreate_f(grp, dset_name, h5t_neko_integer, &
            fspace, dset_id, ierr, dcpl_id = dcpl)
       call h5sclose_f(fspace, ierr)
    end if
 
    call h5dget_space_f(dset_id, fspace, ierr)
    call h5screate_simple_f(1, cnt, mspace, ierr)
    call h5sselect_hyperslab_f(fspace, h5s_select_set_f, offset, cnt, ierr)
 
    call h5dwrite_f(dset_id, h5t_neko_integer, value, cnt, ierr, &
         file_space_id = fspace, mem_space_id = mspace, xfer_prp = xf_id)
 
    call h5sclose_f(mspace, ierr)
    call h5sclose_f(fspace, ierr)
    call h5dclose_f(dset_id, ierr)
  end subroutine vtkhdf_write_numberof
 
  subroutine vtkhdf_write_i8_at(grp_id, name, value, index)
    integer(hid_t), intent(in) :: grp_id
    character(len=*), intent(in) :: name
    integer(kind=i8), intent(in) :: value
    integer, intent(in) :: index
 
    integer :: ierr
    integer(hid_t) :: dset_id, dcpl_id, xf_id, filespace, memspace
    integer(hsize_t), dimension(1) :: dims, maxdims, count, offset, chunkdims
    integer(hid_t) :: H5T_NEKO_INTEGER
    logical :: exists
 
    h5t_neko_integer = h5kind_to_type(i8, h5_integer_kind)
 
    ! Create collective transfer property list
    call h5pcreate_f(h5p_dataset_xfer_f, xf_id, ierr)
    call h5pset_dxpl_mpio_f(xf_id, h5fd_mpio_collective_f, ierr)
 
    call h5lexists_f(grp_id, name, exists, ierr)
    if (exists) then
       call h5dopen_f(grp_id, name, dset_id, ierr)
       call h5dget_space_f(dset_id, filespace, ierr)
       call h5sget_simple_extent_dims_f(filespace, dims, maxdims, ierr)
       call h5sclose_f(filespace, ierr)
 
       if (int(index, hsize_t) .eq. dims(1)) then
          dims(1) = int(index + 1, hsize_t)
          call h5dset_extent_f(dset_id, dims, ierr)
       else if (int(index, hsize_t) .gt. dims(1)) then
          call neko_error("VTKHDF: Values written out of order.")
       end if
    else
       dims = 1_hsize_t
       maxdims = h5s_unlimited_f
       chunkdims = 1_hsize_t
 
       call h5screate_simple_f(1, dims, filespace, ierr, maxdims)
       call h5pcreate_f(h5p_dataset_create_f, dcpl_id, ierr)
       call h5pset_chunk_f(dcpl_id, 1, chunkdims, ierr)
       call h5dcreate_f(grp_id, name, h5t_neko_integer, &
            filespace, dset_id, ierr, dcpl_id = dcpl_id)
       call h5sclose_f(filespace, ierr)
       call h5pclose_f(dcpl_id, ierr)
    end if
 
    count = 1_hsize_t
    offset = int(index, hsize_t)
 
    call h5dget_space_f(dset_id, filespace, ierr)
    call h5screate_simple_f(1, count, memspace, ierr)
    call h5sselect_hyperslab_f(filespace, h5s_select_set_f, offset, count, ierr)
    call h5dwrite_f(dset_id, h5t_neko_integer, value, count, ierr, &
         file_space_id = filespace, mem_space_id = memspace, xfer_prp = xf_id)
 
    call h5sclose_f(memspace, ierr)
    call h5sclose_f(filespace, ierr)
    call h5dclose_f(dset_id, ierr)
    call h5pclose_f(xf_id, ierr)
 
  end subroutine vtkhdf_write_i8_at
 
  subroutine write_scalar_field(hdf_root, name, x, n_local, &
       precision, n_total, offset)
    integer, intent(in) :: n_local
    integer(hid_t), intent(in) :: hdf_root
    character(len=*), intent(in) :: name
    real(kind=rp), dimension(n_local), intent(in) :: x
    integer, intent(in), optional :: precision
    integer, intent(in), optional :: n_total, offset
 
    integer(hsize_t), dimension(1) :: dims, dcount, doffset
    integer(hid_t) :: xf_id, dset_id, filespace, memspace, precision_hdf
    integer :: i, ierr, precision_local, n_tot, off
 
    ! Setup data sizes, offsets and precision
    dcount = int(n_local, hsize_t)
 
    if (present(n_total)) then
       dims = int(n_total, hsize_t)
    else
       call mpi_allreduce(n_local, n_tot, 1, mpi_integer, mpi_sum, neko_comm, &
            ierr)
       dims = int(n_tot, hsize_t)
    end if
 
    if (present(offset)) then
       doffset = int(offset, hsize_t)
    else
       call mpi_exscan(n_local, off, 1, mpi_integer, mpi_sum, neko_comm, &
            ierr)
       doffset = int(off, hsize_t)
    end if
 
    if (present(precision)) then
       precision_local = precision
    else
       precision_local = rp
    end if
    precision_hdf = h5kind_to_type(precision_local, h5_real_kind)
 
    ! Prepare memory and filespaces
    call h5pcreate_f(h5p_dataset_xfer_f, xf_id, ierr)
    call h5pset_dxpl_mpio_f(xf_id, h5fd_mpio_collective_f, ierr)
 
    call h5screate_simple_f(1, dims, filespace, ierr)
    call h5screate_simple_f(1, dcount, memspace, ierr)
    call h5sselect_hyperslab_f(filespace, h5s_select_set_f, &
         doffset, dcount, ierr)
 
    call h5dcreate_f(hdf_root, trim(name), precision_hdf, &
         filespace, dset_id, ierr)
 
    if (precision_local .eq. rp) then
       call h5dwrite_f(dset_id, precision_hdf, x, dcount, ierr, &
            file_space_id = filespace, mem_space_id = memspace, &
            xfer_prp = xf_id)
 
    else if (precision_local .eq. sp) then
       block
         real(kind=sp), allocatable :: x_sp(:)
 
         allocate(x_sp(n_local))
         do concurrent(i = 1:n_local)
            x_sp(i) = real(x(i), sp)
         end do
 
         call h5dwrite_f(dset_id, precision_hdf, x_sp, dcount, ierr, &
              file_space_id = filespace, mem_space_id = memspace, &
              xfer_prp = xf_id)
         deallocate(x_sp)
       end block
 
    else if (precision_local .eq. dp) then
       block
         real(kind=dp), allocatable :: x_dp(:)
 
         allocate(x_dp(n_local))
         do concurrent(i = 1:n_local)
            x_dp(i) = real(x(i), dp)
         end do
 
         call h5dwrite_f(dset_id, precision_hdf, x_dp, dcount, ierr, &
              file_space_id = filespace, mem_space_id = memspace, &
              xfer_prp = xf_id)
         deallocate(x_dp)
       end block
 
    else if (precision_local .eq. qp) then
       block
         real(kind=qp), allocatable :: x_qp(:)
 
         allocate(x_qp(n_local))
         do concurrent(i = 1:n_local)
            x_qp(i) = real(x(i), qp)
         end do
 
         call h5dwrite_f(dset_id, precision_hdf, x_qp, dcount, ierr, &
              file_space_id = filespace, mem_space_id = memspace, &
              xfer_prp = xf_id)
         deallocate(x_qp)
       end block
 
    else
       call neko_error("Unsupported precision in HDF5 write_scalar_field")
    end if
 
    call h5sclose_f(filespace, ierr)
    call h5sclose_f(memspace, ierr)
    call h5dclose_f(dset_id, ierr)
    call h5pclose_f(xf_id, ierr)
  end subroutine write_scalar_field
 
  subroutine write_vector_field(hdf_root, name, u, v, w, n_local, &
       precision, n_total, offset)
    integer, intent(in) :: n_local
    integer(hid_t), intent(in) :: hdf_root
    character(len=*), intent(in) :: name
    real(kind=rp), dimension(n_local), intent(in) :: u, v, w
    integer, intent(in), optional :: precision
    integer, intent(in), optional :: n_total, offset
 
    integer(hsize_t), dimension(2) :: dims, dcount, doffset
    integer(hid_t) :: xf_id, dset_id, filespace, memspace, precision_hdf
    integer :: i, ierr, precision_local, n_tot, off
 
    ! Setup data sizes, offsets and precision
    dcount = [3_hsize_t, int(n_local, hsize_t)]
 
    if (present(n_total)) then
       dims = [3_hsize_t, int(n_total, hsize_t)]
    else
       call mpi_allreduce(n_local, n_tot, 1, mpi_integer, mpi_sum, neko_comm, &
            ierr)
       dims = [3_hsize_t, int(n_tot, hsize_t)]
    end if
 
    if (present(offset)) then
       doffset = [0_hsize_t, int(offset, hsize_t)]
    else
       call mpi_exscan(n_local, off, 1, mpi_integer, mpi_sum, neko_comm, &
            ierr)
       doffset = [0_hsize_t, int(off, hsize_t)]
    end if
 
    if (present(precision)) then
       precision_local = precision
    else
       precision_local = rp
    end if
    precision_hdf = h5kind_to_type(precision_local, h5_real_kind)
 
    ! Prepare memory and filespaces
    call h5pcreate_f(h5p_dataset_xfer_f, xf_id, ierr)
    call h5pset_dxpl_mpio_f(xf_id, h5fd_mpio_collective_f, ierr)
 
    call h5screate_simple_f(2, dims, filespace, ierr)
    call h5screate_simple_f(2, dcount, memspace, ierr)
    call h5sselect_hyperslab_f(filespace, h5s_select_set_f, &
         doffset, dcount, ierr)
 
    call h5dcreate_f(hdf_root, trim(name), precision_hdf, filespace, dset_id, &
         ierr)
 
    if (precision_local .eq. sp) then
       block
         real(kind=sp), allocatable :: f(:,:)
 
         allocate(f(3, n_local))
         do concurrent(i = 1:n_local)
            f(1, i) = real(u(i), sp)
            f(2, i) = real(v(i), sp)
            f(3, i) = real(w(i), sp)
         end do
 
         call h5dwrite_f(dset_id, precision_hdf, f, dcount, ierr, &
              file_space_id = filespace, mem_space_id = memspace, &
              xfer_prp = xf_id)
         deallocate(f)
       end block
 
    else if (precision_local .eq. dp) then
       block
         real(kind=dp), allocatable :: f(:,:)
 
         allocate(f(3, n_local))
         do concurrent(i = 1:n_local)
            f(1, i) = real(u(i), dp)
            f(2, i) = real(v(i), dp)
            f(3, i) = real(w(i), dp)
         end do
 
         call h5dwrite_f(dset_id, precision_hdf, f, dcount, ierr, &
              file_space_id = filespace, mem_space_id = memspace, &
              xfer_prp = xf_id)
         deallocate(f)
       end block
 
    else if (precision_local .eq. qp) then
       block
         real(kind=qp), allocatable :: f(:,:)
 
         allocate(f(3, n_local))
         do concurrent(i = 1:n_local)
            f(1, i) = real(u(i), qp)
            f(2, i) = real(v(i), qp)
            f(3, i) = real(w(i), qp)
         end do
 
         call h5dwrite_f(dset_id, precision_hdf, f, dcount, ierr, &
              file_space_id = filespace, mem_space_id = memspace, &
              xfer_prp = xf_id)
         deallocate(f)
       end block
 
    else
       call neko_error("Unsupported precision in HDF5 write_vector_field")
    end if
 
    call h5sclose_f(filespace, ierr)
    call h5sclose_f(memspace, ierr)
    call h5dclose_f(dset_id, ierr)
    call h5pclose_f(xf_id, ierr)
  end subroutine write_vector_field
 
  ! -------------------------------------------------------------------------- !
  ! Reader functions and routines
 
  subroutine vtkhdf_file_read(this, data)
    class(vtkhdf_file_t) :: this
    class(*), target, intent(inout) :: data
    character(len=1024) :: fname
    integer(hid_t) :: plist_id, file_id, vtkhdf_grp
    integer :: mpi_info, mpi_comm, counter, ierr
    type(field_list_t) :: fields
    integer :: i
 
    ! Open the file
    fname = trim(this%get_base_fname())
    counter = this%get_counter() - this%get_start_counter()
    if (counter .lt. 0) counter = 0
 
    mpi_info = mpi_info_null%mpi_val
    mpi_comm = neko_comm%mpi_val
 
    call h5open_f(ierr)
    call h5pcreate_f(h5p_file_access_f, plist_id, ierr)
    call h5pset_fapl_mpio_f(plist_id, mpi_comm, mpi_info, ierr)
 
    call h5fopen_f(fname, h5f_acc_rdonly_f, file_id, ierr, &
         access_prp = plist_id)
    if (ierr .ne. 0) then
       call neko_error('Error opening VTKHDF file: ' // trim(fname))
    end if
 
    call h5gopen_f(file_id, "VTKHDF", vtkhdf_grp, ierr)
    if (ierr .ne. 0) then
       call h5fclose_f(file_id, ierr)
       call h5pclose_f(plist_id, ierr)
       call neko_error('VTKHDF group not found in file: ' // trim(fname))
    end if
 
    select type (data)
    type is (field_t)
       call fields%init(1)
       ! Specific bindings instead of the assign generic (see
       ! vtkhdf_file_write); Fujitsu frt cannot resolve the generic
       call fields%assign_to_field(1, data)
    type is (field_list_t)
       call fields%assign_to_list(data)
    class default
       call neko_error("Unsupported data type in vtkhdf_file_read")
    end select
 
    do i = 1, fields%size()
       call vtkhdf_read_field(vtkhdf_grp, fname, counter, fields%get(i))
    end do
    call fields%copy_from(host_to_device, .true.)
 
    call h5gclose_f(vtkhdf_grp, ierr)
    call h5fclose_f(file_id, ierr)
    call h5pclose_f(plist_id, ierr)
    call h5close_f(ierr)
 
    call fields%free()
 
  end subroutine vtkhdf_file_read
 
  subroutine vtkhdf_read_field(vtkhdf_grp, fname, counter, fld)
    integer(hid_t), intent(in) :: vtkhdf_grp
    character(len=*), intent(in) :: fname
    integer, intent(in) :: counter
    type(field_t), intent(inout) :: fld
 
    character(len=:), allocatable :: field_name
    integer :: ierr, local_points, total_points, point_offset, component
    integer(hid_t) :: pointdata_grp, dset_id, filespace, memspace, xf_id
    integer(hid_t) :: dcpl_id, ext_file_id, ext_plist_id, attr_id
    integer(hid_t) :: H5T_NEKO_REAL
    integer(hsize_t), dimension(1) :: dcount1, doffset1
    integer(hsize_t), dimension(2) :: dcount2, doffset2
 
    integer(size_t) :: vds_count
    character(len=128) :: dset_name
    character(len=1024) :: vds_src_file, ext_fname, main_path, main_name, main_suffix
    logical :: exists, is_vds
    real(kind=rp), allocatable :: vec_component(:,:)
    integer :: mpi_info, mpi_comm, pct_pos, nsteps
    character(len=256) :: error_message
 
    field_name = trim(fld%name)
 
    ! Validate counter against file's NSteps if temporal data exists
    call h5lexists_f(vtkhdf_grp, "Steps", exists, ierr)
    if (exists) then
       call h5aopen_by_name_f(vtkhdf_grp, "Steps", "NSteps", attr_id, ierr)
       call h5aread_f(attr_id, h5t_native_integer, nsteps, [1_hsize_t], ierr)
       call h5aclose_f(attr_id, ierr)
 
       if (counter .ge. nsteps) then
          write(error_message, '(A,I0,A,I0)') &
               'VTKHDF read: counter ', counter, ' >= NSteps ', nsteps
          call neko_error(trim(error_message))
       end if
    end if
 
    mpi_info = mpi_info_null%mpi_val
    mpi_comm = neko_comm%mpi_val
 
    local_points = fld%dof%size()
    total_points = fld%dof%global_size()
    point_offset = 0
    call mpi_exscan(local_points, point_offset, 1, mpi_integer, &
         mpi_sum, neko_comm, ierr)
 
    h5t_neko_real = h5kind_to_type(rp, h5_real_kind)
 
    ! Open PointData group
    call h5gopen_f(vtkhdf_grp, "PointData", pointdata_grp, ierr)
 
    ! Determine dataset name with standard mappings
    dset_name = trim(field_name)
    if (trim(dset_name) .eq. 'p') dset_name = 'Pressure'
 
    component = -1
    call h5lexists_f(pointdata_grp, trim(dset_name), exists, ierr)
 
    if (.not. exists) then
       if (trim(field_name) .eq. 'u' .or. trim(field_name) .eq. 'v' .or. &
            trim(field_name) .eq. 'w') then
          call h5lexists_f(pointdata_grp, "Velocity", exists, ierr)
          if (exists) then
             dset_name = 'Velocity'
             select case (trim(field_name))
             case ('u')
                component = 0
             case ('v')
                component = 1
             case ('w')
                component = 2
             end select
          end if
       end if
    end if
 
    if (.not. exists) then
       call h5gclose_f(pointdata_grp, ierr)
       call neko_error('VTKHDF PointData field not found: ' // trim(dset_name))
    end if
 
    ! Open the dataset and check if it's a VDS
    call h5dopen_f(pointdata_grp, trim(dset_name), dset_id, ierr)
    call h5dget_create_plist_f(dset_id, dcpl_id, ierr)
    call h5pget_virtual_count_f(dcpl_id, vds_count, ierr)
    is_vds = (vds_count .gt. 0_size_t)
 
    if (is_vds) then
       call filename_split(fname, main_path, main_name, main_suffix)
 
       ! Query mapping 0 to detect whether this is a printf-style pattern VDS
       ! (Neko's own format uses a single mapping with %b as a block counter
       ! placeholder) or a multi-mapping VDS where each timestep has its own
       ! literal source file.
       call h5pget_virtual_filename_f(dcpl_id, 0_size_t, vds_src_file, ierr)
       pct_pos = index(vds_src_file, '%b')
 
       if (pct_pos .gt. 0) then
          ! Pattern-based VDS: replace %b with the timestep counter.
          ! If the stored path is absolute, use it as-is; otherwise prepend the
          ! main file's directory so the path is resolved correctly regardless
          ! of the process working directory.
          if (vds_src_file(1:1) .eq. '/') then
             write(ext_fname, '(A,I0,A)') &
                  vds_src_file(1:pct_pos-1), counter, &
                  trim(vds_src_file(pct_pos+2:))
          else
             write(ext_fname, '(A,A,I0,A)') &
                  trim(main_path), vds_src_file(1:pct_pos-1), counter, &
                  trim(vds_src_file(pct_pos+2:))
          end if
       else
          ! Multi-mapping VDS: one mapping per timestep. Query the mapping
          ! that corresponds directly to the requested counter instead of
          ! always using mapping 0, which would silently read the wrong file.
          if (int(counter, size_t) .ge. vds_count) then
             call h5pclose_f(dcpl_id, ierr)
             call h5dclose_f(dset_id, ierr)
             call h5gclose_f(pointdata_grp, ierr)
             write(error_message, '(A,I0,A,I0)') &
                  'VTKHDF: VDS counter ', counter, &
                  ' is out of range, number of mappings: ', int(vds_count)
             call neko_error(trim(error_message))
          end if
          call h5pget_virtual_filename_f(dcpl_id, int(counter, size_t), &
               vds_src_file, ierr)
          if (vds_src_file(1:1) .eq. '/') then
             ext_fname = trim(vds_src_file)
          else
             write(ext_fname, '(A,A)') trim(main_path), trim(vds_src_file)
          end if
       end if
 
       call h5pclose_f(dcpl_id, ierr)
       call h5dclose_f(dset_id, ierr)
       call h5gclose_f(pointdata_grp, ierr)
 
       call h5pcreate_f(h5p_file_access_f, ext_plist_id, ierr)
       call h5pset_fapl_mpio_f(ext_plist_id, mpi_comm, mpi_info, ierr)
       call h5fopen_f(trim(ext_fname), h5f_acc_rdonly_f, ext_file_id, ierr, &
            access_prp = ext_plist_id)
 
       if (ierr .ne. 0) then
          call neko_error('VTKHDF: Cannot open VDS source file: ' // &
               trim(ext_fname))
       end if
 
       ! Open the dataset from the external file
       call h5dopen_f(ext_file_id, trim(dset_name), dset_id, ierr)
 
       call h5pclose_f(ext_plist_id, ierr)
    else
       call h5pclose_f(dcpl_id, ierr)
       ext_file_id = -1
    end if
 
    ! Now read from the dataset (either original or from external file)
    call h5pcreate_f(h5p_dataset_xfer_f, xf_id, ierr)
    call h5pset_dxpl_mpio_f(xf_id, h5fd_mpio_collective_f, ierr)
    call h5dget_space_f(dset_id, filespace, ierr)
 
    if (trim(dset_name) .eq. 'Velocity' .and. component .ge. 0) then
       ! Velocity is stored as (3, n_points), read one component
       dcount2 = [1_hsize_t, int(local_points, hsize_t)]
       doffset2 = [int(component, hsize_t), int(point_offset, hsize_t)]
 
       call h5screate_simple_f(2, dcount2, memspace, ierr)
       call h5sselect_hyperslab_f(filespace, h5s_select_set_f, &
            doffset2, dcount2, ierr)
 
       allocate(vec_component(1, local_points))
       call h5dread_f(dset_id, h5t_neko_real, vec_component, dcount2, ierr, &
            file_space_id = filespace, mem_space_id = memspace, &
            xfer_prp = xf_id)
       fld%x = reshape(vec_component, shape(fld%x))
       deallocate(vec_component)
    else
       ! Scalar field: 1D array
       dcount1 = int(local_points, hsize_t)
       doffset1 = int(point_offset, hsize_t)
 
       call h5screate_simple_f(1, dcount1, memspace, ierr)
       call h5sselect_hyperslab_f(filespace, h5s_select_set_f, &
            doffset1, dcount1, ierr)
 
       call h5dread_f(dset_id, h5t_neko_real, fld%x(1,1,1,1), dcount1, ierr, &
            file_space_id = filespace, mem_space_id = memspace, &
            xfer_prp = xf_id)
    end if
 
    call h5sclose_f(memspace, ierr)
    call h5sclose_f(filespace, ierr)
    call h5dclose_f(dset_id, ierr)
    call h5pclose_f(xf_id, ierr)
 
    if (is_vds) then
       call h5fclose_f(ext_file_id, ierr)
    else
       call h5gclose_f(pointdata_grp, ierr)
    end if
 
  end subroutine vtkhdf_read_field
 
#else
  ! -------------------------------------------------------------------------- !
  ! Dummy functions and subroutines
 
  subroutine vtkhdf_file_write(this, data, t)
    class(vtkhdf_file_t), intent(inout) :: this
    class(*), target, intent(in) :: data
    real(kind=rp), intent(in), optional :: t
    call neko_error('Neko needs to be built with HDF5 support')
  end subroutine vtkhdf_file_write
 
  subroutine vtkhdf_file_read(this, data)
    class(vtkhdf_file_t) :: this
    class(*), target, intent(inout) :: data
    call neko_error('Neko needs to be built with HDF5 support')
  end subroutine vtkhdf_file_read
 
#endif
 
  ! -------------------------------------------------------------------------- !
  ! Sub-cell node ordering functions for VTK compatibility
 
  pure function subdivide_to_hex_ordering(lx, ly, lz) result(node_order)
    integer, intent(in) :: lx, ly, lz
    integer :: node_order(8 * (lx - 1) * (ly - 1) * (lz - 1))
    integer :: ii, jj, kk, idx
 
    idx = 0
 
    do ii = 1, lx - 1
       do jj = 1, ly - 1
          do kk = 1, lz - 1
             node_order(idx + 1) = (kk - 1) * lx * ly + (jj - 1) * lx + ii - 1
             node_order(idx + 2) = (kk - 1) * lx * ly + (jj - 1) * lx + ii
             node_order(idx + 3) = (kk - 1) * lx * ly + jj * lx + ii
             node_order(idx + 4) = (kk - 1) * lx * ly + jj * lx + ii - 1
             node_order(idx + 5) = kk * lx * ly + (jj - 1) * lx + ii - 1
             node_order(idx + 6) = kk * lx * ly + (jj - 1) * lx + ii
             node_order(idx + 7) = kk * lx * ly + jj * lx + ii
             node_order(idx + 8) = kk * lx * ly + jj * lx + ii - 1
             idx = idx + 8
          end do
       end do
    end do
 
  end function subdivide_to_hex_ordering
 
  pure function subdivide_to_quad_ordering(lx, ly) result(node_order)
    integer, intent(in) :: lx, ly
    integer :: node_order(4 * (lx - 1) * (ly - 1))
    integer :: ii, jj, idx
 
    idx = 0
 
    do jj = 1, ly - 1
       do ii = 1, lx - 1
          node_order(idx + 1) = (jj - 1) * lx + ii - 1
          node_order(idx + 2) = (jj - 1) * lx + ii
          node_order(idx + 3) = jj * lx + ii
          node_order(idx + 4) = jj * lx + ii - 1
          idx = idx + 4
       end do
    end do
 
  end function subdivide_to_quad_ordering
 
end module vtkhdf_file