9

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

Peculiar error when solving the Poisson equation

boyfarrell
  • 5,409
  • 3
  • 35
  • 67
  • 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 05:12
  • What does the reaction term look like? – Jan Sep 23 '13 at 05: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 09: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, http://scicomp.stackexchange.com/questions/8538/applying-dirichlet-boundary-conditions-to-the-poisson-equation-with-finite-volum – boyfarrell Sep 23 '13 at 12:42
  • @Jan I do not make any approximation to the integral of the source term, I just use the un-integrated value. I do this because in my discretisation I divide the flux by the cell volume to get a "concentration" per cell, this makes sense for the advection-diffusion-reaction equation from which the code is originally based. If we integrate the source term over the cell we get $\int_{\Omega_j} s(x) dx = s_j h_j$, dividing by the cell volume (in 1D just the cell width) yields $s_j$. If I increase the number of cells the problems remains, is this what you mean by insufficient sampling? – boyfarrell Sep 23 '13 at 12:47
  • Yes, however, I would call this quadrature using the midpoint rule. Do you divide every equation by the volume of the cell it is posed on? If you use this approximation of the rhs, the value you take to approximate $\rho$ in the "large" cell is ambiguous as it is directly at the jump. If you make the cell just a bit larger, you should take $0$ instead of $1$ and your model changes severely. However, if you refine the discretization in the region of the discontinuity, these effects should vanish. – Jan Sep 23 '13 at 13:04
  • Yes every equation is divided by the volume of the cell it is posed on. Yes I see, in the case where $\Omega_9$ is removed the cell straddles the discontinuity. In other cases the discontinuity is at the cell face so the integration is valid. Is there a general way around this problem in terms of computing the source term approximation in a different way? Or is the answer to always have a cell face where the source term has a discontinuity? The latter option seems a little bit unsatisfying. Ideally I would like the model to be stable for all source function profiles. – boyfarrell Sep 23 '13 at 14:16
  • Yes, and the general approach would be an adaptive refinement of the mesh. I do think that your scheme is stable. If you refine the mesh, this error should decrease. – Jan Sep 23 '13 at 14:20
  • Yes the magnitude of the error does decrease as more cells are added but even with 2000 cells it cause 5% error (approx.), so the effect is quite noticeable. I will make sure that any discontinuity of the source term happens on a cell face from now on. – boyfarrell Sep 23 '13 at 15:01
  • Stable for all source function profiles may be a little difficult given the way you impose your BCs. If you have the discontinuity in source at cell faces, that should alleviate a lot, but maybe not all cases. – Jesse Chan Sep 23 '13 at 17:50

2 Answers2

9

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.

Jesse Chan
  • 3,132
  • 1
  • 13
  • 17
  • 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 08:45
  • An relevant article, although I don't quite follow all the details, http://www.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 04 '13 at 10:16
  • Thanks! That's a problem, I don't honestly know. My immediate thought would be to steal from FEM DG methods and try a jump penalty (these restrict how discontinuous the cell averages can be). See slide 5 in http://maths.dur.ac.uk/~dma0mpj/summer_school/IPHO.pdf and set $\tau = 0$, this would be the basic interior penalty method. However, I don't know how to interpret this in terms of a FV method. The other approach would be to look more closely at Suli's paper and see if you can get a stable FV method from that (can't yet get it as it's behind a paywall at the moment). – Jesse Chan Oct 04 '13 at 12:14
  • Message me, I have a copy of the paper. This comment thread is getting very long, so maybe we should move to email anyway :) – boyfarrell Oct 05 '13 at 04:39
0

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?

Guido Kanschat
  • 1,124
  • 7
  • 11
  • 1
    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