Timeline for Convergence stall when solving 2D Poisson PDE with pure Neumann boundaries (finite differences)
Current License: CC BY-SA 4.0
15 events
| when toggle format | what | by | license | comment | |
|---|---|---|---|---|---|
| Sep 24, 2022 at 15:22 | comment | added | Akhaim | @ConvexHull, OK, I managed to fix NaN issue, as I didn't set the C parameter everywhere in the same way when I added epsilon. However, with this correction convergence stall was back on despite epsilon presence. Is there a closed-form way to estimate epsilon order of magnitude for given discretization parameters (unit cell sizes dx and dy), or do we just guess it as a small positive value? Currently, C=2(1/dx/dx+1/dy/dy+epsilon). If epsilon is 1e-6 (I went as high as 1e-3), and dx=dy=0.1, they completely overshadow epsilon, as C becomes 400.001 (epsilon = 1e-3) or 400.000001 (epsilon 1e-6) | |
| Sep 24, 2022 at 15:09 | comment | added | ConvexHull | Generally this should work for pure Neumann BCs. It is simply a workaround to make your matrix invertible. Try to vary the size. If this does not work you may have another issue with your code. | |
| Sep 24, 2022 at 14:31 | comment | added | Akhaim | @ConvexHull, you are correct about SOR being iterative matrix inversion method, and I have possibly expressed myself wrong. My thoughts were along the lines that perhaps due to the way how SOR is implemented that your first proposed solution (from your original answer) needs to be added in some-kind of a "SOR-friendly" way. Unfortunately, I still do not quite understand what would be the proper way of implementing such solution (Tikhonov regularization) in my code. I've tried directly, just to add some small epsilon (1e-6) in the term C (in my notation), but the solution just diverged to NaN | |
| Sep 24, 2022 at 13:58 | comment | added | ConvexHull | Yes, the proposed methods are quite general and applicable to other kinds of PDE's. Don't get me wrong, however, I think you have a wrong understanding what SOR does. The SOR algorithm is nothing more than a iterative matrix inversion. To be more correct, the SOR solves your system in an iterative sense in order to get $x$. | |
| Sep 24, 2022 at 13:45 | comment | added | Akhaim |
@ConvexHull, the proposed ideas from your answer, are they general for any kind of solver that deal with problems of type Ax=b, or are the specific solver families the only ones benefiting from these approaches? My impression is that adding small positive constant to the diagonal will resolve the issue when one needs to find an inverse of a matrix. In SOR, method I use (which is mentioned in my question, if you read it) there is no matrix inversion, at least explicit. In fact, after trying this modification, I obtain all NaNs as a solution, which means it is not in fact a generic fix.
|
|
| Sep 23, 2022 at 18:41 | comment | added | ConvexHull | See my related answer: mathematica.stackexchange.com/a/271682/87543 | |
| Sep 23, 2022 at 12:44 | comment | added | Akhaim | @helloworld922 for the Dirichlet's BC (DBC), what I did was to select a single corner point and force DBC to be zero there, while preserving all Neumann boundaries, as they were. The resulting solution was then dominated by the DBC, and it's profile was overall incorrect. If neccessary, I could document these findings in a separate answer to this question | |
| Sep 23, 2022 at 12:41 | comment | added | helloworld922 | I mean you need to pick one equation (say the equation for the BC for the ghost node -1,0), and replace it with a row of all 1's and the RHS is 0. This is equivalent to applying a BC $\int \phi = 0$ | |
| Sep 23, 2022 at 12:33 | comment | added | Akhaim | @helloworld922 do you mean, for instance, to replace boundary at e.g. entire left side with a value that is equal to a mean of the remaining boundaries (and possibly with inverse sign)? If so, won't that mess with the expected boundaries? My implementation was originally to modify the RHS prior to the input to the solver as RHS-mean(RHS), but that did messed up with the boundaries. | |
| Sep 23, 2022 at 12:30 | comment | added | helloworld922 | A similar thing applies to fixing a single node to a Dirichlet value; you need to modify the A matrix by replacing one of the Neumann BCs with a Dirichlet BC. | |
| Sep 23, 2022 at 12:27 | comment | added | helloworld922 | For the second solution did you remove one of the BC equations in your system and replace it with the integration condition? I don't think you can just perform an iteration of SOR on the original system then remove the mean value. | |
| Sep 23, 2022 at 12:20 | comment | added | Akhaim | @helloworld922, the second option you proposed, I have already mentioned in my question to not yield any results. The first option, though not mentioned explicitly in my question, is in the linked Stack Exchange post. That option was tested and it did not work, too, as I have also stated in my original post. It is for this reason I have went to great lengths to provide so many implementation details from my code and my discretization scheme. So please, if you can provide a bit more insightful comment, I would appreciate it greatly. | |
| Sep 23, 2022 at 11:07 | comment | added | helloworld922 | The general way to deal with the non-convergence is to add some extra constraint. For example, you could fix a single node in a corner to a Dirichlet value. Another option is to require the integral of $\phi$ in your domain to be a certain value (0 is a popular choice). | |
| Sep 23, 2022 at 6:53 | comment | added | Akhaim | While this is correct generally speaking, from a practical point of view a different set of Neumann boundary conditions,, e.g. G0,1=H0,1=1, converges. So, it seems that the way the Neumann boundaries are defined (or possibly their mean value) impacts the solution. After all, there is a bunch of commercial solvers which can resolve pure Neumann boundaries, and I would like to understand if there is a generic way of dealing with them. | |
| Sep 23, 2022 at 3:37 | history | answered | Maxim Umansky | CC BY-SA 4.0 |