How are integral constraints handled in Dolfinx?

[fea]

I am on the latest version of Dolfinx. I have a basic 2D Poisson Equation on a unit square with homogeneous boundary conditions
\nabla^2 u = -f \ \ \ \ \forall (x,y) \in \Omega
u = 0 \ \ \ \ \forall (x,y) \in \partial\Omega,

and I want my solution to satisfy another constraint:
\int_\Omega u dx= c

where c is some constant. How would an integral constraint like this one be enforced in Dolfinx?

If it is not too much trouble, below is an example I would like to see this implemented for:

---

(I do not have an exact example grounded in physics, so here is one in which I set c to a purely arbitrary constant)
\nabla^2 u = -1
u = 0 \ \ \ \ \forall (x,y) \in \partial\Omega
\int_\Omega u = - 0.035

from mpi4py import MPI
import dolfinx as df
import ufl
from petsc4py import PETSc
from dolfinx.fem.petsc import LinearProblem

Nx,Ny = 10,10 # Number of points in x and y
Lx, Ly = 1,1
mesh = df.mesh.create_unit_square(MPI.COMM_WORLD,Nx, Ny)
tdim = mesh.topology.dim
fdim = tdim-1
mesh.topology.create_connectivity(tdim-1, tdim)

# functionspace (scalar space)
S = df.fem.functionspace(mesh, ("CG", 1)) 
# boundary conditions
exterior_facets = df.mesh.exterior_facet_indices(mesh.topology)
exterior_dofs = df.fem.locate_dofs_topological(S,fdim,exterior_facets)
bc = df.fem.dirichletbc(df.fem.Constant(mesh,df.default_scalar_type(0)),exterior_dofs,S)

# ---- Solve Poisson Equation ----
f = df.fem.Constant(mesh, df.default_scalar_type(-1))
u = ufl.TrialFunction(S)
v = ufl.TestFunction(S)
a = ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx
L = f * v * ufl.dx
problem = LinearProblem(
    a,
    L,
    bcs=[bc],
    petsc_options={"ksp_type": "preonly", "pc_type": "lu"},
    petsc_options_prefix="Poisson",
)
uh = problem.solve()

df.fem.assemble_scalar(df.fem.form(uh * ufl.dx)) returns -0.034029666046715173. I want to see how to enforce \int_\Omega u_h = -0.035 so that df.fem.assemble_scalar(df.fem.form(uh * ufl.dx)) will instead return -0.035.

[dokken]

See Real function spaces — scifem

[fea]

Unfortunately, this only half worked on my end. The integral constraint I set was satisfied, but the homogeneous Dirichlet boundary condition did not apply which can be seen based on observation of the plot of my solution:

where we can see that the boundary values are not 0. My code followed yours almost exactly, except I have created a Dirichlet boundary condition using Dolfinx to apply lifting, and the boundary integral was excluded from the weak form:

from mpi4py import MPI
import dolfinx as df
import ufl
from petsc4py import PETSc
import numpy as np
from scifem import create_real_functionspace
from scifem.petsc import apply_lifting_and_set_bc

M = 20
mesh = df.mesh.create_unit_square(
    MPI.COMM_WORLD, M, M, df.mesh.CellType.triangle, dtype=np.float64
)
V = df.fem.functionspace(mesh, ("Lagrange", 1))

#---- define boundary condition ----
fdim, tdim = mesh.topology.dim-1, mesh.topology.dim
mesh.topology.create_connectivity(fdim,tdim)
exterior_facets = df.mesh.exterior_facet_indices(mesh.topology)
exterior_dofs = df.fem.locate_dofs_topological(V,fdim,exterior_facets)
bc = df.fem.dirichletbc(df.fem.Constant(mesh,df.default_scalar_type(0)),exterior_dofs,V)

#---- Integral constraint ----
h = -0.035

# Forcing term
f = df.fem.Constant(mesh, df.default_scalar_type(-1))
# Real space
R = create_real_functionspace(mesh)
# Mixed space
W = ufl.MixedFunctionSpace(V, R)
u, lmbda = ufl.TrialFunctions(W)
du, dl = ufl.TestFunctions(W)

zero = df.fem.Constant(mesh, df.default_scalar_type(0.0))

#--------------------------------------------------
#---- Build Matrix --------------------------------
#--------------------------------------------------

# ---- Build components: bilinear form ----
a00 = ufl.inner(ufl.grad(u), ufl.grad(du)) * ufl.dx
a01 = ufl.inner(lmbda, du) * ufl.dx
a10 = ufl.inner(u, dl) * ufl.dx

# ---- Build components: linear form ----
L0 = ufl.inner(f, du) * ufl.dx
L1 = ufl.inner(zero, dl) * ufl.dx

# ---- Build matrix: A ----
a = [[a00, a01], [a10, None]]
a_compiled = df.fem.form(a)
A = df.fem.petsc.assemble_matrix(a_compiled)
A.assemble()

# ---- Build matrix: L ----
L = [L0, L1]
L_compiled = df.fem.form(L)

#--------------------------------------------------
#---- Apply boundary conditions -------------------
#--------------------------------------------------
bcs = [bc]
b = df.fem.petsc.assemble_vector(L_compiled, kind="mpi")
apply_lifting_and_set_bc(b, a_compiled, bcs=bcs)

#--------------------------------------------------
#---- Solve ---------------------------------------
#--------------------------------------------------
uh = df.fem.Function(V, name="u")

#---- Insert constraint into R ----
rh = df.fem.Function(R)
rh.x.array[0] = h
# Next we need to add this value to the existing right hand side vector.
# Therefore we create assign 0s to the primal space
b_real_space = b.duplicate()
uh.x.array[:] = 0
# Transfer the data to `b_real_space`
df.fem.petsc.assign([uh, rh], b_real_space)
# And accumulate the values in the right hand side vector
b.axpy(1, b_real_space)
# We destroy the temporary work vector after usage
b_real_space.destroy()

#---- KSP solver ----
ksp = PETSc.KSP().create(mesh.comm)
ksp.setOperators(A)
ksp.setType("preonly")
pc = ksp.getPC()
pc.setType("lu")
pc.setFactorSolverType("mumps")
xh = b.duplicate()

ksp.solve(b, xh)
xh.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)

uh = df.fem.Function(V, name="u")
df.fem.petsc.assign(xh, [uh, rh])

b.destroy()
xh.destroy()
A.destroy()
ksp.destroy()

When I applied the exact same homogeneous Dirichlet boundary condition to your scifem tutorial example, I notice that while it did change the solution, it however did not force the boundary values to 0. Am I using apply_lifting_and_set_bc incorrectly?

[spcvanschie]

You are over-constraining your finite element model. As is covered on this Wikipedia page, the Poisson equation with Dirichlet boundary conditions (i.e. your model minus the integral constraint) has a unique solution. Therefore you can’t just impose an extra integral constraint and hope to satisfy both that constraint and your original boundary conditions at the same time.
That same Wikipedia page mentions that the Poisson equation with Neumann boundary conditions has a unique solution up to a constant. So you should be able to impose the integral constraint at the same time as a set of Neumann boundary conditions. Which is exactly the case in the tutorial case that @dokken linked to.

[fea]

So if it has Dirichlet boundary conditions, then there is only one solution (ergo a unique solution) which can satisfy the problem i.e. there is no room for constraints? If that is the case, then this approach might not be the approach I was looking for for handling a problem which has been stonewalling me for a while now. The problem was

\nabla^2 \psi_I = \hat{\textbf{k}}\cdot \nabla \times \textbf{V} =curl(\textbf{V})
\psi_I = 0 \ \ \ \ \forall (x,y)\in \partial \Omega
where I was hoping I could add the constraint
\text{constraint}(\psi_I):=\oint_{\partial \Omega} \nabla \psi_I \cdot \hat{\textbf{s}} = -\oint_{\partial\Omega} \textbf{V}\cdot \hat{\textbf{n}}ds

The paper I am reading from stated that if this constraint was not satisfied, then modifications to the solution would need to be made to the solution until it did; my thought process was to modify the Poisson equation itself rather than the solution.

[spcvanschie]

So the problem setup that you’re talking about here is very different from your initial test problem. It seems like both \psi_I and \mathbf{V} are now unknowns? Since the integral constraint relates \psi_I to \mathbf{V}, the integral constraint doesn’t actually seem to constraint \psi_I, but rather acts as a constraint on \mathbf{V} (since again, the Poisson equation with Dirichlet boundary conditions leads to a unique solution for \psi_I). I don’t know whether \mathbf{V} is actually an unknown vector. I would assume that what’s meant in the paper that you mentioned, is that \mathbf{V} would need to be modified in order to satisfy the boundary integral constraint. If \mathbf{V} is also an unknown vector, you would likely need to supply some boundary conditions for it as well.

[fea]

\textbf{V} is known. It’s a Poisson Equation where the curl of the vector \textbf{V} is the forcing term where, again, \textbf{V} is known. Sorry if that was not explicitly clear! So the gist of the problem is that we want \psi_I to satisfy the Poisson equation and additionally an integral constraint such that the net circulation of the gradient of the solution along the boundary of the domain is the same (with opposite sign) as the net flux of the known vector \textbf{V}. The issue is that my solution was not satisfying this condition, so I tried to force a constraint into the weak form. The paper I am reading from recommended postprocessing it instead, and it is probably because of the reason you stated; the solution is unique so any additional constraints is over determining the problem.

[Stein]

@spcvanschie

the Poisson equation with Dirichlet boundary conditions (i.e. your model minus the integral constraint) has a unique solution. Therefore you can’t just impose an extra integral constraint and hope to satisfy both that constraint and your original boundary conditions at the same time.

This is certainly true, but often “adding a constraint” is understood to go hand-in-hand with the neccessity of incorporating a Lagrange multiplier that acts like the forcing to satisfy said constraint. Indeed, per your comment on solution uniqueness, the new solution (with constraint and Lagrange multiplier) does naturally not satisfy the PDE without the Lagrange multiplier forcing.

Applying that to @fea’s problem, you’d have to realize that what you’d actually be solving when you’re adding the constraint is:

\nabla^2 u = -1 + \lambda
u = 0 \ \ \ \ \forall (x,y) \in \partial\Omega
\int_\Omega u = - 0.035

With some constant \lambda.

If that is acceptable, then here follows the fixed version of the code. You were only missing the bcs in the assembly of the matrix:

from mpi4py import MPI
import dolfinx as df
import ufl
from petsc4py import PETSc
import numpy as np
from scifem import create_real_functionspace
from scifem.petsc import apply_lifting_and_set_bc

M = 20
mesh = df.mesh.create_unit_square(
    MPI.COMM_WORLD, M, M, df.mesh.CellType.triangle, dtype=np.float64
)
V = df.fem.functionspace(mesh, ("Lagrange", 1))

#---- define boundary condition ----
fdim, tdim = mesh.topology.dim-1, mesh.topology.dim
mesh.topology.create_connectivity(fdim,tdim)
exterior_facets = df.mesh.exterior_facet_indices(mesh.topology)
exterior_dofs = df.fem.locate_dofs_topological(V,fdim,exterior_facets)
bc = df.fem.dirichletbc(df.fem.Constant(mesh,df.default_scalar_type(0)),exterior_dofs,V)

#---- Integral constraint ----
h = -0.035

# Forcing term
f = df.fem.Constant(mesh, df.default_scalar_type(-1))
# Real space
R = create_real_functionspace(mesh)
# Mixed space
W = ufl.MixedFunctionSpace(V, R)
u, lmbda = ufl.TrialFunctions(W)
du, dl = ufl.TestFunctions(W)

zero = df.fem.Constant(mesh, df.default_scalar_type(0.0))

#--------------------------------------------------
#---- Build Matrix --------------------------------
#--------------------------------------------------

# ---- Build components: bilinear form ----
a00 = ufl.inner(ufl.grad(u), ufl.grad(du)) * ufl.dx
a01 = ufl.inner(lmbda, du) * ufl.dx
a10 = ufl.inner(u, dl) * ufl.dx

# ---- Build components: linear form ----
L0 = ufl.inner(f, du) * ufl.dx
L1 = ufl.inner(zero, dl) * ufl.dx

# ---- Build matrix: A ----

bcs = [bc]
a = [[a00, a01], [a10, None]]
a_compiled = df.fem.form(a)
A = df.fem.petsc.assemble_matrix(a_compiled, bcs=bcs)
A.assemble()

# ---- Build matrix: L ----
L = [L0, L1]
L_compiled = df.fem.form(L)

#--------------------------------------------------
#---- Apply boundary conditions -------------------
#--------------------------------------------------
b = df.fem.petsc.assemble_vector(L_compiled, kind="mpi")
apply_lifting_and_set_bc(b, a_compiled, bcs=bcs)

#--------------------------------------------------
#---- Solve ---------------------------------------
#--------------------------------------------------
uh = df.fem.Function(V, name="u")

#---- Insert constraint into R ----
rh = df.fem.Function(R)
rh.x.array[0] = h
# Next we need to add this value to the existing right hand side vector.
# Therefore we create assign 0s to the primal space
b_real_space = b.duplicate()
uh.x.array[:] = 0
# Transfer the data to `b_real_space`
df.fem.petsc.assign([uh, rh], b_real_space)
# And accumulate the values in the right hand side vector
b.axpy(1, b_real_space)
# We destroy the temporary work vector after usage
b_real_space.destroy()

#---- KSP solver ----
ksp = PETSc.KSP().create(mesh.comm)
ksp.setOperators(A)
ksp.setType("preonly")
pc = ksp.getPC()
pc.setType("lu")
pc.setFactorSolverType("mumps")
xh = b.duplicate()

ksp.solve(b, xh)
xh.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)

uh = df.fem.Function(V, name="u")
df.fem.petsc.assign(xh, [uh, rh])

import pyvista4dolfinx as p4d
p4d.plot(uh, show=True)

b.destroy()
xh.destroy()
A.destroy()
ksp.destroy()

I’ve added

import pyvista4dolfinx as p4d
p4d.plot(uh, show=True)

For plotting, yielding:

[fea]

Oh so to make it work, I have to add a correction onto the forcing term. Something like that is done in one of the papers I am looking at which explains that the forcing term needs to have a constant added or subtracted off of it if certain conditions are not being met, and that might translate to this. And yes, I did miss that one bit in the code, now it applies.

It’s very convenient that the variational form of the Poisson Equation modified in this way stayed the same, and all that needed to be done was ensuring that the boundary conditions were applied.

Thank you Stein!

[Stein]

@fea

Oh so to make it work, I have to add a correction onto the forcing term

Yes indeed, but that is precisely what this term in your weak form does for you:

@Stein

a01 = ufl.inner(lmbda, du) * ufl.dx

So no further alteration is required.

(Minor comment: I’m noticing a sign mistake in your notation. You’re actually solving
-\nabla^2 u + \lambda = -1 , i.e., with a negative sign in front of the Laplacian. Or, equiv \nabla^2 u = 1 + \lambda )