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

Question

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).

Answer 1

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.

Answer 2

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?