Gradients and divergences in fenicsx

[Gianna]

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

a = form([[inner(grad(u), grad(v)) * dx, inner(p, div(v)) * dx],
          [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():
    b_sub.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)

# 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")
ksp.getPC().setFieldSplitType(PETSc.PC.CompositeType.ADDITIVE)

# 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)
div_u = ufl.grad(u)
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):
    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 = 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))

[dokken]

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

[nate]

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

Finite element method for the Navier-Stokes equation with local conservation properties

If you’re solving the incompressible Navier-Stokes equations, I’d recommend looking into divergence conforming Discontinuous Galerkin (DG) methods. These methods are non-conforming and are similar in some ways to finite volume methods. Upwinding can be incorporated naturally through the specification of the numerical flux. Provided the finite element spaces are chosen such that the divergence operator maps the velocity space into the pressure space, mass is conserved exactly. Exact mass conserva…

Consider the posts regarding solenoidal approximations:

How to create a divergence-free VectorFunctionSpace

Usually the divergence-free constraint is enforced in a discrete sense, i.e., with respect to some finite dimensional set of pressure test functions: (\nabla\cdot\mathbf{u}^h,q^h)_{L^2} = 0~~~\forall q^h\in\mathcal{Q}^h\text{ ,} where \mathcal{Q}^h is the discrete pressure test space. A velocity satisfying such a condition can be solved for alongside the corresponding pressure by using a mixed element, as in the Stokes demo here: https://bitbucket.org/fenics-project/dolfin/src/master/python…

Projection of div-free vector field on crossed mesh is not div-free

Hi all, I feel like this should have come up already, but I was not able to find anything through the search function. Here is my MWE (2019.1.0): from dolfin import * n = 16 mesh_L = UnitSquareMesh( n, n, 'left') mesh_C = UnitSquareMesh( n, n, 'crossed') V_L = VectorFunctionSpace(mesh_L, 'CG', 1) V_C = VectorFunctionSpace(mesh_C, 'CG', 1) ic = Expression( ('sin(x[1])', 'sin(x[0])'), degree=1) v_L = project(ic, V_L) v_C = project(ic, V_C) print("div v_L_L2 =", assemble(div(v_L)*div(v_L)*dx…

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

[Leopart, Geopart] HDG Stokes velocity interpolation

It looks like you’re observing suboptimal convergence of the Stokes system in a non-convex domain. This is a pretty well known result and is typically noticed by singularities in the pressure approximation. Consider my modification of your example. Let \Vert \mathbf{v} \Vert_\text{noslip} = \sqrt{\int_{\partial\Omega_\text{wall}} \mathbf{v}^2 \mathrm{d} s} such that we measure \Vert \mathbf{u}_h \Vert_\text{noslip} and \Vert \bar{\mathbf{u}}_h \Vert_\text{noslip}. You should see that \Vert \…

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

[Gianna]

@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!

[dokken]

@Gianna

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

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

@Gianna

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 −∇2B=0-\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 ∇⋅B=0\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!

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

[Gianna]

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.