DOLFIN

Laplace equation with von Neumann boundary condition

Asked by Peter Franek

I tried to solve an equation Laplace(u)=0 with a boundary condition on derivative, \partial_n u=u_0. This equation has a solution iff \int u_0 ds=0. However, the solution is determined only up to a constant. How to do it in FEniCS? I tried to formulate it as a variational problem with no Dirichlet boundary condition, but it computed a result that contained very large numbers (like E13). Is there a way to "tell him" to look for the solution not in the space ("Lagrange", 1) but in a space of functions with zero integral, or so?

Question information

Language:
English Edit question
Status:
Solved
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Solved by:
Marie Rognes
Solved:
Last query:
Last reply:
Revision history for this message
Best Marie Rognes (meg-simula) said : 		#1

Yes, see

  demo/undocumented/neumann-poisson/python/demo_neumann-poisson.py

(or the corresponding cpp version) for an example

Revision history for this message
Jan Blechta (blechta) said : 		#2

Or consider less costly option of imposing Dirichlet condition in "one point", for example
    DirichletBC(V, 0.0, "x[0]< DOLFIN_EPS && x[1] < DOLFIN_EPS", "pointwise")

Revision history for this message
Peter Franek (peter-franek) said : 		#3

Thanks Marie Rognes, that solved my question.

Revision history for this message
Nico Schlömer (nschloe) said : 		#4

I think the example should be revised a little bit:
If the system is consistent then any Krylov method should converge. The problem space is effectively reduced by one dimension.

An advantage of not adding an extra condition is that the structure of the equation system isn't touched. AMG preconditioners can perform well as ever. One has to make sure though that the preconditioner preserves the null space, so you can't use ILU, not even as a coarse-grid solver in AMG.

Here's the mod on the demo problem that does away with the augmentation:

================ *snip* ================
from dolfin import *

# Create mesh and define function space
mesh = UnitSquareMesh(64, 64)
V = FunctionSpace(mesh, "CG", 1)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)")
g = Expression("-sin(5*x[0])")
a = (inner(grad(u), grad(v)))*dx
L = f*v*dx + g*v*ds

# Compute solution
w = Function(V)
A = assemble(a)
b = assemble(L)
null = Function(V)
null.interpolate(Constant(1.0))
nullvec = null.vector()
nullvec /= norm(nullvec)
alpha = b.inner(nullvec)
normB = norm(b)
if alpha/normB > 1.0e-10:
    warning('System not consistent! <b,n>/||b|| = %g. Projecting out inconsistent component.'
            % (alpha/normB)
           )
    b -= alpha * nullvec

prec = PETScPreconditioner('hypre_amg')
prec.parameters['hypre']['BoomerAMG']['relax_type_coarse'] = 'jacobi'
solver = PETScKrylovSolver('cg', prec)
solver.parameters['absolute_tolerance'] = 0.0
solver.parameters['relative_tolerance'] = 1.0e-10
solver.parameters['maximum_iterations'] = 100
solver.parameters['monitor_convergence'] = True
# Create solver and solve system
A_petsc = as_backend_type(A)
b_petsc = as_backend_type(b)
w_petsc = as_backend_type(w.vector())
solver.set_operator(A_petsc)
solver.solve(w_petsc, b_petsc)

# Plot solution
plot(w, interactive=True)
================ *snap* ================