I’m trying to solve a PDE Lf=0 (for some operator L) with the mass constraint \int f dx=1. I suspect there is a non-trivial solution f\not\equiv 0 but I know for sure that f\equiv 0 is also a solution for my model. I wrote the below implementation with a Lagrange multiplier, but when I test the residual against f it shows that the PDE residual (excluding the Lagrange multiplier part) is of the same order (10e-1) as that of the Lagrange multiplier, whereas the total residual is 10e-14 after using solve. How can I ensure that both residuals are small, i.e. the PDE residual and the total residual are both small?
from dolfin import *
import matplotlib.pyplot as plt
# Sub domain for Dirichlet boundary condition
class NofluxBoundary(SubDomain):
def inside(self, x, on_boundary):
return bool((x[0] < DOLFIN_EPS or x[0] > (1.0 - DOLFIN_EPS)) \
and on_boundary)
# Sub domain for Periodic boundary conditional
class PeriodicBoundaryC(SubDomain):
# Left boundary is "target domain" G
def inside(self, x, on_boundary):
return bool(x[1] < DOLFIN_EPS and x[1] > -DOLFIN_EPS and on_boundary)
# Map right boundary (H) to left boundary (G)
def map(self, x, y):
y[1] = x[1] - 1.0
y[0] = x[0]
class PeriodicBoundary(SubDomain):
# Left boundary is "target domain" G
def inside(self, x, on_boundary):
return bool(x[1] < DOLFIN_EPS and x[1] > -DOLFIN_EPS and on_boundary)
# Create mesh and finite element
mesh = UnitSquareMesh(32, 32)
pbc = PeriodicBoundaryC()
nofluxb = NofluxBoundary()
pb = PeriodicBoundary()
Q = FiniteElement("CG", triangle, 1)
LM = FiniteElement("Real", triangle, 0)
M = FunctionSpace(mesh, MixedElement([Q, LM]), constrained_domain=pbc)
# Create domain and boundaries
domains = MeshFunction("size_t", mesh, mesh.topology().dim())
domains.set_all(0)
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
boundaries.set_all(0)
pb.mark(boundaries, 1)
nofluxb.mark(boundaries, 2)
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)
drc = Expression('-sin(pi*x[0])',degree=2)
# Create form
u = Function(M)
f, lm = split(u)
wt = TestFunction(M)
v, lmt = split(wt)
r = Expression('x[0]',degree=2)
cospsi = Expression('cos(2*pi*x[1])',degree=2)
sinpsi = Expression('sin(2*pi*x[1])',degree=2)
B = r*cospsi*grad(f)[0]*v*dx - sinpsi*grad(f)[1]*v*dx \
+ r*sinpsi*drc*f*grad(v)[1]*dx \
- r*grad(f)[0]*v*ds(1) \
+ lm*v*dx + lmt*f*dx - lmt*1*dx
solve(B == 0, u)
f,lm = u.split()
# Check residuals testing against f itself
Bf_tot = r*cospsi*grad(f)[0]*f*dx - sinpsi*grad(f)[1]*f*dx \
+ r*sinpsi*drc*f*grad(f)[1]*dx \
- r*grad(f)[0]*f*ds(1) \
+ lm*f*dx + lm*f*dx - lm*1*dx
Bf = r*cospsi*grad(f)[0]*f*dx - sinpsi*grad(f)[1]*f*dx \
+ r*sinpsi*drc*f*grad(f)[1]*dx \
- r*grad(f)[0]*f*ds(1)
print("f tot res = {:.2e}".format(abs(assemble(Bf_tot))))
print("f PDE res = {:.2e}".format(abs(assemble(Bf))))
# lm is equal to a real value everywhere on the domain, here we take (0.5,0.5)
print("lm = {:.2e}".format(lm(0.5,0.5)))