# Peculiar error when solving the Poisson equation on a non-uniform mesh (1D only) finite volume method

I have been trying to debug this error the last few days I wondered if anybody has advice on how to proceed.

I am solving the Poisson equation for a step charge distribution (a common problem in electrostatics/semiconductor physics) on a non-uniform finite volume mesh where the unknown are defined on cell centres and the fluxes on the cell faces.

$$0 = (\phi_x)_x + \rho(x)$$

the charge profile (the source term) is given by,

$$\rho(x)= \begin{cases} -1,& \text{if } -1 \leq x \leq 0\\ 1,& \text{if } 0 \leq x \leq 1\\ 0, & \text{otherwise} \end{cases}$$

and the boundary conditions are,

$$\phi(x_L)=0 \\ \frac{\partial\phi}{\partial x}\bigg|_{x_R}=0$$

and the domain is $[-10,10]$.

I am using code developed to solve the advection-diffusion-reaction equation (I have written myself see my notes here, http://danieljfarrell.github.io/FVM). The advection-diffusion-reaction equation is a more general case of the Poisson equation. Indeed the Poisson equation can be recovered by setting the advection velocity to zero and removing the transient term.

The code has been tested against a number of situations for uniform, nonuniform and random grids and always produces a reasonable solutions (http://danieljfarrell.github.io/FVM/examples.html) for the advection-diffusion-reaction equation.

To show where the code breaks down I have made the following example. I setup a uniform mesh of 20 cells and then make it nonuniform by removing a single cell. In the left figure I have removed cell $\Omega_8$ and in the right $\Omega_9$ has been removed. The 9th cell covers the region where the source term (i.e. the charge) changes sign. The bug appears when the grid is nonuniform in a region where the reaction term changes sign. As you can see below.

Any ideas what could possibility be causing this issue? Let me know if more information regarding the discretisation would be helpful (I didn't want to pack too much detail into this question).

• Can you specify how you impose the Dirichlet condition at $x=0$, and what you mean by $\rho = -1$ as an initial condition (isn't the equation you specified steady state)? – Jesse Chan Sep 23 '13 at 5:12
• What does the reaction term look like? – Jan Sep 23 '13 at 5:43
• What scheme do you use to approximate the integrals of the source term? This behavior might also be caused by insufficient sampling of the source. (What, probably, is the same mechanism mentioned in @JLC 's answer.) – Jan Sep 23 '13 at 9:39
• I have updated the question to use standard terminology. I have a source term ($\rho$) not a reaction term because as you pointed out we just need the steady-state value. The correct spatial dependence of $\rho$ is now given (initial value was incorrect). – boyfarrell Sep 23 '13 at 12:26
• @JLC the Dirichlet BCs are imposed using a ghost cell approach (my notes online are out of date regarding this implementation detail), see here for how I do it, scicomp.stackexchange.com/questions/8538/… – boyfarrell Sep 23 '13 at 12:42

Just as an aside, your github documentation is fantastic.

This is just a guess from DG methods, which can have similar issues if numerical fluxes aren't chosen carefully (I figure FV methods are a subset of DG methods). If you're using interpolation from cell centers to define your fluxes, then this should be equivalent to using the average as a numerical flux in DG and using piecewise constant basis. For standard DG methods for Poisson, this leads to numerically non-unique solutions - you can get a non-trivial null space for the discrete operator, which I think is what's causing your issues in the 2nd example. See this DG paper for their theory on it from the DG side.

I'll try to mock up an example for FV which shows how this comes into play.

Edit: so here's a small example of what's going on. Consider cells 1-9 and 11-20 in which $\rho(x) = 0$. From the right side (11-20), we have $f(x_{20}) = 0$ due to the Neumann condition, which tells us from conservation for that cell that $f(x_{19}) = \ldots = f(x_{11}) = 0$. Since the flux is the average of cell values, this tells us that $\phi(x)$ is constant over all these cells.

From the left side (1-9), we have $f(x_{i+1})-f(x_i) = 0$. If $f(-10) = 0$ and we use ghost cells, then $f(-10) = \phi_{9.5} - \phi_{\rm ghost} = \phi_{9.5}$. Conservation over the next few cells gives that $f(x_i) = f(-10) = \phi_{9.5}$ (i.e. constant slope). However, note that this can be any slope, just constant.

The issue comes about in the middle cell. Like Jan mentioned, you undersample the forcing in the second mesh. This throws off the balance equations at that point, gives you an error in $f(10)$, which then propagates backwards and messes up both the slope in the left half of the domain as well as the value of $\phi(9.5)$.

This sensitivity to errors in forcing is what's problematic - unlike FEM or FD methods which explicitly enforce the Dirchlet condition at $x=-10$, FV enforces it weakly using ghost nodes. Intuitively, ghost node weak imposition is like setting a Neumann condition at your left boundary as well. If you have two Neumann conditions for a diffusion problem, your problem is ill posed and has a nonunique solution (you can add any constant to that problem and still have a solution). You don't quite get that at the discrete level here, but you do get very sensitive and mesh-dependent behavior as you see in your experiments.

• From experimenting I can demonstrate that the FVM method is only stable when the cells either side of the discontinuity (sign change) of the source function have equal volumes. Would your analysis agree with this? This means that I must pay more attention to generating a sensible grid of my problems that I have done before. Maybe I should consider learning the FEM method next? – boyfarrell Sep 24 '13 at 8:45
• An relevant article, although I don't quite follow all the details, jstor.org/discover/10.2307/2157873 – boyfarrell Sep 24 '13 at 14:22
• The FVM method is only stable in this case when the grid is aligned with the source function somehow. If your source function changes, then you'll have to tune your grid again. I don't think generating a sensible grid is the right approach to this problem - you have an unstable method. – Jesse Chan Sep 24 '13 at 14:23
• That's a good find. Suli's a solid analyst. I would say learning FEM might be fun, but FD should also work for any elliptic 1D problems. You might also see what FV folks do (maybe augment their fluxes with penalty terms) to get convergence for 2nd order elliptic problems on general grids. Mathematical folk wisdom usually says that FV/upwinded FD is great for hyperbolic problems while FEM/central diff FD is great for elliptic. – Jesse Chan Sep 24 '13 at 14:26
• I'm revising this problem. Re-reading your answer I must say it's fantastic! I see your point that the method should change because that is the root of the problem (not the grid). Do you have any suggestions or things I could follow (that are acceable to a non expert) on how to better approximate the flux in this case. I.e in a way that might make it more stable. If possible I would like to find a better FVM for this equation. – boyfarrell Oct 4 '13 at 10:16

The first thing to notice are your boundary conditions. Since you can change the slope and the value, you have neither Dirichlet, nor Neumann conditions.

Then, every straight line is a solution where the right hand side is zero. You got that part.

Your fluxes are probably dependent on $h$. Do you use the correct $h$ where you eliminate a cell?

• No, that's not correct. The problem is well posed. For the case $\rho \equiv 0$ only $\phi \equiv 0$ is a solution, there is no other linear function that is zero at one point and has zero slope at a second point. – Jan Sep 23 '13 at 19:06