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:
- 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.
- I am not sure if I am saving the functions correctly. (I intend to use the u,p functions later for interpolation).
- 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!