Confused about singular poisson example

Question

Hello,

I'm struggling to understand the singular-poisson example.

As I understand, specifying only Neumann boundary conditions means that the solution is only determined up to a constant term. In other words, I can find two solutions u1 and u2 so that a(u1, v) = a(u2, v) = L(v) for all v.

I see that this implies that A has a non-zero null space. However, I do not understand why this causes a problem that needs to be handled in practice.

It seems to me that if b is calculated by evaluating L(v) for all v, it automatically lies in the range of A, and there is no need to subtract anything. As a matter of fact, when I coded some simple examples ignoring the singularity of A, FEniCS calculated the proper solution just fine. If we can find a b not in the range of A, wouldn't this mean that there is no solution to the original PDE with the given boundary conditions?

In addition to that, the example code...

# Create vector that spans the null space
null_vec = Vector(u.vector())
V.dofmap().set(null_vec, 1.0)
null_vec *= 1.0/null_vec.norm("l2")
[...]
# Orthogonalize b with respect to the null space (this gurantees that
# a solution exists)
null_space.orthogonalize(b);

...seems to use the null space of A to orthogonalize b. Shouldn't this be the nullspace of transpose(A)? (cf. issue 186)

Finally, I also do not understand how ensuring that b is in the range of A fixes the constant term of the solution. Even with the new b, there will still be two solutions u1 and u2 such that A.u1 = A.u2 = b.

Can someone explain what I'm missing here?

Answer 1

With some help from the mailing list, I was able sort out most of my questions:

• As long as the solution is in A's range, A having a nullspace is not a problem per se. However, some solvers may struggle if they do not know about the nullspace.
• If the solution vector b is not in the range of A, then the original PDE indeed does not have a solution with the given boundary conditions. However, boundary condition that satisfy the compatibility condition analytically may no longer satisfy it after discretization. In this case it is desirable to project out the error.
• The code to orthogonalize b is correct because A is symmetric, so the nullspace of A and transpose(A) is the same.
• This example does not fix the constant term in the solution.

Answer 2

You can fix one point's (on boundary) value as 0

bc = DirichletBC(Th, value, Position, "pointwise")

then you can solve you PDE without errors.

Another way is you can change your PDE with
$$-\Delta u +0.0000001u=f$$
after this you will get that A is full rank
You will find $\int_{\Omega} u d \Omega =0$

When you solve Stokes equation you will find the pressure need this trick

Comments

Thanks! However, I am actually already able to solve the PDE without errors (as far as I can tell). What I am trying to understand is what is being done (and why it is being done) in the singular-poisson example.