# Gradients and divergences in fenicsx

Hi all,

The following code is a very slightly modified code taken from the Stokes equations tutorial. I want to check how “well” the constrain \nabla \cdot \mathbf{u}=\partial _xu_x+\partial _yu_y=0 holds. I need it for a 3d example I am trying to solve…

How can I calculate this quantity in fenicsx? I tried it with no success I guess due to the type of the finite elements ufl.VectorElement("Lagrange", msh.ufl_cell(), 2)?

Any help will be greatly appreciated!!!

import numpy as np
import matplotlib.pyplot as plt
import ufl
import dolfinx
from dolfinx import cpp as _cpp
from dolfinx import fem
from dolfinx.fem import (Constant, Function, FunctionSpace, dirichletbc,
extract_function_spaces, form,
locate_dofs_geometrical, locate_dofs_topological)
from dolfinx.io import XDMFFile
from dolfinx.mesh import (CellType, GhostMode, create_rectangle,
locate_entities_boundary)
from ufl import div, dx, grad, inner

from mpi4py import MPI
from petsc4py import PETSc
import math

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

# Generate experimental data
mesh=msh
# Create the function space
P2 = ufl.VectorElement("Lagrange", msh.ufl_cell(), 2)
P1 = ufl.FiniteElement("Lagrange", msh.ufl_cell(), 1)
V, Q = FunctionSpace(msh, P2), FunctionSpace(msh, P1)

mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
boundary_facets = dolfinx.mesh.exterior_facet_indices(mesh.topology)
boundary_dofs = dolfinx.fem.locate_dofs_topological(V, mesh.topology.dim-1, boundary_facets)
dof_coordinates = V.tabulate_dof_coordinates()[boundary_dofs]
points1 = mesh.geometry.x
# Export the mesh points and the boundary points
length=len(points1)
length_BC=len(boundary_dofs)
x_points=np.zeros(length)
y_points=np.zeros(length)
z_points=np.zeros(length)
x_points_boundary=np.zeros(length_BC)
y_points_boundary=np.zeros(length_BC)
z_points_boundary=np.zeros(length_BC)
x_points = [points1[i][0] for i in range(length)]
y_points = [points1[i][1] for i in range(length)]
z_points = [points1[i][2] for i in range(length)]
x_points_boundary = [dof_coordinates[i][0] for i in range(length_BC)]
y_points_boundary = [dof_coordinates[i][1] for i in range(length_BC)]
z_points_boundary = [dof_coordinates[i][2] for i in range(length_BC)]
points_boundary=np.zeros((length_BC,3))
points_boundary=[(x_points_boundary[i],y_points_boundary[i],z_points_boundary[i]) for i in range(length_BC)]
points = np.zeros((2, len(points1)))
points=[[x_points[i],y_points[i],z_points[i]] for i in range(length)]

# Function to mark x = 0, x = 1 and y = 0
def noslip_boundary(x):
return np.logical_or(np.logical_or(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])))

# No-slip boundary condition for velocity field (V) on boundaries
# where x = 0, x = 1, and y = 0
noslip = np.zeros(msh.geometry.dim, dtype=PETSc.ScalarType)
facets = locate_entities_boundary(msh, 1, noslip_boundary)
bc0 = dirichletbc(noslip, locate_dofs_topological(V, 1, facets), V)

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

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

# Define variational problem
(u, p) = ufl.TrialFunction(V), ufl.TrialFunction(Q)
(v, q) = ufl.TestFunction(V), ufl.TestFunction(Q)
f = Constant(msh, (PETSc.ScalarType(0), PETSc.ScalarType(0)))

[inner(div(u), q) * dx, None]])
L = form([inner(f, v) * dx, inner(Constant(msh, PETSc.ScalarType(0)), q) * dx])

a_p11 = form(inner(p, q) * dx)
a_p = [[a[0][0], None],
[None, a_p11]]

########## SOLUTION ##########

A = fem.petsc.assemble_matrix_nest(a, bcs=bcs)
A.assemble()
P11 = fem.petsc.assemble_matrix(a_p11, [])
P = PETSc.Mat().createNest([[A.getNestSubMatrix(0, 0), None], [None, P11]])
P.assemble()
b = fem.petsc.assemble_vector_nest(L)

# Modify ('lift') the RHS for Dirichlet boundary conditions
fem.petsc.apply_lifting_nest(b, a, bcs=bcs)

# Sum contributions from ghost entries on the owner
for b_sub in b.getNestSubVecs():

# Set Dirichlet boundary condition values in the RHS
bcs0 = fem.bcs_by_block(extract_function_spaces(L), bcs)
fem.petsc.set_bc_nest(b, bcs0)

# Create nullspace vector
null_vec = fem.petsc.create_vector_nest(L)

# Set velocity part to zero and the pressure part to a non-zero constant
null_vecs = null_vec.getNestSubVecs()
null_vecs[0].set(0.0), null_vecs[1].set(1.0)

# Normalize the vector, create a nullspace object, and attach it to the
# matrix
null_vec.normalize()
nsp = PETSc.NullSpace().create(vectors=[null_vec])
assert nsp.test(A)
A.setNullSpace(nsp)

ksp = PETSc.KSP().create(msh.comm)
ksp.setOperators(A, P)
ksp.setType("minres")
ksp.setTolerances(rtol=1e-9)
ksp.getPC().setType("fieldsplit")

# Define the matrix blocks in the preconditioner with the velocity and
# pressure matrix index sets
nested_IS = P.getNestISs()
ksp.getPC().setFieldSplitIS(
("u", nested_IS[0][0]),
("p", nested_IS[0][1]))

# Set the preconditioners for each block
ksp_u, ksp_p = ksp.getPC().getFieldSplitSubKSP()
ksp_u.setType("preonly")
ksp_u.getPC().setType("gamg")
ksp_p.setType("preonly")
ksp_p.getPC().setType("jacobi")

# Monitor the convergence of the KSP
ksp.setFromOptions()

u, p = Function(V), Function(Q)
x = PETSc.Vec().createNest([_cpp.la.petsc.create_vector_wrap(u.x), _cpp.la.petsc.create_vector_wrap(p.x)])
ksp.solve(b, x)
print(np.shape(points))

u_values = []
p_values = []
div_u_values = []
from dolfinx import geometry
bb_tree = geometry.BoundingBoxTree(msh, msh.topology.dim)
cells = []
points_on_proc = []
# Find cells whose bounding-box collide with the the points
cell_candidates = geometry.compute_collisions(bb_tree, points)
# Choose one of the cells that contains the point
colliding_cells = geometry.compute_colliding_cells(msh, cell_candidates, points)
for i, point in enumerate(points):
points_on_proc.append(point)
points_on_proc = np.array(points_on_proc, dtype=np.float64)
u_values = u.eval(points_on_proc, cells)
p_values = p.eval(points_on_proc, cells)
#div_u_values = div_u.eval(points_on_proc, cells)

uu=np.array(u_values)
p1=np.array(p_values)

ux=np.zeros(length)
uy=np.zeros(length)
pp=np.zeros(length)
ux=[uu[i][0] for i in range(length)]
uy=[uu[i][1] for i in range(length)]
pp=[p1[i]    for i in range(length)]

norm_u_0 = u.x.norm()
norm_p_0 = p.x.norm()
if MPI.COMM_WORLD.rank == 0:
print("(A) Norm of velocity coefficient vector (nested, iterative): {}".format(norm_u_0))
print("(A) Norm of pressure coefficient vector (nested, iterative): {}".format(norm_p_0))

You can:

• Interpolate ufl.div(u) into an appropriate function space using dolfinx.fem.Expression.
• Project it into a suitable function space.
• Assemble the error of the diverence over any region of the mesh.

See Electromagnetics example — FEniCSx tutorial this for an example where this is done for the curl operator

2 Likes

Additional to @dokken’s comment, this is an area of active research and you should be careful about how you measure the divergence of the velocity approximation.

@jpdean has a small collection of papers here

Consider the posts regarding solenoidal approximations:

Somewhat related regarding pressure robust approximations and their impact on the divergence free constraint of the approximation

And of course I’ll advertise mine and coworkers’ works and results with divergence free fields:
https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2020GC009349
https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2021GC009922

2 Likes

@dokken
Thanks for the suggestion, I somehow did not think to try the dolfinx.fem.Expression, I tried something like this

WW = FunctionSpace(mesh, ("DG", 2))
divuu = Function(WW)
u_expr = dolfinx.fem.Expression(ufl.div(u), WW.element.interpolation_points())
divuu.interpolate(u_expr)
divu=divuu.x.array

Plotting shows some points in the boundary where the divu is not zero though.

@nate Thanks for bringing those to my attention! To be honest I did not know how far the rabbit hole went.

Now my problem is slightly different. I want to calculate the Magnetic Field in some 3D geometry. Now there is some experimental input for the boundary. Long story short the boundaries have some sort of uncertainty in their values and by solving the 3D Poisson equation -\nabla ^2\mathbf{B}=0 I get a family of solutions and I need to keep the physical one, i.e. the one for which \nabla \cdot \mathbf{B}=0. Introducing a Lagrange multiplier to impose the constraint you are left with a system of PDEs that are similar to the Stokes ones, .i.e. the magnetic field is the velocity and the Lagrange multiplier is the pressure. The BC are different though, so since I used fenics a lot in the past I tried to solve my problem with fenicsx using the Poisson and stokes tutorial as a basis.

Seems like I stumbled in a non trivial problem, so thanks again for the help!

A Lagrange space is not divergence free, so this is expected. (Ref. @nate’s references)

Is there a question in this? Is something not working as expected, or are you getting an error message?

1 Like

No more questions, I just felt like giving some more context as to why imposing the Divergenceless condition is a priority for me. Using the code I posted works, even though, now that I see nate’s references I understand why the problem is difficult.

Thanks again, I may follow up with another question in the future and I hope this topic helps some more people.

1 Like