Integral constraint to replace Neumann Boundary Condition

variational formulation
1

I’d like to replace a Neumann boundary condition with an integral constraint. A simple example is the Laplace equation \Delta c=0 on the unit interval.

I’d like to solve for c with a Dirichlet condition on the right (say c(1)=1.5) and setting the average value of c to 1.0. This should give the solution c=0.5 + x.

I’ve found a couple of examples here and here that I’ve tried to follow. My idea was to add an equation:

\int_x{(c-C_0)\ r\ dx}=0

using an additional test function, r. Based on the examples above, I used a Real finite element for the additional test function. Sample code is given below that does not converge. I believe it is because of the implicit Neumann boundary condition at x=0. Any suggestions would be much appreciated!

from dolfin import *
import matplotlib.pyplot as pl

# Parameters
C0= 1.0  # conserved mean concentration

# Mesh
mesh = IntervalMesh(10,0,1)

# Function space
P2 = FiniteElement('Lagrange', mesh.ufl_cell(), 2)
R = FiniteElement('Real', mesh.ufl_cell(), 0)
element = MixedElement([P2, R])
V = FunctionSpace(mesh, element)

# Functions
u = Function(V)
u.assign(Constant([C0,C0]))
c,m = split(u)
v, r = TestFunctions(V)
f = Constant(0)

# Conservation term
Fc = (c - Constant(C0)) * r * dx

# Laplace Term
Fl = dot(grad(c), grad(v))*dx

# Combined terms
F= Fc+Fl

# Boundary Conditions
def left_boundary(x, on_boundary):
    return on_boundary and near(x[0],0,1e-6)
    
def right_boundary(x, on_boundary):
    return on_boundary and near(x[0],1.0,1e-6)
    
bc0=DirichletBC(V.sub(0), Constant(0.5), left_boundary)
bc1=DirichletBC(V.sub(0), Constant(1.5), right_boundary)
bc=[bc0, bc1]

# Solution with boring DirichletBC
solve(Fl == 0, u, bc)   # this works
plot(c)
pl.show()

# Solution with integral constraint
solve(F == 0, u, bc1)   # this does not converge

plot(c)
pl.show()
2

You can get rid of the implicit Neumann BC at x=0 by including the boundary term from integration by parts in your weak form:

# Laplace Term
n = FacetNormal(mesh)
Fl = dot(grad(c), grad(v))*dx - dot(grad(c),n)*v*ds

That runs and gives the expected solution of c(x) = 0.5 + x.

3

Thanks, for the quick reply. Your solution makes sense but it only works if the problem with Dirichlet BCs is solved first. If I comment out the “boring DirichletBC” section and run the revised code it still does not converge. Any other suggestions?

Here is the revised code based on your suggestion:

from dolfin import *
import matplotlib.pyplot as pl

# Parameters
C0= 1.0  # conserved mean concentration

# Mesh
mesh = IntervalMesh(10,0,1)

# Function space
P2 = FiniteElement('Lagrange', mesh.ufl_cell(), 2)
R = FiniteElement('Real', mesh.ufl_cell(), 0)
element = MixedElement([P2, R])
V = FunctionSpace(mesh, element)

# Functions
u = Function(V)
u.assign(Constant([C0,C0]))
c,m = split(u)
v, r = TestFunctions(V)
f = Constant(0)

# Conservation term
Fc = (c - Constant(C0)) * r * dx

# Laplace Term
Fl0 = dot(grad(c), grad(v))*dx

n = FacetNormal(mesh)
Fl = dot(grad(c), grad(v))*dx - dot(grad(c),n)*v*ds

# Combined terms
F= Fc + Fl

# Boundary Conditions
def left_boundary(x, on_boundary):
    return on_boundary and near(x[0],0,1e-6)
    
def right_boundary(x, on_boundary):
    return on_boundary and near(x[0],1.0,1e-6)
    
bc0=DirichletBC(V.sub(0), Constant(0.5), left_boundary)
bc1=DirichletBC(V.sub(0), Constant(1.5), right_boundary)
bc=[bc0, bc1]

# Solution with boring DirichletBC
# solve(Fl0 == 0, u, bc)   # this works
# 
# plot(c)
# pl.show()

# Solution with integral constraint
solve(F == 0, u, bc1)   # this does not converge

plot(c)
pl.show()
4

Okay, looking a little closer, one thing to notice is that the Lagrange multiplier m for the constraint is absent from the variational form F, so the constraint has no way to influence the PDE problem (weighted against v). Using

# Conservation term
Fc = (c - Constant(C0)) * r * dx + m*v*ds

produces the expected solution in this case, although that looks weird to me, because it doesn’t derive from a Lagrangian. (I would expect something like derivative(m*(c-Constant(C0))*dx,u), which would have + m*v*dx, although that enforces the constraint through a volume term, which leads to a different solution.)

EDIT: Also, if you want m to have the interpretation of a flux, you can drop the boundary term from Fl, although the solution for c stays the same.

5

Yes, thanks that works. I thought I could get away without calculating m, since I could with the Dirichlet BCs.

6

Hello,

I have the same problem. I want to solve the Navier-Stokes equation in a 2D rectangular domain. I would like to have a fully-developed flow at the inlet so I want to impose a constraint regarding the total volumetric flow rate. I followed the steps proposed here but I get errors. The unknowns of the problem are the velocity U(Ux,Uy), pressure P and I also defined a “Real” element for flow rate Q. At the wall I impose no-slip condition and at the inlet/outlet I include the term from the integration by parts. The code is converged for Dirichlet BCs at the inlet. Here is a sample of my code:

*# Define function spaces *
V           = VectorElement('CG'  ,mesh.ufl_cell(),2)
P           = FiniteElement('CG'  ,mesh.ufl_cell(),1)
Q           = FiniteElement('Real',mesh.ufl_cell(),0)
VPQ         = FunctionSpace(mesh,MixedElement([V,P,Q]))
phi,eta,ksi = TestFunctions(VPQ)


       	 
*# Set datum pressure*
def PinPoint(x, on_boundary ):
	tol = 1.e-4
	pointPres = near( x[0], 0.995833    , tol) and near(x[1],  0.005  , tol) 

	if (pointPres is True):
		return (pointPres)


*#BCs*
 Pref = Constant(0.0)
ns2  = Constant((0.0,0.0))


vFluidWall  = DirichletBC(VPQ.sub(0), ns2 ,boundaries, wall  ) 

bcp = DirichletBC(VPQ.sub(1),Pref,PinPoint,'pointwise')

bcs     = [bcp , vFluidWall ]
bcs_    = [bcp , vFluidWall ]



*# Create functions*
W         = Function(VPQ)
[bc.apply(W.vector()) for bc in bcs]
u,p,q     = split(W) 



*#	stresses*
    def tau(X):
	nu = Constant(0.1)
	return nu*(grad(X)+grad(X).T)

def cauchy(X,P):
	return -P*Identity(len(X)) + tau(X)



#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
# 						*weak forms*
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-


dx = Measure("dx")(subdomain_data=subdomains)
ds = Measure("ds")(subdomain_data=boundaries) # outer boundary
dS = Measure("dS")(subdomain_data=boundaries) # inner boundary


n       =   FacetNormal(mesh)
Diff    =   inner( cauchy(u,p) , grad(phi))*dx() 
outFBC  = - inner( dot( cauchy(u,p),n ),phi )*ds(outlet) 
inFBC   = - inner( dot( cauchy(u,p),n ),phi )*ds(inlet)


Cont  = eta*div(u)*dx()
Conv  = inner(grad(u)*u, phi)*dx()


*# Conservation term*
Co = 0.01
Fc = (u - Constant(C0)) * ksi * dx + q*phi*ds(inlet) 


Res    = Diff + Cont + Conv + outFBC + inFBC + Fc
Jac    = derivative(Res,W)
Corr   = Function(VPQ)



*# Solve*
NewtonRaphson(Res,Jac,bcs_,W,Corr,rank)

Thank you in advance