If the constraint does not involve a derivative, like say \int_\Omega u dx = C, the constraint is satisfied to within machine precision, however this is not the case if a derivative is involved. I could always just increase the order of the functionspace, but that is highly expensive, and depending on the use case, that can make it highly impractical. So, I want to know if there are any ways to get around this issue.
Now, let’s get into the PDE and provide the code. Our Poisson Equation with its integral constraint is
\nabla^2 u = -1
u = 0 \ \ \ \ \forall (x,y)\in \partial \Omega
\oint \nabla u\cdot\hat{\textbf{n}}ds = h = 1.05
(As user @spcvanschie pointed out in my previous post, this is a Poisson equation with Dirichlet boundary conditions, ergo its solution is unique, so to add an integral constraint without overdetermining the solution, the forcing term will need to be modified by our Lagrange multiplier \nabla^2 u = -1 +\lambda.) (This specific Poisson Equation I did not come up with using any physical problem as a reference, it is one I came up with using arbitrarily values for everything.) We can now set up the variational form of our constrained Poisson Equation using the Lagrange Multiplier method:
\int_\Omega \nabla u \cdot \nabla v dx + \int_\Omega \lambda v dx = \int_\Omega fv dx
\int_\Omega 0\cdot u \delta\lambda dx + \int_\Omega \lambda \delta\lambda dx = \int_\Omega f \delta\lambda dx + h \delta\lambda.
Setting up the matrices
A = \pmatrix{\int_\Omega \nabla u \cdot \nabla v dx & \int_\Omega \lambda v dx \\
\int_\Omega 0\cdot u \delta\lambda dx & \int_\Omega \lambda \delta\lambda dx}
L = \pmatrix{\int_\Omega fv dx \\ \int_\Omega f \delta\lambda dx}
and we build the constrained mixed variational form and compare the constrained solution to the original in (the latest version of) Dolfinx (where I have primarily taken code from this tutorial to help me along)
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
from dolfinx.fem.petsc import LinearProblem
import pyvista4dolfinx as p4d
Nx,Ny = 20,20
mesh = df.mesh.create_unit_square(
MPI.COMM_WORLD, Nx, Ny, df.mesh.CellType.triangle, dtype=np.float64
)
n = ufl.FacetNormal(mesh) # unit normal vector (used later)
t = ufl.perp(n)
order = 2
# Real space
V = df.fem.functionspace(mesh, ("CG", order))
R = create_real_functionspace(mesh)
# Mixed space
W = ufl.MixedFunctionSpace(V, R)
u, lmbda = ufl.TrialFunctions(W)
du, dl = ufl.TestFunctions(W)
#---- 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 = 1.05
# Forcing term
f = df.fem.Constant(mesh, df.default_scalar_type(-1))
#--------------------------------------------------
#---- 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 = df.fem.Constant(mesh, df.default_scalar_type(0)) * u * dl * ufl.dx
a11 = ufl.inner(lmbda,dl) * ufl.dx
# ---- Build components: linear form ----
L0 = ufl.inner(f, du) * ufl.dx
L1 = ufl.inner(f,dl) * ufl.dx
# ---- Build matrix: A ----
bcs = [bc]
a = a00 + a01 + a10 + a11
a_blocks = ufl.extract_blocks(a)
a_compiled = df.fem.form(a_blocks)
A = df.fem.petsc.assemble_matrix(a_compiled, bcs=bcs)
A.assemble()
# ---- Build matrix: L ----
L = L0 + L1
L_blocks = ufl.extract_blocks(L)
L_compiled = df.fem.form(L_blocks)
#--------------------------------------------------
#---- 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()
#---- Solve ----
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()
#--------------------------------------------------
#---- Compare -------------------------------------
#--------------------------------------------------
# ---- Uncontstrained problem ----
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
a = ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx
L = f * v * ufl.dx
problem2 = LinearProblem(
a,
L,
bcs=[bc],
petsc_options={"ksp_type": "preonly", "pc_type": "lu"},
petsc_options_prefix="Poisson",
)
psi_I = problem2.solve()
# ---- Compare fluxes ----
print("-"*40,"\nNo constraint")
print("-"*40)
print("∫∇u ⋅ n ds = ",df.fem.assemble_scalar(df.fem.form(ufl.dot(n,ufl.grad(psi_I)) * ufl.ds)),"\n")
print("-"*40,f"\nWith constraint C(u):= ∫∇u⋅n ds = {h}")
print("-"*40)
print(f"∫∇u_c ⋅ n ds (should be approximately {h}) = ",df.fem.assemble_scalar(df.fem.form(ufl.dot(n,ufl.grad(uh)) * ufl.ds )))
Additionally, the constrained potential preserved desired properties which the unconstrained one satisfied. I also produced plots seen below which show that their gradients are nearly identical (this is a property I also wanted to see in my constrained solution, so this is a good thing)
The only issue is that for the enforced condition, if we take a look at what is printed when we run the code block, it does not satisfy it to within machine precision like want. The result is 1.0471563447846715 instead of 1.05. I want to know if there are any known ways to get around this unwanted error that arises when a derivative is involved with the constraint. It is highly important that integral constraints are as closely satisfied as possible for what I am doing.
