Poisson solution reaches minimum residual after 1 iteration, then increases

1

I am a complete newbie to FEM and FEniCS. I’m trying to solve a Poisson equation with Neumann BC. The variational statement of my problem is equation (A.1) here with I_n = 0, \rho_e = 0, and I_m = 0.1

The boundary \Gamma is a circle within a box \partial\Omega_N, and the space between the circle and the box is the domain \Omega_e (see plot of solution below). Since my BC is zero everywhere except the circle, I have implemented it by writing a UserExpression boundary which returns 0 everywhere except the circle, and I integrate this over the entire boundary \Gamma \cup \partial\Omega_N

Here is my FEniCS implementation:

from fenics import *
from mshr import *

# Generate geometry, mesh
CIR_X, CIR_Y, CIR_R = 0.5, 0.3, 0.1
box = Rectangle(Point(0, 0), Point(1, 1))
circle = Circle(Point(CIR_X, CIR_Y), CIR_R)
domain = box - circle 
mesh = generate_mesh(domain, 256)

# Define the boundary term
class BoundaryVals(UserExpression):
    def eval(self, value, x):
        rad_sq = (x[0] - CIR_X)**2 + (x[1] - CIR_Y)**2
        if near(rad_sq, CIR_R**2):
            value[0] = 0.1
        else:
            value[0] = 0
boundary = BoundaryVals(degree=2)

# Define the problem
Omega = FunctionSpace(mesh, 'CG', 2)
sigma = 1
Theta = Function(Omega)
v = TestFunction(Omega)
LHS = sigma * inner(grad(Theta), grad(v)) * dx - boundary * v * ds

solve(LHS == 0, Theta)

Here is the output before it reaches the max # of iterations:

WARNING: user expression has not supplied value_shape method or an element. Assuming scalar element.
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) = 8.132e-04 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
  Newton iteration 1: r (abs) = 4.875e-05 (tol = 1.000e-10) r (rel) = 5.994e-02 (tol = 1.000e-09)
  Newton iteration 2: r (abs) = 4.921e-05 (tol = 1.000e-10) r (rel) = 6.051e-02 (tol = 1.000e-09)
  Newton iteration 3: r (abs) = 5.313e-05 (tol = 1.000e-10) r (rel) = 6.534e-02 (tol = 1.000e-09)
  Newton iteration 4: r (abs) = 4.697e-05 (tol = 1.000e-10) r (rel) = 5.776e-02 (tol = 1.000e-09)
  Newton iteration 5: r (abs) = 5.392e-05 (tol = 1.000e-10) r (rel) = 6.631e-02 (tol = 1.000e-09)
  Newton iteration 6: r (abs) = 5.692e-05 (tol = 1.000e-10) r (rel) = 6.999e-02 (tol = 1.000e-09)
  Newton iteration 7: r (abs) = 5.317e-05 (tol = 1.000e-10) r (rel) = 6.538e-02 (tol = 1.000e-09)
  Newton iteration 8: r (abs) = 4.929e-05 (tol = 1.000e-10) r (rel) = 6.061e-02 (tol = 1.000e-09)
  Newton iteration 9: r (abs) = 6.195e-05 (tol = 1.000e-10) r (rel) = 7.618e-02 (tol = 1.000e-09)
  Newton iteration 10: r (abs) = 5.791e-05 (tol = 1.000e-10) r (rel) = 7.121e-02 (tol = 1.000e-09)
  Newton iteration 11: r (abs) = 5.737e-05 (tol = 1.000e-10) r (rel) = 7.054e-02 (tol = 1.000e-09)
  Newton iteration 12: r (abs) = 7.151e-05 (tol = 1.000e-10) r (rel) = 8.794e-02 (tol = 1.000e-09)
  Newton iteration 13: r (abs) = 6.114e-05 (tol = 1.000e-10) r (rel) = 7.519e-02 (tol = 1.000e-09)
  Newton iteration 14: r (abs) = 6.960e-05 (tol = 1.000e-10) r (rel) = 8.559e-02 (tol = 1.000e-09)
  Newton iteration 15: r (abs) = 6.560e-05 (tol = 1.000e-10) r (rel) = 8.067e-02 (tol = 1.000e-09)
  Newton iteration 16: r (abs) = 6.360e-05 (tol = 1.000e-10) r (rel) = 7.821e-02 (tol = 1.000e-09)
  Newton iteration 17: r (abs) = 7.650e-05 (tol = 1.000e-10) r (rel) = 9.407e-02 (tol = 1.000e-09)
  Newton iteration 18: r (abs) = 7.797e-05 (tol = 1.000e-10) r (rel) = 9.588e-02 (tol = 1.000e-09)
  Newton iteration 19: r (abs) = 8.459e-05 (tol = 1.000e-10) r (rel) = 1.040e-01 (tol = 1.000e-09)
  Newton iteration 20: r (abs) = 7.426e-05 (tol = 1.000e-10) r (rel) = 9.131e-02 (tol = 1.000e-09)
  Newton iteration 21: r (abs) = 8.411e-05 (tol = 1.000e-10) r (rel) = 1.034e-01 (tol = 1.000e-09)
  Newton iteration 22: r (abs) = 7.860e-05 (tol = 1.000e-10) r (rel) = 9.666e-02 (tol = 1.000e-09)
  Newton iteration 23: r (abs) = 8.671e-05 (tol = 1.000e-10) r (rel) = 1.066e-01 (tol = 1.000e-09)
  Newton iteration 24: r (abs) = 1.008e-04 (tol = 1.000e-10) r (rel) = 1.240e-01 (tol = 1.000e-09)
  Newton iteration 25: r (abs) = 9.876e-05 (tol = 1.000e-10) r (rel) = 1.214e-01 (tol = 1.000e-09)
  Newton iteration 26: r (abs) = 9.744e-05 (tol = 1.000e-10) r (rel) = 1.198e-01 (tol = 1.000e-09)
  Newton iteration 27: r (abs) = 1.002e-04 (tol = 1.000e-10) r (rel) = 1.232e-01 (tol = 1.000e-09)
  Newton iteration 28: r (abs) = 1.088e-04 (tol = 1.000e-10) r (rel) = 1.337e-01 (tol = 1.000e-09)
  Newton iteration 29: r (abs) = 1.028e-04 (tol = 1.000e-10) r (rel) = 1.265e-01 (tol = 1.000e-09)
  Newton iteration 30: r (abs) = 1.083e-04 (tol = 1.000e-10) r (rel) = 1.331e-01 (tol = 1.000e-09)
  Newton iteration 31: r (abs) = 1.092e-04 (tol = 1.000e-10) r (rel) = 1.343e-01 (tol = 1.000e-09)
  Newton iteration 32: r (abs) = 1.028e-04 (tol = 1.000e-10) r (rel) = 1.264e-01 (tol = 1.000e-09)
  Newton iteration 33: r (abs) = 1.210e-04 (tol = 1.000e-10) r (rel) = 1.488e-01 (tol = 1.000e-09)
  Newton iteration 34: r (abs) = 1.055e-04 (tol = 1.000e-10) r (rel) = 1.297e-01 (tol = 1.000e-09)
  Newton iteration 35: r (abs) = 1.022e-04 (tol = 1.000e-10) r (rel) = 1.257e-01 (tol = 1.000e-09)
  Newton iteration 36: r (abs) = 1.334e-04 (tol = 1.000e-10) r (rel) = 1.641e-01 (tol = 1.000e-09)
  Newton iteration 37: r (abs) = 1.436e-04 (tol = 1.000e-10) r (rel) = 1.766e-01 (tol = 1.000e-09)
  Newton iteration 38: r (abs) = 1.275e-04 (tol = 1.000e-10) r (rel) = 1.567e-01 (tol = 1.000e-09)
  Newton iteration 39: r (abs) = 1.311e-04 (tol = 1.000e-10) r (rel) = 1.612e-01 (tol = 1.000e-09)
  Newton iteration 40: r (abs) = 1.299e-04 (tol = 1.000e-10) r (rel) = 1.598e-01 (tol = 1.000e-09)
  Newton iteration 41: r (abs) = 1.444e-04 (tol = 1.000e-10) r (rel) = 1.776e-01 (tol = 1.000e-09)
  Newton iteration 42: r (abs) = 1.346e-04 (tol = 1.000e-10) r (rel) = 1.655e-01 (tol = 1.000e-09)
  Newton iteration 43: r (abs) = 1.356e-04 (tol = 1.000e-10) r (rel) = 1.667e-01 (tol = 1.000e-09)
  Newton iteration 44: r (abs) = 1.306e-04 (tol = 1.000e-10) r (rel) = 1.606e-01 (tol = 1.000e-09)
  Newton iteration 45: r (abs) = 1.503e-04 (tol = 1.000e-10) r (rel) = 1.849e-01 (tol = 1.000e-09)
  Newton iteration 46: r (abs) = 1.306e-04 (tol = 1.000e-10) r (rel) = 1.606e-01 (tol = 1.000e-09)
  Newton iteration 47: r (abs) = 1.780e-04 (tol = 1.000e-10) r (rel) = 2.189e-01 (tol = 1.000e-09)
  Newton iteration 48: r (abs) = 1.980e-04 (tol = 1.000e-10) r (rel) = 2.435e-01 (tol = 1.000e-09)
  Newton iteration 49: r (abs) = 2.151e-04 (tol = 1.000e-10) r (rel) = 2.645e-01 (tol = 1.000e-09)
  Newton iteration 50: r (abs) = 2.053e-04 (tol = 1.000e-10) r (rel) = 2.524e-01 (tol = 1.000e-09)

If I stop after the first iteration by manually setting the tolerance like this:

solve(LHS == 0, Theta, solver_parameters={"newton_solver": {"absolute_tolerance": 5.1e-5}})

The solution after 1 iteration looks reasonable:

862×826

I’m surprised this solution was found after just 1 iteration. Should I take that as a sign that the solution is probably not accurate? Or is it just a consequence of my simple geometry?

Other questions:

  • How can I compute the normal derivatives along the boundaries in order to verify the solution? (I know I can “call” the solution to approximate the derivative like this: [Theta(CIR_X+CIR_R+h, CIR_Y) - Theta(CIR_X+CIR_R, CIR_Y)] / h but is there another way?)
  • What is the meaning of “relative tolerance”? Relative to what?
  • Why does the residual increase for all most iterations after the first? (I have assumed that the r it reports after each iteration means “residual”)
  • How is r computed?
  • When I start solving the problem on more complicated geometries, how should I think about choosing a tolerance?

Any insight you can provide will likely be helpful because this is all new to me. Thanks in advance.

2

For a Poisson problem, the solution is not unique if you only have Neumann conditions, see for instance:
Poisson with pure Neumann conditions.

3

Thanks. I added the constraint that \int \Theta = 0 as described in that tutorial, and I am able to compute a solution without the iterative solver (solve() doesn’t report the results of an iterative solver, so I assume it’s not being used). I have to admit I don’t quite understand exactly what’s going on (what did the iterative solver do? Was it an attempt at finding the best solution given the non-uniqueness?) but I got a similar solution and I think it’s working now