Reading back the xdmf file

General
1

Hello,
In the code bellow, I create an XDMF file which saves the strain tensor in all of the time steps:

import dolfinx as fe
from dolfinx.fem.petsc import NonlinearProblem, LinearProblem, assemble_vector, assemble_matrix, create_vector, apply_lifting, set_bc
from dolfinx.nls.petsc import NewtonSolver
import numpy as np
import os
from mpi4py import MPI
import ufl
from petsc4py import PETSc


def main(linear = True, time_dependent = True):
    
    L = 1
    W = 0.2
    mesh = fe.mesh.create_box(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([L, W, W])],
                         [20, 6, 6], cell_type=fe.mesh.CellType.hexahedron)
    element = ufl.VectorElement("Lagrange", mesh.ufl_cell(), 1)
    V = fe.fem.functionspace(mesh, element)

    dt = fe.default_scalar_type(0.1)
    num_steps = int(1/float(dt))
    maximum_pressure = 15
    minimum_pressure = 0
    pressure = np.append(np.linspace(minimum_pressure, maximum_pressure, num=int(num_steps/2)),\
                         np.linspace(maximum_pressure, minimum_pressure, num=int(num_steps/2)))
    p = fe.fem.Constant(mesh, pressure[0])

    def clamped_boundary(x):
        return np.isclose(x[0], 0)

    fdim = mesh.topology.dim - 1
    clamped_boundary_facets = fe.mesh.locate_entities_boundary(mesh, fdim, clamped_boundary)
    u_D = np.array([0, 0, 0], dtype=fe.default_scalar_type)
    bc = fe.fem.dirichletbc(u_D, fe.fem.locate_dofs_topological(V, fdim, clamped_boundary_facets), V)

    boundaries = [(1, lambda x: np.isclose(x[0], L))]

    facet_indices, facet_markers = [], []
    fdim = mesh.topology.dim - 1
    for (marker, locator) in boundaries:
        facets = fe.mesh.locate_entities(mesh, fdim, locator)
        facet_indices.append(facets)
        facet_markers.append(np.full_like(facets, marker))
    facet_indices = np.hstack(facet_indices).astype(np.int32)
    facet_markers = np.hstack(facet_markers).astype(np.int32)
    sorted_facets = np.argsort(facet_indices)
    facet_tag = fe.mesh.meshtags(mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])


    ######################################################################################################
    #########################################Formulation##################################################
    u = ufl.TrialFunction(V)            
    v  = ufl.TestFunction(V)          
    old_u  = fe.fem.Function(V)
    old_velocity  = fe.fem.Function(V)
    old_acceleration  = fe.fem.Function(V)

    d = len(u)
    I = ufl.variable(ufl.Identity(d))             # Identity tensor
    F = ufl.variable(I + ufl.grad(u))             # Deformation gradient
    C = ufl.variable(F.T*F)                   # Right Cauchy-Green tensor
    J = ufl.variable(ufl.det(F))

    ####Check out for the metadata
    metadata = {"quadrature_degree": 4}
    ds = ufl.Measure("ds", domain=mesh, subdomain_data=facet_tag, metadata=metadata)
    dx = ufl.Measure("dx", domain=mesh, metadata=metadata)

    # Generalized-alpha method parameters
    alpha_m = fe.fem.Constant(mesh, 0.2)
    alpha_f = fe.fem.Constant(mesh, 0.4)
    gamma   = fe.default_scalar_type(0.5+alpha_f-alpha_m)
    beta    = fe.default_scalar_type((gamma+0.5)**2/4.)

    # Update formula for acceleration
    # a = 1/(2*beta)*((u - u0 - v0*dt)/(0.5*dt*dt) - (1-2*beta)*a0)
    def update_a(u, u_old, v_old, a_old, ufl=True):
        if ufl:
            dt_ = dt
            beta_ = beta
        else:
            dt_ = float(dt)
            beta_ = float(beta)
        return (u-u_old-dt_*v_old)/beta_/dt_**2 - (1-2*beta_)/2/beta_*a_old

    # Update formula for velocity
    # v = dt * ((1-gamma)*a0 + gamma*a) + v0
    def update_v(a, u_old, v_old, a_old, ufl=True):
        if ufl:
            dt_ = dt
            gamma_ = gamma
        else:
            dt_ = float(dt)
            gamma_ = float(gamma)
        return v_old + dt_*((1-gamma_)*a_old + gamma_*a)

    def update_fields(u_new, u_old, v_old, a_old):
        """Update fields at the end of each time step."""
        # Get vectors (references)
        u_vec, u0_vec  = u_new.x.array[:], u_old.x.array[:]
        v0_vec, a0_vec = v_old.x.array[:], a_old.x.array[:]

        # use update functions using vector arguments
        a_vec = update_a(u_vec, u0_vec, v0_vec, a0_vec, ufl=False)
        v_vec = update_v(a_vec, u0_vec, v0_vec, a0_vec, ufl=False)

        # Update (u_old <- u)
        v_old.x.array[:], a_old.x.array[:] = v_vec, a_vec
        u_old.x.array[:] = u_new.x.array[:]

    
    def avg(x_old, x_new, alpha):
        return alpha*x_old + (1-alpha)*x_new

    normal = -ufl.FacetNormal(mesh)
    
    rho = 1
    E = fe.default_scalar_type(1.0e4)
    nu = fe.default_scalar_type(0.3)
    mu = fe.fem.Constant(mesh, E / (2 * (1 + nu)))
    lmbda = fe.fem.Constant(mesh, E * nu / ((1 + nu) * (1 - 2 * nu)))
    
    def S(u):
        return 2.0 * mu * ufl.sym(ufl.grad(u)) + lmbda * ufl.tr(ufl.sym(ufl.grad(u))) * I

    acceleration = update_a(u, old_u, old_velocity, old_acceleration, ufl=False)
    velocity = update_v(acceleration, old_u, old_velocity, old_acceleration, ufl=True)

    formulation = rho * ufl.inner(alpha_m*old_acceleration + (1-alpha_m)* acceleration, v) * dx \
      + ufl.inner(ufl.grad(v), S(avg(old_u, u, alpha_f))) * dx \
      - ufl.inner(v, p * normal) * ds(1)

    bilinear_form = fe.fem.form(ufl.lhs(formulation))
    linear_form = fe.fem.form(ufl.rhs(formulation))

    ######################################################################################################
    ######################################################################################################

    ######################################################################################################
    ###############################################Solving################################################
    A = assemble_matrix(bilinear_form, bcs=[bc])
    A.assemble()
    b = create_vector(linear_form)

    solver = PETSc.KSP().create(mesh.comm)
    solver.setInitialGuessNonzero(True)
    solver.setOperators(A)
    solver.getPC().setType(PETSc.PC.Type.SOR)
    solver.view()

    displacement = fe.fem.Function(V)

    ###############Strain#################
    el_strain = ufl.TensorElement("Lagrange", mesh.ufl_cell(), 1, shape=(3, 3))
    Q_strain = fe.fem.FunctionSpace(mesh, el_strain)
    lagrange_finite_strain = ufl.dot((ufl.Identity(d) + ufl.grad(displacement)).T, ufl.Identity(d) + ufl.grad(displacement)) - ufl.Identity(d)
    strain_expr = fe.fem.Expression(lagrange_finite_strain, Q_strain.element.interpolation_points())
    strain = fe.fem.Function(Q_strain)     #### 2E = C - I
    strain.name = "lagrange_finite_strain"

    xdmf = fe.io.XDMFFile(mesh.comm, os.path.join("strain_minimal.xdmf"), "w")
    xdmf.write_mesh(mesh)

    for i in range(num_steps):
        p.value = pressure[i]

        # Update the right hand side reusing the initial vector
        with b.localForm() as loc_b:
            loc_b.set(0)
        assemble_vector(b, linear_form)

        # Apply Dirichlet boundary condition to the vector
        apply_lifting(b, [bilinear_form], [[bc]])
        b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
        set_bc(b, [bc])

        # Solve linear problem
        solver.solve(b, displacement.vector)
        displacement.x.scatter_forward()

        # Update old fields with new quantities
        update_fields(displacement, old_u, old_velocity, old_acceleration)

        strain.interpolate(strain_expr)
        xdmf.write_function(strain, i)

    xdmf.close()


if __name__ == "__main__":
    main()

And later, when I want to read back the data:

import os
import sys

import PetscBinaryIO

io = PetscBinaryIO.PetscBinaryIO(complexscalars=True)
f = io.readBinaryFile(os.path.join("strain_minimal.xdmf"))
print(f)

I get this error:
Traceback (most recent call last):
File “/Users/korosh/anaconda3/pkgs/petsc-3.20.6-real_he4ff8f5_100/lib/petsc/bin/PetscBinaryIO.py”, line 410, in readObjectType
objecttype = self._classid[header]
~~~~~~~~~~~~~^^^^^^^^
KeyError: 1010792557

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
File “/Users/korosh/Desktop/Korosh/cardiac_simulation/myo_simulation/minimal_example_reading_file.py”, line 11, in
f = io.readBinaryFile(os.path.join(“strain_minimal.xdmf”))
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File “/Users/korosh/anaconda3/pkgs/petsc-3.20.6-real_he4ff8f5_100/lib/petsc/bin/PetscBinaryIO.py”, line 109, in decorated_f
result = f(self, *args, **kwargs)
^^^^^^^^^^^^^^^^^^^^^^^^
File “/Users/korosh/anaconda3/pkgs/petsc-3.20.6-real_he4ff8f5_100/lib/petsc/bin/PetscBinaryIO.py”, line 443, in readBinaryFile
objecttype = self.readObjectType(fid)
^^^^^^^^^^^^^^^^^^^^^^^^
File “/Users/korosh/anaconda3/pkgs/petsc-3.20.6-real_he4ff8f5_100/lib/petsc/bin/PetscBinaryIO.py”, line 109, in decorated_f
result = f(self, *args, **kwargs)
^^^^^^^^^^^^^^^^^^^^^^^^
File “/Users/korosh/anaconda3/pkgs/petsc-3.20.6-real_he4ff8f5_100/lib/petsc/bin/PetscBinaryIO.py”, line 412, in readObjectType
raise IOError(‘Invalid PetscObject CLASSID or object not implemented for python’)
OSError: Invalid PetscObject CLASSID or object not implemented for python

Thanks in advance.

2

Not sure if PETSc really does support reading in a file generated by another software (in this case, dolfinx.io).

If your goal is to restart the simulation at a later stage, please consider GitHub - jorgensd/adios4dolfinx: Extending DOLFINx with checkpointing functionality instead.

If your goal is to dump the data in some tabular format, then use PETSc I/O to dump the PETSc vector stored inside the dolfinx.fem.Function, accessible with dolfinx.fem.Function.vector (dolfinx < 0.8.0) or dolfinx.fem.Function.petsc_vec (dolfinx >= 0.8.0).

3

Thank you for your attention. For the PETSc I/O case, do I need different files to write dolfinx.fem.Function corresponding to different time steps or is there a way to have all of them in a single file as it is the case with dolfinx.io?

4

You’ll have to look into PETSc documentation for this. There seem to be an append mode at PetscViewerFileSetMode — PETSc 3.21.2 documentation , but I’ve personally never used it (I typically save one vector per file).

1 Like