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

dolfinx
,
1

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!

2

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

3

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?

4

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

5

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.])
6

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)
7

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!

8

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.

9

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.