Iterative Newton solver for Poisson eqn fails with Neumann boundary

Question

The nonlinear-poisson demo gives an example of using an iterative Newton solver to solve a nonlinear Poisson equation with Dirichlet boundary conditions. I want to use the the solver with Neumann conditions, since the iterative solver solves for Function u rather than TrialFunction, which seems to make it easier to setup a nonlinear source function.

I can get a Neumann solution using direct solvers (via TrialFunction). But when I try Neumann conditions with a Newton solver, it immediately diverges:

No Jacobian form specified for nonlinear variational problem.
Differentiating residual form F to obtain Jacobian J = F'.
Solving nonlinear variational problem.
  Newton iteration 0: r (abs) = 9.896e-01 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-06)
  Newton iteration 1: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-06)
  Newton iteration 2: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-06)

I've tried variations such as solving the Dirichlet problem first, then using that solution as an initial guess for the Neumann problem. But that gives the same kind of divergence, even if I set the Neumann condition using the gradient from the Dirichlet solution.

The simplest test case, reduced to 1D, is

from dolfin import *

# Sub domain for Dirichlet boundary condition
class BoundaryRight(SubDomain):
    def inside(self, x, on_boundary):
        return abs(x[0] - 1.0) < DOLFIN_EPS and on_boundary

# Create mesh and define function space
mesh = UnitIntervalMesh(32)

V = FunctionSpace(mesh, "CG", 1)

boundaries = FacetFunction("size_t", mesh)
boundaries.set_all(0)
boundaryRight = BoundaryRight()
boundaryRight.mark(boundaries, 1)
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)

g = Constant(1.0)

# Define variational problem
u = Function(V)
v = TestFunction(V)
f = Expression("x[0]", degree=2)
F = inner(grad(u), grad(v))*dx - f*v*dx + g*v*ds(1)

# Compute solution
solve(F == 0, u, [], solver_parameters={"newton_solver": {"relative_tolerance": 1e-6}})

plot(u, title="Solution Neumann")
plot(grad(u), title="Solution Neumann gradient")
interactive()

Would you expect this code to converge?

Comments

WIthout a dirichlet boundary condition, your problem is not well-posed.

---

In addition to nate's comment I would add that the post here, whilst relating to the linear Poisson problem with pure Neumann BC's, has some valuable information and discusses different approaches which may help in your case.

---

Thanks guys, that link makes sense. The Neumann demo discussess it too at https://fenicsproject.org/documentation/dolfin/2016.1.0/python/demo/documented/neumann-poisson/python/documentation.html

I had been thinking the condition ∫Ωfdx=−∫∂Ωgds was met automatically, but apparently it needs to be imposed explicitly. The strange thing is that I thought I did have a Neumann boundary working with a direct solver (NonlinearVariationalProblem with NonlinearVariationalSolver). I'll have to go back and check what I was doing in that case.

---

I've identified why I'm confused. My code sample above does converge successfully if I use a nonlinear source term

f = exp(-u)

in place of

f = Expression("x[0]", degree=2)

Does the nonlinear source term necessarily make the problem well-posed, or is it a phantom solution?

Answer 1

I think the issue is your original problem was not well posed. For a solution to exist one must have
$$ -g(1)+g(0) = \int_{0}^{1} f dx $$
which is obtained by integrating your original problem
$$
-\frac{\partial^{2} u}{\partial x^{2}} = f
$$
with respect to x and applying the boundary conditions.
In your original problem with $f=x$, $g(0)=0$ and $g(1)=1$, this condition is not satisfied and so no solutions exist. Changing $g(1)$ to $-1/2$ would ensure solutions exist but there are now infinitely many of them and any solver used would need to be able to deal with this, e.g. a Krylov solver would work as it is a linear problem. A quick test indicates that Newton's method does not work here.

Determining whether solutions exist and are unique for an $f$ depending on $u$ is more difficult. It would appear a unique solution exists in the case $f=e^{-u}$ to which Newton's method converges.

Comments

Thanks Brenden. That's in interesting point about the solutions for g(1) = −1/2. I should have a closer look at the Krylov solvers, I haven't worked with them yet.