6
$\begingroup$

I am trying to solve a linear equation system $\textbf{A}\textbf{x}=\textbf{b}$, e.g. a Poisson equation discretized in strong form, using biCGstab method.

Since there are only natural Neumann boundary conditions, at least one Dirichlet condition must be given. I simply add condition $\sum_{j} x=0$ to row $i$ in $\textbf{A}$. The result converges, but the error in the gradient field $\nabla x$ at $i$ is very large.

What is the reason for this, and how should I imposition this condition?

Improve this question
$\endgroup$
2
  • 2
    $\begingroup$ Look at Bochev and Lehoucq, On the Finite Element Solution of the Pure Neumann Problem, SIAM Review 2005, epubs.siam.org/doi/abs/10.1137/S0036144503426074. $\endgroup$
    – Christian Clason
    Commented Jun 1, 2015 at 14:05
  • 1
    $\begingroup$ Depending on your preconditioner, you may get away with not specifying such a condition. Krylov solvers will, in some sense, pick a pressure for you. $\endgroup$
    – Bill Barth
    Commented Jun 1, 2015 at 14:31

3 Answers 3

Reset to default
7
$\begingroup$

Pure Neumann problem is unique up to a constant. My two favourite solution strategies:

  1. Modifying the equation $-\Delta u = f$ to $-\Delta u + \varepsilon u = f$ for some small $\varepsilon>0$. If you perform reduced integration this corresponds roughly to adding $\varepsilon$ to the diagonal of your system matrix.

  2. Imposing $\int_\Omega u \,\mathrm{d}x=0$ using Lagrange multiplier. Then your system would be $$\begin{bmatrix} A & 1 \\ 1^T & 0 \end{bmatrix} \begin{bmatrix} x \\ \lambda \end{bmatrix} = \begin{bmatrix} f \\ 0 \end{bmatrix}.$$

Improve this answer
$\endgroup$
4
  • $\begingroup$ Just a note: BiCGstab won't care, only CG requires symmetry and positive definiteness. (On the other hand, CG will work for the unmodified system, as long as the (discrete) system is consistent -- at least in exact arithmetic, and up to a constant.) $\endgroup$
    – Christian Clason
    Commented Jan 28, 2016 at 19:33
  • $\begingroup$ Thanks for the clarification. I made a faulty assumption that this property is shared by CG and BiCGstab although I don't know the latter method well. $\endgroup$
    – knl
    Commented Jan 28, 2016 at 20:12
  • $\begingroup$ @knl I tried to follow your approach in my last question, however it's not apparently working $\endgroup$
    – FEGirl
    Commented Oct 3, 2021 at 10:52
  • $\begingroup$ The vector $1^{\text{T}}$ is only valid if each DOF has the same weight. A more general solution would be to consider the correct quadrature weights. Moreover, if the domain is not Cartesian one would also have to consider the correct Jacobian weights. Same has to be considered for method 1 using a diagonal "jitter". $\endgroup$
    – ConvexHull
    Commented Jul 27, 2022 at 13:16
4
$\begingroup$

Including even one Dirichelt condition changes the problem you are trying to solve and will not give you the correct solution! You must fulfil the Discrete Compatibility Criteria, see e.g. first and second pages of: http://eprints.ma.man.ac.uk/894/2/covered/MIMS_ep2007_156_Sample_Chapter.pdf

It basically states that the summation of each element of b must be equal to zero for all-Neumann problems, otherwise your solution will drift as there is nothing holding it at the boundaries.

You simply need to add a line stating that b = b - mean(b) in your code before solving, if you are solving in double precision. In single precision it might also be necessary to ensure, at every other iteration, that the residual in biCGstab also meets the Discrete Compatibility Criteria.

I attempted to solve an all-Neumann problem using CG by imposing a Dirichlet condition at a single point, as one is usually told to do anecdotally by e.g. some Prof., and then compared it to using a Discrete Cosine Transform approach. It does not yield the same result, using the Discrete Compatibility Criteria instead does.

Improve this answer
$\endgroup$
1
  • $\begingroup$ Could you clarify what you mean by: "In single precision it might also be necessary to ensure, at every other iteration, that the residual in biCGstab also meets the Discrete Compatibility Criteria."? $\endgroup$
    – lightxbulb
    Commented May 16, 2021 at 15:06
0
$\begingroup$

I remember that the matrices from the Poisson equations with pure Neumann boundary conditions is singular. Is it right for your case ? If so, you cannot use the BiCGSTAB to solve the resulting linear systems.

Improve this answer
$\endgroup$
7
  • $\begingroup$ Yes, it is singular, and that is why I have to add one Dirichlet bc. Because I only care about the gradient field. Maybe I can use QR solvers, but they are too slow. $\endgroup$
    – Kozuki
    Commented Jun 2, 2015 at 4:59
  • $\begingroup$ Well, if the matrix is singular, maybe I suggest you use some regularization methods, such as Tikhonov regularization. $\endgroup$
    – Hsien-Ming Ku
    Commented Jun 2, 2015 at 7:20
  • $\begingroup$ The resulting matrix will still be singular, even with the additional condition added. $\endgroup$
    – shuhalo
    Commented Jun 3, 2015 at 7:44
  • 1
    $\begingroup$ I case someone still stumbles upon this answer, it is incorrect and Conjugate Gradients works for pure Neumann problems, see e.g. cs.sandia.gov/~pbboche/papers_pdf/2005SIREV.pdf $\endgroup$
    – DrHansGruber
    Commented Jan 21, 2017 at 23:54
  • 1
    $\begingroup$ @DrHansGruber The link is dead. Could you tell me the title of the paper? Nvm, used the wayback machine, it's: ON THE FINITE ELEMENT SOLUTION OF THE PURE NEUMANN PROBLEM by Pavel Bochev. $\endgroup$
    – lightxbulb
    Commented May 15, 2021 at 23:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.