How to read functions saved as bp file with adios4dolfinx and solved using mixed_element

,

I am using dolfinx V0.8.0.
Here’s the MWE modified from demo_stokes.py

from mpi4py import MPI
from petsc4py import PETSc

import numpy as np

import ufl
from basix.ufl import element, mixed_element
from dolfinx import fem, la
from dolfinx.fem import (
    Function,
    dirichletbc,
    form,
    functionspace,
    locate_dofs_topological,
)
from dolfinx.fem.petsc import assemble_matrix_block, assemble_vector_block
from dolfinx.io import XDMFFile
from dolfinx.mesh import CellType, create_rectangle, locate_entities_boundary
from ufl import div, dx, grad, inner

# We create a {py:class}`Mesh <dolfinx.mesh.Mesh>`, define functions for
# locating geometrically subsets of the boundary, and define a function
# for the  velocity on the lid:

# +
# Create mesh
msh = create_rectangle(
    MPI.COMM_WORLD, [np.array([0, 0]), np.array([1, 1])], [32, 32], CellType.triangle
)


# Function to mark x = 0, x = 1 and y = 0
def noslip_boundary(x):
    return np.isclose(x[0], 0.0) | np.isclose(x[0], 1.0) | np.isclose(x[1], 0.0)


# Function to mark the lid (y = 1)
def lid(x):
    return np.isclose(x[1], 1.0)


# Lid velocity
def lid_velocity_expression(x):
    return np.stack((np.ones(x.shape[1]), np.zeros(x.shape[1])))


# -

# Two {py:class}`function spaces <dolfinx.fem.FunctionSpace>` are
# defined using different finite elements. `P2` corresponds to a
# continuous piecewise quadratic basis (vector) and `P1` to a continuous
# piecewise linear basis (scalar).


P2 = element("Lagrange", msh.basix_cell(), 2, shape=(msh.geometry.dim,))
P1 = element("Lagrange", msh.basix_cell(), 1)

def mixed_direct():
    # Create the Taylot-Hood function space
    TH = mixed_element([P2, P1])
    W = functionspace(msh, TH)

    # No slip boundary condition
    W0 = W.sub(0)
    Q, _ = W0.collapse()
    noslip = Function(Q)
    facets = locate_entities_boundary(msh, 1, noslip_boundary)
    dofs = locate_dofs_topological((W0, Q), 1, facets)
    bc0 = dirichletbc(noslip, dofs, W0)

    # Driving velocity condition u = (1, 0) on top boundary (y = 1)
    lid_velocity = Function(Q)
    lid_velocity.interpolate(lid_velocity_expression)
    facets = locate_entities_boundary(msh, 1, lid)
    dofs = locate_dofs_topological((W0, Q), 1, facets)
    bc1 = dirichletbc(lid_velocity, dofs, W0)

    # Collect Dirichlet boundary conditions
    bcs = [bc0, bc1]

    # Define variational problem
    (u, p) = ufl.TrialFunctions(W)
    (v, q) = ufl.TestFunctions(W)
    f = Function(Q)
    a = form((inner(grad(u), grad(v)) + inner(p, div(v)) + inner(div(u), q)) * dx)
    L = form(inner(f, v) * dx)

    # Assemble LHS matrix and RHS vector
    A = fem.petsc.assemble_matrix(a, bcs=bcs)
    A.assemble()
    b = fem.petsc.assemble_vector(L)

    fem.petsc.apply_lifting(b, [a], bcs=[bcs])
    b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)

    # Set Dirichlet boundary condition values in the RHS
    fem.petsc.set_bc(b, bcs)

    # Create and configure solver
    ksp = PETSc.KSP().create(msh.comm)
    ksp.setOperators(A)
    ksp.setType("preonly")

    # Configure MUMPS to handle pressure nullspace
    pc = ksp.getPC()
    pc.setType("lu")
    pc.setFactorSolverType("mumps")
    pc.setFactorSetUpSolverType()
    pc.getFactorMatrix().setMumpsIcntl(icntl=24, ival=1)
    pc.getFactorMatrix().setMumpsIcntl(icntl=25, ival=0)

    # Compute the solution
    U = Function(W)
    try:
        ksp.solve(b, U.x.petsc_vec)
    except PETSc.Error as e:
        if e.ierr == 92:
            print("The required PETSc solver/preconditioner is not available. Exiting.")
            print(e)
            exit(0)
        else:
            raise e

    # Split the mixed solution and collapse
    u, p = U.sub(0).collapse(), U.sub(1).collapse()

    return u,p

import adios4dolfinx
from pathlib import Path
u,p = mixed_direct()
filename = Path(f"function_checkpoint.bp")
adios4dolfinx.write_mesh(filename, msh)
adios4dolfinx.write_function(filename, u, time=0.0, name="velocity")
adios4dolfinx.write_function(filename, p, time=0.0, name="pressure")

My main purpose is being able to save the functions (velocity and pressure) as files and read them back in later so I can use u.eval().
I have a few questions:

  1. if I am using only one process, will create_rectangle always give me the same mesh given that the input parameters remain the same. In other words, is it worth it to save the mesh in bp.
  2. I am not sure if I am saving the functions correctly. (I intend to use the u,p functions later for interpolation).
  3. I don’t how I am supposed to read the functions as they are solved in mixed_element. Can I just construct the corresponding functionspace for velocity and pressure with P2 and P1 elements or do I still need the mixed_element in order to read them back in.

Thanks in advance!

Consider the following demo: Checkpoint on input mesh — ADIOS2Wrappers
As DOLFINx doesn’t treat built in meshes in a specialized way, there is internal re-ordering of cells and vertices due to data-locality, and one would have to use the following specific routine to do a checkpoint on the original mesh (i.e. the one made through create_rectangle).

See reply to 1.

I would collapse each function, as you do, and then just read in on a P2 space for u and P1 space for p

Thanks for your reply @dokken ! I have updated the MWE accordingly by saving the mesh as xdmf and write functions using write_function_on_input_mesh ,but it stills does not give the expected result:

from mpi4py import MPI
from petsc4py import PETSc

import numpy as np

import ufl
from basix.ufl import element, mixed_element
from dolfinx import fem, la
from dolfinx.fem import (
    Function,
    dirichletbc,
    form,
    functionspace,
    locate_dofs_topological,
)
from dolfinx.fem.petsc import assemble_matrix_block, assemble_vector_block
from dolfinx.io import XDMFFile
from dolfinx.mesh import CellType, create_rectangle, locate_entities_boundary
from ufl import div, dx, grad, inner

# Create mesh
msh = create_rectangle(
    MPI.COMM_WORLD, [np.array([0, 0]), np.array([1, 1])], [32, 32], CellType.triangle
)
# Save mesh
with XDMFFile(MPI.COMM_WORLD, "TestMesh.xdmf", "w") as xdmf:
        xdmf.write_mesh(msh)
# Function to mark x = 0, x = 1 and y = 0
def noslip_boundary(x):
    return np.isclose(x[0], 0.0) | np.isclose(x[0], 1.0) | np.isclose(x[1], 0.0)


# Function to mark the lid (y = 1)
def lid(x):
    return np.isclose(x[1], 1.0)


# Lid velocity
def lid_velocity_expression(x):
    return np.stack((np.ones(x.shape[1]), np.zeros(x.shape[1])))


# -

# Two {py:class}`function spaces <dolfinx.fem.FunctionSpace>` are
# defined using different finite elements. `P2` corresponds to a
# continuous piecewise quadratic basis (vector) and `P1` to a continuous
# piecewise linear basis (scalar).


P2 = element("Lagrange", msh.basix_cell(), 2, shape=(msh.geometry.dim,))
P1 = element("Lagrange", msh.basix_cell(), 1)

def mixed_direct():
    # Create the Taylot-Hood function space
    TH = mixed_element([P2, P1])
    W = functionspace(msh, TH)

    # No slip boundary condition
    W0 = W.sub(0)
    Q, _ = W0.collapse()
    noslip = Function(Q)
    facets = locate_entities_boundary(msh, 1, noslip_boundary)
    dofs = locate_dofs_topological((W0, Q), 1, facets)
    bc0 = dirichletbc(noslip, dofs, W0)

    # Driving velocity condition u = (1, 0) on top boundary (y = 1)
    lid_velocity = Function(Q)
    lid_velocity.interpolate(lid_velocity_expression)
    facets = locate_entities_boundary(msh, 1, lid)
    dofs = locate_dofs_topological((W0, Q), 1, facets)
    bc1 = dirichletbc(lid_velocity, dofs, W0)

    # Collect Dirichlet boundary conditions
    bcs = [bc0, bc1]

    # Define variational problem
    (u, p) = ufl.TrialFunctions(W)
    (v, q) = ufl.TestFunctions(W)
    f = Function(Q)
    a = form((inner(grad(u), grad(v)) + inner(p, div(v)) + inner(div(u), q)) * dx)
    L = form(inner(f, v) * dx)

    # Assemble LHS matrix and RHS vector
    A = fem.petsc.assemble_matrix(a, bcs=bcs)
    A.assemble()
    b = fem.petsc.assemble_vector(L)

    fem.petsc.apply_lifting(b, [a], bcs=[bcs])
    b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)

    # Set Dirichlet boundary condition values in the RHS
    fem.petsc.set_bc(b, bcs)

    # Create and configure solver
    ksp = PETSc.KSP().create(msh.comm)
    ksp.setOperators(A)
    ksp.setType("preonly")

    # Configure MUMPS to handle pressure nullspace
    pc = ksp.getPC()
    pc.setType("lu")
    pc.setFactorSolverType("mumps")
    pc.setFactorSetUpSolverType()
    pc.getFactorMatrix().setMumpsIcntl(icntl=24, ival=1)
    pc.getFactorMatrix().setMumpsIcntl(icntl=25, ival=0)

    # Compute the solution
    U = Function(W)
    try:
        ksp.solve(b, U.x.petsc_vec)
    except PETSc.Error as e:
        if e.ierr == 92:
            print("The required PETSc solver/preconditioner is not available. Exiting.")
            print(e)
            exit(0)
        else:
            raise e

    # Split the mixed solution and collapse
    u, p = U.sub(0).collapse(), U.sub(1).collapse()

    return u,p

u,p = mixed_direct()

# save solutions
import adios4dolfinx

adios4dolfinx.write_function_on_input_mesh(
        "TestFunction.bp",
        u,
        mode=adios4dolfinx.adios2_helpers.adios2.Mode.Write,
        time=0.0,
        name="velocity",
    )


# read solutions
timestamp = 0
with XDMFFile(MPI.COMM_WORLD, "TestMesh.xdmf", "r") as xdmf:
    in_mesh = xdmf.read_mesh()
# Define function spaces based on the mesh
# P2 element for velocity (quadratic)
P2 = element("Lagrange", in_mesh.basix_cell(), 2, shape=(in_mesh.geometry.dim,))
# P1 element for pressure (linear)
P1 = element("Lagrange", in_mesh.basix_cell(), 1)

V = functionspace(in_mesh,P2)
velocity = Function(V)
adios4dolfinx.read_function("TestFunction.bp", velocity, time=timestamp, name="velocity")
np.testing.assert_allclose(u.x.array, velocity.x.array, atol=1e-14)

which gives

AssertionError: 
Not equal to tolerance rtol=1e-07, atol=1e-14

Mismatched elements: 8407 / 8450 (99.5%)
Max absolute difference: 1.19787861
Max relative difference: 6017470.16147278
 x: array([0., 0., 0., ..., 0., 0., 0.])
 y: array([ 0.      ,  0.      ,  0.      , ...,  0.      ,  0.276687,
       -0.09824 ])

Could you spot where the problem is?

But now you are not using the “input mesh”, as you are reading the partitioned and re-ordered mesh from an XDMF file:

You should use
in_msh = create_rectangle( MPI.COMM_WORLD, [np.array([0, 0]), np.array([1, 1])], [32, 32], CellType.triangle or store the mesh along with the checkpoint, as shown in the other demos at: Writing a function checkpoint — ADIOS2Wrappers

Thanks for the input! My apologies I still cannot get it working. Following Writing a function checkpoint — ADIOS2Wrappers was the first thing I attempted before sending the post and didn’t work, here’s the MWE:

from mpi4py import MPI
from petsc4py import PETSc

import numpy as np

import ufl
from basix.ufl import element, mixed_element
from dolfinx import fem, la
from dolfinx.fem import (
    Function,
    dirichletbc,
    form,
    functionspace,
    locate_dofs_topological,
)
from dolfinx.fem.petsc import assemble_matrix_block, assemble_vector_block
from dolfinx.io import XDMFFile
from dolfinx.mesh import CellType, create_rectangle, locate_entities_boundary
from ufl import div, dx, grad, inner

# We create a {py:class}`Mesh <dolfinx.mesh.Mesh>`, define functions for
# locating geometrically subsets of the boundary, and define a function
# for the  velocity on the lid:

# +
# Create mesh
msh = create_rectangle(
    MPI.COMM_WORLD, [np.array([0, 0]), np.array([1, 1])], [32, 32], CellType.triangle
)
# Function to mark x = 0, x = 1 and y = 0
def noslip_boundary(x):
    return np.isclose(x[0], 0.0) | np.isclose(x[0], 1.0) | np.isclose(x[1], 0.0)


# Function to mark the lid (y = 1)
def lid(x):
    return np.isclose(x[1], 1.0)


# Lid velocity
def lid_velocity_expression(x):
    return np.stack((np.ones(x.shape[1]), np.zeros(x.shape[1])))


# -

# Two {py:class}`function spaces <dolfinx.fem.FunctionSpace>` are
# defined using different finite elements. `P2` corresponds to a
# continuous piecewise quadratic basis (vector) and `P1` to a continuous
# piecewise linear basis (scalar).


P2 = element("Lagrange", msh.basix_cell(), 2, shape=(msh.geometry.dim,))
P1 = element("Lagrange", msh.basix_cell(), 1)

def mixed_direct():
    # Create the Taylot-Hood function space
    TH = mixed_element([P2, P1])
    W = functionspace(msh, TH)

    # No slip boundary condition
    W0 = W.sub(0)
    Q, _ = W0.collapse()
    noslip = Function(Q)
    facets = locate_entities_boundary(msh, 1, noslip_boundary)
    dofs = locate_dofs_topological((W0, Q), 1, facets)
    bc0 = dirichletbc(noslip, dofs, W0)

    # Driving velocity condition u = (1, 0) on top boundary (y = 1)
    lid_velocity = Function(Q)
    lid_velocity.interpolate(lid_velocity_expression)
    facets = locate_entities_boundary(msh, 1, lid)
    dofs = locate_dofs_topological((W0, Q), 1, facets)
    bc1 = dirichletbc(lid_velocity, dofs, W0)

    # Collect Dirichlet boundary conditions
    bcs = [bc0, bc1]

    # Define variational problem
    (u, p) = ufl.TrialFunctions(W)
    (v, q) = ufl.TestFunctions(W)
    f = Function(Q)
    a = form((inner(grad(u), grad(v)) + inner(p, div(v)) + inner(div(u), q)) * dx)
    L = form(inner(f, v) * dx)

    # Assemble LHS matrix and RHS vector
    A = fem.petsc.assemble_matrix(a, bcs=bcs)
    A.assemble()
    b = fem.petsc.assemble_vector(L)

    fem.petsc.apply_lifting(b, [a], bcs=[bcs])
    b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)

    # Set Dirichlet boundary condition values in the RHS
    fem.petsc.set_bc(b, bcs)

    # Create and configure solver
    ksp = PETSc.KSP().create(msh.comm)
    ksp.setOperators(A)
    ksp.setType("preonly")

    # Configure MUMPS to handle pressure nullspace
    pc = ksp.getPC()
    pc.setType("lu")
    pc.setFactorSolverType("mumps")
    pc.setFactorSetUpSolverType()
    pc.getFactorMatrix().setMumpsIcntl(icntl=24, ival=1)
    pc.getFactorMatrix().setMumpsIcntl(icntl=25, ival=0)

    # Compute the solution
    U = Function(W)
    try:
        ksp.solve(b, U.x.petsc_vec)
    except PETSc.Error as e:
        if e.ierr == 92:
            print("The required PETSc solver/preconditioner is not available. Exiting.")
            print(e)
            exit(0)
        else:
            raise e

    # Split the mixed solution and collapse
    u, p = U.sub(0).collapse(), U.sub(1).collapse()

    return u,p

u,p = mixed_direct()

# save solutions
import adios4dolfinx
from pathlib import Path
# save mesh and function
filename = Path("function_checkpoint.bp")
adios4dolfinx.write_mesh(filename, msh)
adios4dolfinx.write_function(filename, u, time=0.0, name="velocity")
# read solutions
timestamp = 0.0
in_mesh = adios4dolfinx.read_mesh(filename, MPI.COMM_WORLD)
# Define function spaces based on the mesh
# P2 element for velocity (quadratic)
P2 = element("Lagrange", in_mesh.basix_cell(), 2, shape=(in_mesh.geometry.dim,))
# P1 element for pressure (linear)
P1 = element("Lagrange", in_mesh.basix_cell(), 1)
V = functionspace(in_mesh,P2)
velocity = Function(V)
adios4dolfinx.read_function("function_checkpoint.bp", velocity, time=timestamp, name="velocity")
np.testing.assert_allclose(u.x.array, velocity.x.array, atol=1e-14)

that gives

AssertionError: 
Not equal to tolerance rtol=1e-07, atol=1e-14

Mismatched elements: 8277 / 8450 (98%)
Max absolute difference: 1.05305187
Max relative difference: 1.07032121e+08
 x: array([0., 0., 0., ..., 0., 0., 0.])
 y: array([0., 0., 0., ..., 0., 0., 0.])

Please read through the demos that I have sent to you. Here is a minimal reproducible example, that illustrates the two different strategies:

from mpi4py import MPI

import numpy as np

from basix.ufl import element, mixed_element
import dolfinx


def f(x):
    return x[0], x[1]

def g(x):
    return np.sin(x[1])


def mixed_direct(msh: dolfinx.mesh.Mesh):
    P2 = element("Lagrange", msh.basix_cell(), 2, shape=(msh.geometry.dim,))
    P1 = element("Lagrange", msh.basix_cell(), 1)

    # Create the Taylot-Hood function space
    TH = mixed_element([P2, P1])
    W = dolfinx.fem.functionspace(msh, TH)
    U = dolfinx.fem.Function(W)
    U.sub(0).interpolate(lambda x: f(x))
    U.sub(1).interpolate(lambda x: g(x))

    # Split the mixed solution and collapse
    u, p = U.sub(0).collapse(), U.sub(1).collapse()

    return u,p

if __name__ == "__main__":
    import adios4dolfinx
    from pathlib import Path

    msh0 = dolfinx.mesh.create_rectangle(
    MPI.COMM_WORLD, [np.array([0, 0]), np.array([1, 1])], [32, 32], dolfinx.mesh.CellType.triangle)

    u,p = mixed_direct(msh0)
    u.name = "velocity"
    # METHOD 1: Save mesh to file, then save function to file
    filename = Path("function_checkpoint.bp")
    timestamp = 0.0
    adios4dolfinx.write_mesh(filename, msh0)
    adios4dolfinx.write_function(filename, u, time=timestamp)

    # Mesh is read in then function is read in
    in_mesh = adios4dolfinx.read_mesh(filename, MPI.COMM_WORLD)
    P2 = element("Lagrange", in_mesh.basix_cell(), 2, shape=(in_mesh.geometry.dim,))
    V = dolfinx.fem.functionspace(in_mesh, P2)
    velocity = dolfinx.fem.Function(V)
    velocity.name = u.name
    adios4dolfinx.read_function(filename, velocity, time=timestamp)

    # Compare with exact solution on this current grid
    velocity_reference = dolfinx.fem.Function(V)
    velocity_reference.interpolate(lambda x: f(x))

    np.testing.assert_allclose(velocity_reference.x.array, velocity.x.array, atol=1e-14)
  
    # METHOD 2: Only save function to file
    new_file = Path("function_checkpoint2.bp")
    adios4dolfinx.write_function_on_input_mesh(
        new_file,
        u,
        mode=adios4dolfinx.adios2_helpers.adios2.Mode.Write,
        time=timestamp,
    )

    # Read in function on msh0
    v = dolfinx.fem.Function(u.function_space)
    adios4dolfinx.read_function(new_file, v, time=timestamp, name=u.name)
    np.testing.assert_allclose(u.x.array, v.x.array, atol=1e-14)

Thank you so much for the MWE! I think I might have expressed my request in a wrong way. All I wanted to achieve is to save the mesh and the collapsed solution u,p in the dolfinx function format(correct me if I am wrong) of result from my MWE:

u,p = mixed_direct()

in a file (.bp using adios4dolfinx), so that I can read them back in later and use in Making curve plot throughout the domain section of the tutorial, namely:

from dolfinx import geometry
bb_tree = geometry.bb_tree(domain, domain.topology.dim)
cells = []
points_on_proc = []
# Find cells whose bounding-box collide with the the points
cell_candidates = geometry.compute_collisions_points(bb_tree, points.T)
# Choose one of the cells that contains the point
colliding_cells = geometry.compute_colliding_cells(domain, cell_candidates, points.T)
for i, point in enumerate(points.T):
    if len(colliding_cells.links(i)) > 0:
        points_on_proc.append(point)
        cells.append(colliding_cells.links(i)[0])
points_on_proc = np.array(points_on_proc, dtype=np.float64)
u_values = uh.eval(points_on_proc, cells)
p_values = pressure.eval(points_on_proc, cells)

where domain is the msh in my case.

Because it seemed that I was doing exactly what you did for method1 except that you did:

# Compare with exact solution on this current grid
velocity_reference = dolfinx.fem.Function(V)
velocity_reference.interpolate(lambda x: f(x))

np.testing.assert_allclose(velocity_reference.x.array, velocity.x.array, atol=1e-14)

and when I change the last line to (intend to check the velocity read from file is equivalent to computed u)

np.testing.assert_allclose(u.x.array, velocity.x.array, atol=1e-14)

it throws a similar error. Probably this is a wrong way to compared them?
Apologies for such a long message!

as they are on different meshes, you cannot compare them in this way. That is the whole point of the two different methods (as I’ve now tried to explain 2-3 times).

I use velocity_reference to verify that the data has been read in in the correct way on the mesh.

Now I get what you mean, thanks a lot for you patience and again apologies for my misunderstanding! I thought:

msh = create_rectangle(
    MPI.COMM_WORLD, [np.array([0, 0]), np.array([1, 1])], [32, 32], CellType.triangle
)
adios4dolfinx.write_mesh(filename, msh)
in_mesh = adios4dolfinx.read_mesh(filename, MPI.COMM_WORLD)

are going to give me the same mesh so they are just different. So I was doing the correct thing (sort of) but verifying it in a wrong way.