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

I think the problem here is that you are not collapsing the sub-function into the sub space (split gives you a view of the function, but not an actual sub function). What I would do is:

w = Function(W)
V, V_to_W = W.sub(0).collapse()
u = w.sub(0).collapse()
v = Function(V)
import adios4dolfinx
adios4dolfinx.write_function(u, "my_file")
adios4dolfinx.read_function(v, "my_file")

w.x.array[V_to_W] = v.x.array

Thanks for you’re response, and I have no error now. But now I can’t understand why this lines are not working:

w = Function(W)
V, V_to_W = W.sub(0).collapse()
Q, Q_to_W = W.sub(1).collapse()

# modifying my solution
w.x.array[V_to_W[0]] = 5.0

# splitting it
u = w.sub(0).collapse()
p = w.sub(1).collapse()

# check that my modification occured (prints 5.0 as expected)
print(u.x.array[0])

# creating new solution
a = Function(W)
v = Function(V)
q = Function(Q)

import adios4dolfinx
adios4dolfinx.write_function(u, "velocity_file")
adios4dolfinx.write_function(p, "pressure_file")

adios4dolfinx.read_function(v, "velocity_file")
adios4dolfinx.read_function(q, "pressure_file")

# after reading, assign
a.x.array[V_to_W] = v.x.array
a.x.array[Q_to_W] = q.x.array

# prints 0 but, as I understand it (wrong apparently), should print 5.0
print(a.x.array[V_to_W[0]])

Adding those lines at the end of my mwe, I expect those lines to print 5.0 and 5.0, but they print 5.0 and 0.0.

I covered this earlier on in this post,

I Also covered the reasoning for this and the implementation of a simple checkpointing at FEniCS 23, see:
https://jsdokken.com/checkpointing-presentation/#/
For the presentation and

For the source code.

Hey, I was wondering if now it is possible to read mesh and data from the same xdmf file? I have a xdmf/h5 file, produced by some application, which contains a mesh, a scalar function and a vector field. I would like to import this data into my code in dolfinx, is it possible? Thanks :slight_smile:

This depends on how the data is structured in the XDMFFile (as there is no one standard). Please provide a minimal example of such a file.

Ok! You can access an example of such file here

it has been generated with lifex-fiber: https://zenodo.org/record/5810269

You can read your mesh with:

import dolfinx
from mpi4py import MPI

with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "idealized_LV_ref_0.xdmf", "r") as xdmf:
    mesh = xdmf.read_mesh(name="mesh", xpath="Xdmf/Domain/Grid")

Reading the scalar values and vector values associated with the nodes requires some more work.
I might not have time to get back to you before I go on holiday, so the earliest I might have a reply for you is end of next week.

If you want to explore this, I would start with: Wrong order of dofs after creating mesh from lists of points and cells - #2 by dokken

@rosilho I had some spare time today and made a “quick and dirty” solution to this relying on my adios4dolfinx library (GitHub - jorgensd/adios4dolfinx: Interface of ADIOS2 for DOLFINx) and h5py.

I tested it with the file you supplied (in serial and parallel), and got the expected results:

from pathlib import Path

import adios4dolfinx
import dolfinx
import numpy as np
from mpi4py import MPI

with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "idealized_LV_ref_0.xdmf", "r") as xdmf:
    mesh = xdmf.read_mesh(name="mesh", xpath="Xdmf/Domain/Grid")

import basix
import h5py


def read_mesh_data(file_path:Path, mesh: dolfinx.mesh.Mesh, data_path:str):
    assert file_path.is_file(), f"File {file_path} does not exist"
    infile = h5py.File(file_path, "r", driver="mpio", comm=mesh.comm)
    num_nodes_global = mesh.geometry.index_map().size_global
    assert data_path in infile.keys(), f"Data {data_path} does not exist"
    dataset = infile[data_path]
    shape = dataset.shape
    assert shape[0] == num_nodes_global, f"Got data of shape {shape}, expected {num_nodes_global, shape[1]}"
    dtype = dataset.dtype
    # Read data locally on each process
    local_input_range = adios4dolfinx.comm_helpers.compute_local_range(mesh.comm, num_nodes_global)    
    local_input_data = dataset[local_input_range[0]:local_input_range[1]]

    # Create appropriate function space (based on coordinate map)
    assert len(mesh.geometry.cmaps) == 1, "Mixed cell-type meshes not supported"
    element = basix.ufl.element(
        basix.ElementFamily.P,
        mesh.topology.cell_name(),
        mesh.geometry.cmaps[0].degree,
        mesh.geometry.cmaps[0].variant,
        shape=(shape[1],),
        gdim=mesh.geometry.dim)

    # Assumption: Same doflayout for geometry and function space, cannot test in python    
    V = dolfinx.fem.FunctionSpace(mesh, element)
    uh = dolfinx.fem.Function(V, name=data_path)
    # Assume that mesh is first order for now
    assert mesh.geometry.cmaps[0].degree == 1, "Only linear meshes supported"
    x_dofmap = mesh.geometry.dofmap
    igi = np.array(mesh.geometry.input_global_indices, dtype=np.int64)
    global_geom_input = igi[x_dofmap]
    global_geom_owner = adios4dolfinx.utils.index_owner(mesh.comm, global_geom_input.reshape(-1), num_nodes_global)
    for i in range(shape[1]):
        arr_i = adios4dolfinx.comm_helpers.send_dofs_and_recv_values(global_geom_input.reshape(-1), global_geom_owner, mesh.comm, local_input_data[:,i], local_input_range[0])
        dof_pos = x_dofmap.reshape(-1)*shape[1]+i
        uh.x.array[dof_pos] = arr_i
    infile.close()
    return uh
infile = Path("idealized_LV_ref_0.h5")


for data in ["f0", "n0", "phi", "psi", "s0"]:
    uh = read_mesh_data(infile, mesh, data)
    with dolfinx.io.XDMFFile(mesh.comm, f"{data}.xdmf", "w") as xdmf:
        xdmf.write_mesh(mesh)
        xdmf.write_function(uh)


2 Likes

Awesome! It works, thanks a lot!
I used the docker dolfinx:nightly build, since for :latest I got the following error:

AttributeError: ‘dolfinx.cpp.mesh.Geometry’ object has no attribute ‘cmaps’. Did you mean: ‘cmap’?

Thanks again :slight_smile:

latest is no longer a maintained tag:

nightly is built nightly, so that is up to date with main.

1 Like

Hello everyone, I started with Fenics recently and have been very impressed by its functionality, however am struggling with i/o of data between different source files.

As I understand the best way is currently to use adios4dolfinx. However, I have problem recovering the data written and read this way. I’m even trying to reproduce the example of @hermanmak above but am getting completely different results. For me only the cases 1 and 3 work and otherwise the results are scrambled, Based on the above posts, 3 should be scrambled and 0 should not. I guess I am doing something wrong but cannot figure out what.

Here is the MWE:

import os
import sys
import numpy as np
import ufl
import dolfinx
from dolfinx import fem, io
from dolfinx.fem.petsc import assemble_matrix
from dolfinx.fem import FunctionSpace, Function

from mpi4py import MPI

new_path = r’/usr/local/lib/python3/dist-packages’ #Need to add this path for access to adios2 which is required by adios4dolfinx
sys.path.append(new_path)

import adios4dolfinx as a4dx

comm = MPI.COMM_WORLD
mesh = dolfinx.mesh.create_unit_square(comm, 100,100, cell_type=dolfinx.mesh.CellType.triangle)
filename = f"./{comm.size}proc.bp"

write

v_cg2 = ufl.VectorElement(“CG”, mesh.ufl_cell(), 2)
V = FunctionSpace(mesh, v_cg2)
u = Function(V)
x = ufl.SpatialCoordinate(mesh)
dfnx_expr = fem.Expression(x, V.element.interpolation_points())
u.interpolate(dfnx_expr)
a4dx.write_mesh(mesh, “./mesh/square.bp”)
a4dx.write_function(u, filename)

read

mesh_new = a4dx.read_mesh(comm, “./mesh/square.bp”, “BP4”, dolfinx.mesh.GhostMode.shared_facet)
v_new = ufl.VectorElement(“CG”, mesh_new.ufl_cell(), 2)
W = FunctionSpace(mesh_new, v_new)
w = Function(W)
v = Function(V)
a4dx.read_function(w, filename)
a4dx.read_function(v, filename)
#v.x.array[:] = o.x.array # array sizes not compatible?
with io.XDMFFile(comm, f"./0_{comm.size}proc.xdmf", “w”) as xdmf:
xdmf.write_mesh(mesh_new)
xdmf.write_function(w)
with io.XDMFFile(comm, f"./1_{comm.size}proc.xdmf", “w”) as xdmf:
xdmf.write_mesh(mesh)
xdmf.write_function(u)
with io.XDMFFile(comm, f"./2_{comm.size}proc.xdmf", “w”) as xdmf:
xdmf.write_mesh(mesh)
xdmf.write_function(w)
with io.XDMFFile(comm, f"./3_{comm.size}proc.xdmf", “w”) as xdmf:
xdmf.write_mesh(mesh_new)
xdmf.write_function(u)
with io.XDMFFile(comm, f"./4_{comm.size}proc.xdmf", “w”) as xdmf:
xdmf.write_mesh(mesh)
xdmf.write_function(v)

I am running this in Jupyterlab from the stable Docker image with Adios4dolfinx v0.1.0. I also tried running it on the nightly image with v0.4.0 but there I have some issues with writing the .vtk files. Nevertheless at some point I somehow managed to plot it also there and the results were the same as on stable so I don’t think this is the version issue.

I would appreciate any help. Thank you very much in advance.

What are the issues you are having with writing vtk-files?

In your example, only 0 and 1 should work.
I can reproduce the expected behavior with the following script:



import os
import sys
import numpy as np
import ufl
import dolfinx
from dolfinx import fem, io
from dolfinx.fem.petsc import assemble_matrix
from dolfinx.fem import FunctionSpace, Function

from mpi4py import MPI

new_path = r"/usr/local/lib/python3/dist-packages" #Need to add this path for access to adios2 which is required by adios4dolfinx
sys.path.append(new_path)

import adios4dolfinx as a4dx

comm = MPI.COMM_WORLD
mesh = dolfinx.mesh.create_unit_square(comm, 100,100, cell_type=dolfinx.mesh.CellType.triangle)
filename = f"./{comm.size}proc.bp"

v_cg2 = ufl.VectorElement("CG", mesh.ufl_cell(), 2)
V = FunctionSpace(mesh, v_cg2)
u = Function(V)
x = ufl.SpatialCoordinate(mesh)
dfnx_expr = fem.Expression(x, V.element.interpolation_points())
u.interpolate(dfnx_expr)
a4dx.write_mesh(mesh, "./mesh/square.bp")
a4dx.write_function(u, filename)

mesh_new = a4dx.read_mesh(comm, "./mesh/square.bp", "BP4", dolfinx.mesh.GhostMode.shared_facet)
v_new = ufl.VectorElement("CG", mesh_new.ufl_cell(), 2)
W = FunctionSpace(mesh_new, v_new)
w = Function(W)

v = Function(V)
a4dx.read_function(w, filename)
a4dx.read_function(v, filename)
#v.x.array[:] = o.x.array # array sizes not compatible?
with io.XDMFFile(comm, f"./0_{comm.size}proc.xdmf", "w") as xdmf:
    # Should work
    xdmf.write_mesh(mesh_new)
    xdmf.write_function(w)

with io.XDMFFile(comm, f"./1_{comm.size}proc.xdmf", "w") as xdmf:
    # Should work
    xdmf.write_mesh(mesh)
    xdmf.write_function(u)

with io.XDMFFile(comm, f"./2_{comm.size}proc.xdmf", "w") as xdmf:
    # Should not work
    xdmf.write_mesh(mesh)
    xdmf.write_function(w)

with io.XDMFFile(comm, f"./3_{comm.size}proc.xdmf", "w") as xdmf:
    # Should not work
    xdmf.write_mesh(mesh_new)
    xdmf.write_function(u)

with io.XDMFFile(comm, f"./4_{comm.size}proc.xdmf", "w") as xdmf:
    # Should not work
    xdmf.write_mesh(mesh)
    xdmf.write_function(v)

with

docker run -ti -v $(pwd):/root/shared -w /root/shared --rm ghcr.io/fenics/dolfinx/dolfinx:v0.6.0-r1
python3 -m pip install -U pip setuptools
python3 -m pip install git+https://github.com/jorgensd/adios4dolfinx@v0.1.0
mpirun -n 3 python3 name_of_script.py 

Thank you, this helped me realize at least the part of the issue. It works for the real build, also in Jupyterlab, but it doesn’t work for the complex one. Also not if I reproduce your approach. Do you have any suggestion for that? Is there maybe any way to use the complex build but write the real and imaginary parts separately with a4dx?

Regarding the vtk files, I am getting the following error with the nightly build:


RuntimeError Traceback (most recent call last)
Cell In[1], line 42
39 with io.XDMFFile(comm, f"./0_{comm.size}proc.xdmf", “w”) as xdmf:
40 # Should work
41 xdmf.write_mesh(mesh_new)
—> 42 xdmf.write_function(w)
44 with io.XDMFFile(comm, f"./1_{comm.size}proc.xdmf", “w”) as xdmf:
45 # Should work
46 xdmf.write_mesh(mesh)

File /usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/io/utils.py:234, in XDMFFile.write_function(self, u, t, mesh_xpath)
218 def write_function(self, u: Function, t: float = 0.0, mesh_xpath=“/Xdmf/Domain/Grid[@GridType=‘Uniform’][1]”):
219 “”“Write function to file for a given time.
220
221 Note:
(…)
232
233 “””
→ 234 super().write_function(getattr(u, “_cpp_object”, u), t, mesh_xpath)

RuntimeError: Degree of output Function must be same as mesh degree. Maybe the Function needs to be interpolated?

This happens for both the real and complex build and it does not happen in stable.

This is because you are writing out a P-2 field on a P-1 grid.
In 0.6.0 XDMFFile implicitly interpolated into P-1.
This has been removed in later versions, and the user has to interpolate the function into an appropriate function space for usage with xdmf (P-1 in your case).

I would have to adapt the source code of adios4dolfinx to work with complex-valued functions. Ive made an issue: Complex valued support · Issue #30 · jorgensd/adios4dolfinx · GitHub

1 Like

That makes sense, thank you.

Thanks, is there any workaround I can do in do meantime?

I’ve made a pull request: Support writing complex valued functions by jorgensd · Pull Request #31 · jorgensd/adios4dolfinx · GitHub

It requires some more testing from my end before merging, but feel free to test it.

wow, I appreciate the quick fix!

However, right now I am getting this error:


ValueError Traceback (most recent call last)
Cell In[1], line 36
33 w = Function(W)
35 v = Function(V)
—> 36 a4dx.read_function(w, filename)
37 a4dx.read_function(v, filename)
38 #v.x.array[:] = o.x.array # array sizes not compatible?

File /usr/local/lib/python3.10/dist-packages/adios4dolfinx/checkpointing.py:190, in read_function(u, filename, engine)
188 dofmap_path = “Dofmap”
189 xdofmap_path = “XDofmap”
→ 190 input_dofmap = read_dofmap(
191 comm, filename, dofmap_path, xdofmap_path, num_cells_global, engine
192 )
193 # Compute owner of dofs in dofmap
194 num_dofs_global = (
195 u.function_space.dofmap.index_map.size_global
196 * u.function_space.dofmap.index_map_bs
197 )

File /usr/local/lib/python3.10/dist-packages/adios4dolfinx/adios2_helpers.py:118, in read_dofmap(comm, filename, dofmap, dofmap_offsets, num_cells_global, engine)
116 io = adios.DeclareIO(io_name)
117 io.SetEngine(engine)
→ 118 infile = io.Open(str(filename), adios2.Mode.Read)
120 # First find step with dofmap offsets, to be able to read in a full row of the dofmap
121 for i in range(infile.Steps()):

ValueError: [Tue Sep 19 07:44:55 2023] [ADIOS2 EXCEPTION] : variable Values already defined in IO dofmap=‘Dofmap’_reader

This happens even with v0.4.0 and even with v0.3.0 now but it was not happening before as far as I remember.

I’m not quite sure what you mean now, as you are talking about versions that does not have this fix in them.

Sorry, what I wanted to say is that it appears in the new fixed version and then I also tested it in all the other mentioned versions and it appears there too. But when I was working with it last week, this error was not occuring in those versions. Could it be some incompatibility with some recent updates in the nightly build and adios4dolfinx?