I/O from XDMF/HDF5 files in dolfin-x

I’ve not seen this error message before.
I’ve merged in the complex support code in the main branch (https://github.com/jorgensd/adios4dolfinx/pull/31).
However, to make sure that it works properly you should update the build of dolfinx to the nightly one (as of last night), as it contains a bug-fix I encountered while extending the code: Add missing templating of inverse transform by jorgensd · Pull Request #2778 · FEniCS/dolfinx · GitHub

The tests now contain tons of example of how to write different functions to file, and I believe it is thoroughly tested (500 test cases, variations of function spaces, polynomial degrees, real/complex and precision of mesh geometry/function).

1 Like

I have no idea what it was but I deleted the old output files manually and reran it and it worked.

I got it and it works beautifully! Thanks a lot again for the fast upgrade!

@marchirschvogel I’ve now made a prototype for reading and writing meshtags as standalone functions:

will put them into adios4dolfinx tomorrow.

1 Like

Thank you, I will check it out!

I still have some problems understanding the need to re-read the mesh though.
So suppose I have a field to read that was written by one process in some preprocessing step (e.g. fiber directions), and subsequently another field that was written from a parallel layout (e.g. a pre-deformation state from another simulation).
In both cases, the mesh indices/node orders would be different, right? So the first field would get invalidated and the functions would need to be re-defined, right? Or do I miss something here?

Thanks, best,
Marc

You need to store all data you want to read with the mesh you want to read.

To read in function data, adios4dolfinx (and legacy dolfin) stores the dofmap.

The dofmap is dependent on the cell ordering.

If you write data with different sources, these cells will be ordered in different ways.

I see the point that data needs to be stored together with some information about the mesh/ordering, yes. But wouldn’t it be beneficial to also have a reader that can map its re-ordered data to the base mesh? Otherwise it won’t be possible to read functions from different sources, even though they have self-contained information (relating a node to a value).
So I currently have this workaround by querying the igi = msh.geometry.input_global_indices and then doing the mapping, which at least works if I don’t change the order of the mesh internally. Since in the end, I will always have to deal with data to read that comes from different sources (possibly even another program that provides nodal values, etc.)

Best,
Marc

How would one do that? It would mean a global search to figure out what nodes that are in each cell and match them (as global ordering is different depending on the number of processes used).

Why not just keep on expanding the checkpoint with more information, ie if you need more data, read it in through the mesh and data with the checkpoint create new data, and save all to a new checkpoint.

If you want to read data depending on nodal values for other programs, you would have to implement a method to match nodes by comparing coordinates.

That is out of scope for the checkpointing. I would probably use nonmatching mesh interpolation to map data between the meshes.

Edit
The capability to write multiple functions and time steps to a single file in adios4dolfinx is being introduced with: Time dependency and multiple functions per file by jorgensd · Pull Request #49 · jorgensd/adios4dolfinx · GitHub

@marchirschvogel I’ve released a new version of adios4dolfinx, v0.7.2 that features multi-function and multi-timestep capabilities, along with MeshTags-support and major performance improvements

@dokken Thanks, I’ll have a look! I guess this certainly will help to store input data in a compact way.

With regards to the other aspect of different functions stored by different processes/programs, I guess the desire would be to at least be able to read any field that is stored in the same ordering logic as the base mesh. So, any nodal data I write has the same ordering as the mesh I provide. Hence, the same logic of re-distributing/ordering the mesh upon read-in can be applied to the field provided (as .xdmf/.h5 file). I don’t know how feasible that would be, it just sounds somehow natural to me - especially when it comes to exchanging input data etc.

Best,
Marc

If you want to

  1. Read a mesh_0 from file
  2. Store a solution to file (f). let us call this mesh mesh_1.
  3. Read data from another file (based on the mesh_0) and use with f from mesh_1

You would have to store more information than I currently do. You would have to store the original cell index and the original geometry indices.

It is not infeasible, but it would slow down storage of checkpoints (and reading them) a bit, as one of the operations would have to deal with this reordering (reading a global array on every process).

The workflow I would suggest is shown in:

@marchirschvogel @Ekrem_Ekici
Functionality for writing checkpoints relating to XDMF input files has now been added to adios4dolfinx: Checkpointing to original mesh by jorgensd · Pull Request #70 · jorgensd/adios4dolfinx · GitHub

2 Likes

Is it possible to read in legacy MeshFunction mesh tags from a legacy XDMF+HDF5 mesh into dolfinx MeshTags?

Legacy meshes and mesh tags are created like in this example:

    from dolfin import UnitSquareMesh, MeshFunction, CompiledSubDomain
    from dolfin import XDMFFile

    nx = 4
    mesh = UnitSquareMesh(nx, nx)
    boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 0)

    left = CompiledSubDomain("near(x[0], 0) && on_boundary")
    right = CompiledSubDomain("near(x[0], 1) && on_boundary")
    bottom = CompiledSubDomain("near(x[1], 0) && on_boundary")
    top = CompiledSubDomain("near(x[1], 1) && on_boundary")

    for i, sd in enumerate((left, right, bottom, top)):
        sd.mark(boundaries, i + 1)

    with XDMFFile("legacy_mesh.xdmf") as xf:
        xf.write(mesh)
        xf.write(boundaries)

(where boundaries is a MeshFunction).

I can read in the mesh with adios4dolfinx:

import dolfinx
from adios4dolfinx import (read_function_from_legacy_h5,
                           read_mesh_from_legacy_h5)
from mpi4py import MPI


mesh = read_mesh_from_legacy_h5(
        "legacy_mesh.h5",
        comm=MPI.COMM_WORLD,
        group="/Mesh/mesh",
    )

But how can I read in the mesh tags?

Thanks!

I think I could make this work if there is a single MeshFunction in the file. As you have currently written it (with mesh and boundaries in the same file), the XDMFFile is not really well defined, as it looks like:

<?xml version="1.0"?>
<!DOCTYPE Xdmf SYSTEM "Xdmf.dtd" []>
<Xdmf Version="3.0" xmlns:xi="http://www.w3.org/2001/XInclude">
  <Domain>
    <Grid Name="mesh" GridType="Uniform">
      <Topology NumberOfElements="32" TopologyType="Triangle" NodesPerElement="3">
        <DataItem Dimensions="32 3" NumberType="UInt" Format="HDF">legacy_mesh.h5:/Mesh/mesh/topology</DataItem>
      </Topology>
      <Geometry GeometryType="XY">
        <DataItem Dimensions="25 2" Format="HDF">legacy_mesh.h5:/Mesh/mesh/geometry</DataItem>
      </Geometry>
    </Grid>
    <Grid Name="mesh" GridType="Uniform">
      <Topology NumberOfElements="56" TopologyType="PolyLine" NodesPerElement="2">
        <DataItem Dimensions="56 2" NumberType="UInt" Format="HDF">legacy_mesh.h5:/MeshFunction/0/mesh/topology</DataItem>
      </Topology>
      <Geometry Reference="XML">/Xdmf/Domain/Grid/Geometry</Geometry>
      <Attribute Name="f" AttributeType="Scalar" Center="Cell">
        <DataItem Dimensions="56 1" Format="HDF">legacy_mesh.h5:/MeshFunction/0/values</DataItem>
      </Attribute>
    </Grid>
  </Domain>
</Xdmf>

i.e. you have two mesh with the same name.
This will be super problematic when one wants to look up the data.

Thanks. I think thanks to the recent release of the X version of GitHub - finsberg/fenicsx-ldrb I don’t have to deal with legacy data anymore

1 Like