4

I have written Poisson solvers using two different methods:

A classic Jacobi scheme and one using the multigrid solver Hypre. I made up a couple of test cases ensuring the validity of those solvers.

For both cases the domain is defined as $(x,y) \in [-1,1]^2$ with periodic boundary conditions. Also, the grids first and last point are the same:

$$p(0,y) = p(N_x-1,y)$$

$$p(x,0) = p(x,N_y-1)$$

Test case 1

  • $f_{rhs} = -8 \pi^2 \sin(2\pi x) \sin(2\pi y)$
  • $p_{exact} = \sin(2\pi x) \sin(2\pi y)$

For both solvers, the solution is 2$^{\mathrm{nd}}$ order accurate is space. No problems here so far.

Test Case 2

  • $p_{rhs} = e^{-10 (x^2 + y^2)}$
  • $p_{exact}:$ No analytical solution, and therefore the numerical solution is differentiated using a high-order compact scheme and compared to $p_{rhs}$

Note that in this case, $\int_V p_{rhs}dV \neq 0$ and therefore the problem is ill defined. Therefore, the $rhs$ must be modified:

$p_{rhs} = e^{-10 (x^2 + y^2)} - \dfrac{\int_D e^{-10 (x^2 + y^2)} dx dy}{V}$

where $V$ is the domain volume. The integral is computed using the Trapezoidal rule.

This is where things get tricky. No matter how fine my grid is, $\left(p_{num}\right)_{xx} + \left(p_{num}\right)_{yy}$ never converge to $p_{rhs}$. When the grid is fine enough, the solvers converge, but the solution is off by ~20% while the overall profiles are relatively correct. When the grid is coarse, the Hypre solver simply diverges.

Question

Have I missed anything? Is my approach inconsistent/wrong?

Community
  • 1
solalito
  • 337
  • 2
  • 11
  • 2
    You say that the solution is unique up to a constant, but I think there exist no solution. But if you correct it to be zero on average, it should exist (and be determined up to a constant). – Vladimir F Героям слава Aug 09 '17 at 15:14
  • @VladimirF Yes you are right. Fixed. – solalito Aug 09 '17 at 15:40
  • 1
    I am not sure you can use the trapezoidal integration to correct the p_rhs. You need the RHS of the discretized linear system to be zero on average. Is it? – Vladimir F Героям слава Aug 10 '17 at 13:36
  • 1
    Fix one node, could be anywhere, for example on corner, and set value at that node to your exact solution, you will get well posed problem. – likask Aug 10 '17 at 19:51
  • Please skip the trapezoidal rule and compute the integral by yourself to the exact value. There is no need to introduce an additional approximation error. Hint: the integral is the product of two one-dimensional integrals. – shuhalo Aug 11 '17 at 17:23
  • @shuhalo That is of no interest to me as I need to integrate the solver to a CFD code where there is no analytical solution of the RHS, hence the trapezoidal rule. – solalito Aug 14 '17 at 07:45
  • @likask You put me on the right track. But fixing one node is not enough, I had to actually impose dirichlet boundary everywhere to a fixed constant. The value of that constant did not matter, the derivatives of the solution were the same in all cases. – solalito Aug 14 '17 at 07:46
  • I’m a bit confused with the rhs subscripts, which one is on the right hand side? The solution you provided may be because exp(...) is NOT periodic, and so you cannot apply periodic BCs. – Charles Feb 10 '18 at 17:17
  • @solalito did you perhaps check this link. This user had a problem with periodic boundaries and ill-defined problem. – Akhaim Dec 21 '22 at 08:01

2 Answers2

1

So after much testing, I have found a solution to this problem. However, I am not quite sure of the mathematical justification for it.

Basically, when using purely periodic BC, iterative solvers would not converge and the solution will keep growing. By fixing Dirichlet BC on all sides of the domain to a fixed constant $C \in \mathbb{R}$, on top of using periodicity, the RHS term can be obtained by differentiating the solution with an order of accuracy of 2 (note that high-order scheme are using for numerical differentiation to ensure the leading error term comes from the Poisson solver).

Furthermore, this got me thinking that there was no need for periodic BC to begin with since only the gradient of the solution is of interest to me:$\; \dfrac{\partial p}{\partial x_i}$. By just enforcing dirichlet BC on all sides, this gradient can be retrieved with the desire order of accuracy. I will be happy to hear more of why that is though.

solalito
  • 337
  • 2
  • 11
0

In your Test Case 1, the solution is periodic. Hence, you had no problems when you imposed periodic boundary conditions.

For the Test Case 2, the solution may not be periodic (I suspect that it is not). By enforcing the periodic conditions, you get some solution which does not actually corresponding to the source term. That is why your solution is off by 20%. This should be due to incompatibility between the actual solution and the obtained solution at the boundaries. If you check carefully, you may find that the maximum error occurs at the boundaries.

You can verify your code with another benchmark example from FEniCS. https://fenicsproject.org/olddocs/dolfin/1.3.0/python/demo/documented/periodic/python/documentation.html

Chenna K
  • 934
  • 5
  • 12
  • Why shouldn't there exist a periodic solution? – shuhalo Aug 11 '17 at 17:15
  • 1
    In order to answer further, I need to know the value you impose on the boundary for the Test Case 2. – Chenna K Aug 11 '17 at 21:56
  • There are periodic boundary conditions. Nothing to impose on the boundary. – shuhalo Aug 12 '17 at 04:07
  • Even for periodic boundaries, you have to specify the boundary values on one of the boundaries in each direction. However, I don't think it really matters whether you use periodic boundary conditions or not. As the Jacobi Poisson solver is used here, it is assumed that the values of all the boundary points are known. If you don't specify anything, then the value is assumed to be zero. – Chenna K Aug 12 '17 at 14:08
  • I tested these two problems using Jacobi solver. For 40x40 grid, Case 1 requires 167 iterations, while Case 2 takes 2674 iterations (~16 times more iterations). I suspect that this is due to incompatible boundary conditions. Even if the periodic solution exists for Case 2, it might not vanish on the boundaries of the domain considered for this problem. So, by enforcing periodic bcs with zero value everywhere on the boundary, you are asking the solver to solve for some unnatural (to the problem) conditions. – Chenna K Aug 12 '17 at 14:12
  • I have never asked for periodic BC with zero value on the boundary. How do solve a Poisson with periodic BC on a line segment? What do you on a circle? – shuhalo Aug 12 '17 at 17:02
  • You haven't. This is not your question. I was providing more insights on the problem. I was trying to explain what could have possibly gone wrong. – Chenna K Aug 13 '17 at 22:18
  • @Shuhalo, To answer question, I solve the Poisson equation with periodic BC like I solve any other PDE with periodic BC. I remove dependent nodes and solve for independent nodes. If I know the value of the field on the boundary, then I use it as Dirichlet BC. If not, I solve for them. One can not apply arbitrary BC without physical significance and without having consequences. Periodic BC can be applied only when the domain and BCs are periodic (repetitive). On a circle, you have either a Dirichlet or a Neumann or a mixed BC. – Chenna K Aug 13 '17 at 22:31
  • Just a comment but I suggest that you'd might want to try to use a direct solver to check if it's the iteretive solver's fault or of the problem is ill defined in some way – CuteCompute Aug 24 '21 at 15:51