Applying Dirichlet boundary conditions to the Poisson equation with finite volume method

I would like to know how Dirichlet conditions are normally applied when using the finite volume method on a cell-centered non-uniform grid,

My current implementation simply imposes the boundary condition my fixing the value of the first cell,

$$\phi_1 = g_D(x_L)$$

where $\phi$ is the solution variable and $g_D(x_L)$ is the Dirichlet boundary condition value at the l.h.s. of the domain (NB $x_L \equiv x_{1/2}$). However this is incorrect because the boundary condition should fix the value of the cell face not the value of the cell itself. What I should really apply is,

$$\phi_{L} = g_D(x_L)$$

For example, lets solve the Poisson equation,

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

with initial condition and boundary conditions,

$$\rho=-1\\ g_D(x_L)=0 \\ g_N(x_R)=0$$

(where $g_N(x_R)$ is a Neumann boundary condition on the right hand side).

Notice how the numerical solution has fixed the value of the cell variable to the boundary condition value ($g_D(x_L)=0$) at the left hand side. This has the affect of shifting the whole solution upwards. The effect can be minimized by using a large number of mesh points but that is not a good solution to the problem.

Question

In what ways are Dirichlet boundary conditions applied when using the finite volume method? I assume I need to fix the value of $\phi_1$ by interpolating or extrapolating using $\phi_0$ (a ghost point) or $\phi_2$ such that the straight line going through these points has the desired value at $x_L$. Can you provide any guidance or an example of how to do this for a non-uniform cell-centered mesh?

Update

Here is my attempt at using a ghost cell approach you suggested, does it look reasonable?

The equation for cell $\Omega_1$ is (where $\mathcal{F}$ represents the flux of $\phi$),

$$\mathcal{F}_{3/2} - \mathcal{F}_{L} = \bar{\rho}$$

We need to write $\mathcal{F}_{L}$ in terms of the boundary condition using a ghost cell $\Omega_0$,

$$\mathcal{F}_{L} = \frac{\phi_1 - \phi_0}{h_{-}} \quad\quad \text{[1]}$$

But we ultimately need to eliminate the $\phi_0$ term from the equation. To do this we write a second equation which is the linear interpolation from the centre of cell $\Omega_0$ to the centre of cell $\Omega_1$. Conveniently this line passes through $x_L$, so this is how the Dirichlet conditions enters the discretistaion (because the value at this point is just $g_D(x_L)$),

$$g_D(x_L) = \frac{h_1}{2h_{-}}\phi_0 + \frac{h_0}{2h_{-}}\phi_1 \quad\quad \text{[2]}$$

Combining equations 1 and 2 we can eliminate $\phi_0$ and find an expression for $\mathcal{F}_L$ in terms of $\phi_1$ and $g_D(x_L)$,

$$\mathcal{F}_L = \frac{1}{h_{{-}}} \left(\phi_{1} - \frac{1}{h_{1}} \left(2 g_{D} h_{{-}} - h_{1} \phi_{1}\right)\right)$$

Assuming that we are free to choose the volume of the ghost cell we can set $h_0 \rightarrow h_1$ to give,

$$\mathcal{F}_L = -\frac{2 g_{D}}{h_{1}} + \frac{2 \phi_{1}}{h_{{-}}}$$

This can be simplified further because if the cells $\Omega_0$ and $\Omega_1$ are the same volume then we can set $h_{-}\rightarrow h_1$ finally giving,

$$\mathcal{F}_L = \frac{2}{h_{1}} \left( \phi_1 - g_D \right)$$

However, this approach has recovered the definition that is unstable so I'm not too sure how to proceed? Did I interpret your advice incorrectly (@Jan)? The strange thing is that is seems to work, see below,

See below, it works,

• Right, your derivation is correct. And it really resembles what I've called (**) in my answer. And, thus, it is proven to be stable. I will add a comment in my answer. – Jan Sep 21 '13 at 9:23
• Also, as a general remark, stability results are typically sufficient conditions. I.e. if a scheme doesn't meet the conditions, in some situations, it may well produce reliable results. – Jan Sep 21 '13 at 9:28

In the stability-analysis of FVM discretizations for elliptic problems with Dirichlet BC, a central assumption is that the inner cells, where you state the PDE, have no intersection with the boundary, i.e. $$\overline \Omega_i \cap \Gamma_D = 0 \quad \quad (*)$$ if seen as a set in $\mathbb R^{n-1}$ if your domain $\Omega \subset \mathbb R^n$, cf., e.g., the book by [Grossmann&Roos, p. 92]

Thus, if in your setup, the approach $$\Bigl(\frac{d\phi}{dx}\Bigr)_{1/2} = \frac{2}{h_1}\Bigl(\phi_1-\phi_{1/2}\Bigr) \quad \quad (**)$$ is unstable, this is no in contradiction to known stability results. EDIT: Using a ghost cell and linear interpolation into it, for a particular choice of volume and distance, one obtains $(**)$ as the flux. Thus, $(**)$ is indeed a stable scheme.

Stability and convergence (of first order in the discrete max-norm) for the poisson problem has been proven by Grossmann&Roos for grids, with distinct boundary cells with their "centers" on the actual boundary as illustrated in my drawing for an 1D case.

Here, the differential quotient on the interface is approximated in a straight-forward manner.

I would say ghost cells are the common approach, because of two reasons.

• They mimick the stable situation described in my drawing but with an interpolated boundary condition
• They are simply attached to the physical boundary. Thus, one can use a triangulation of the domain, what is also advantageous, since one often has also natural BCs that are directly imposed on the interface [Grossmann&Roos, p. 101].

So, I suggest you use ghost cells for the Dirichlet boundary. In your example, this will be adding $\phi_0$ to your system and the condition that an interpolant between $\phi_0$, $\phi_1$, and maybe others equals $g_D$ at the boundary.

• Thank you Jan, that's really interesting. That would certainly mimic my experience with certain approaches being unstable. Am I right, if I use a ghost cell approach I don't need to shift the last cell so that the centre is on the boundary? I also have a problem with the concept of shifting the boundary cell; doesn't it imply that that cell has zero volume? – boyfarrell Sep 20 '13 at 14:42
• Yes, if you introduce a ghost cell, then you need not change the grid of your example picture. Regarding the shift you mentioned to establish the situation of my drawing. No, it is not a degenerated cell! The offset $h_\Gamma$ really enters the equations in so far as this strip does not appear in the integrals, taken, e.g., of the right hand side. – Jan Sep 20 '13 at 14:56
• Would be interesting to see what happens as $h_\Gamma \to 0$, as this limit is the discretization approach (**). (I have assumed you use the linear interpolation between $\phi_1$ and $\phi_0$.) – Jan Sep 20 '13 at 14:58
• Can the dependence on the value of the ghost cell be removed with this approach? I guess is must not be included into the equations but only used a tool to write the boundary conditions. Regarding the "shifted" boundary cell. It looks like that point uses finite difference rather than the finite volume method. Would that be accurate? – boyfarrell Sep 20 '13 at 15:22
• OK I understand! Thank you. There is a typo. in the 2nd paragraph "Thus, if in your setup, the approach [eqn] is unstable, this is no contradiction to known stability results." The "no" should be "in". This flips the meaning of the sentence to mean the opposite of what you want (I think)! – boyfarrell Sep 21 '13 at 9:44

Your "straight line" approach would mean that $\phi_1-\frac{\phi_2-\phi_1}{x_2-x_1}(x_1-x_0)=0$, where $x_0$ is the location of the boundary and $x_i$ are the locations of the places where you define the $\phi_i$. This means that you can eliminate $\phi_1$ in terms of $\phi_2$, in the same way as in your first approach you eliminated $\phi_1$ right away by setting it to zero.

What you're finding here is why finite volumes aren't frequently used for the elliptic equations for which one poses Dirichlet conditions. They are used for conservation laws where the more natural conditions are stated in terms of fluxes.

Let's say that your finite-volume form of the Poisson equation $$\frac{d^2 \phi}{dx^2} = f$$ can be written as $$\Bigl(\frac{d\phi}{dx}\Bigr)_{3/2}-\Bigl(\frac{d\phi}{dx}\Bigr)_{1/2} = \int_{x_{1/2}}^{x_{3/2}}f\, dx$$ Of course, you can write the usual approximation $$\Bigl(\frac{d\phi}{dx}\Bigr)_{3/2} = \frac{\phi_2-\phi_1}{h_+}$$ There are some subtleties related to non-uniform grids that I'm ignoring here.

Now the question is about how to you approximate $(d\phi/dx)_{1/2}$ such that you can incorporate the b.c. $\phi_{1/2}$. One possibility is to change your approximation stencil compared to the usual stencils so that it includes the points $x_{1/2}$, $x_1$ and $x_2$. If my memory serves me correctly, that approximation is (assuming a uniform grid with spacing $h$, so you will need to extend this) $$\Bigl(\frac{d\phi}{dx}\Bigr)_{1/2} = \frac{1}{h}\Bigl(-\frac{1}{3}\phi_2+3\phi_1-\frac{8}{3}\phi_{1/2}\Bigr)$$ but this should be checked. In that approximation, which is second-order accurate, you can incorporate your b.c. Of course, an even simpler and only first order accurate approximation would be $$\Bigl(\frac{d\phi}{dx}\Bigr)_{1/2} = \frac{2}{h_1}\Bigl(\phi_1-\phi_{1/2}\Bigr)$$ In this way, the Dirichlet boundary condition is only enforced "weakly" (meaning through a flux). As already commented by Wolfgang, this is one of (the?) reason why finite volume methods are not used for elliptic problems as much as finite element methods.

Of course, one thing that also needs to be checked is the stability of your discretization with the second order approximation at the boundary. Off the top of my head, I don't know whether it will be stable combined with a centered second order approximation in the interior. A matrix stability analysis will tell you for sure. (I'm virtually certain that the first order approximation at the boundary will be stable.)

You mention the possibility of using ghost points. This leads to the problem that you need to extrapolate from the interior into the ghost point and use the b.c. in the process. I suspect, but have not "proved" it, that at least some ghost point treatments are equivalent to using the kind of approach I've outlined above.

Hope this helps a little bit.

• Hello Brian. I didn't think it was possible to apply Dirichlet boundary conditions using the flux form (i.e. weakly). In fact I asked that question a few months ago, scicomp.stackexchange.com/questions/7777/… I tried to implement something like this back then but, for what ever reason, the implementation was unstable and always failed. Do you know a reference in which Dirichlet conditions are applied to Poisson equation, I am interested to know what is standard? Maybe this is not done for elliptic equations? – boyfarrell Sep 20 '13 at 2:02
• I don't know of a standard, but I cannot imagine that all such implementations are unstable. Did you try the matrix analysis? It should be very simple to carry out in this case. People do solve the Navier-Stokes equations with ghost-point treatments and treatments like the one above. (Of course, there viscous effects do not dominate to such an extent that you can regard the Poisson equation as a good model.) Perhaps these references help: ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/… and nas.nasa.gov/assets/pdf/techreports/1997/nas-97-011.pdf – Brian Zatapatique Sep 20 '13 at 6:22
• Hello Brian. No I didn't try matrix analysis. To be honest I'm not too sure how to do that. I will have time next week to revisit this problem so I may post a new question then! – boyfarrell Sep 20 '13 at 6:38
• My understanding is also that ghost point (quadratic) extrapolation ends up being equivalent to the classic Shortley-Weller finite difference discretization for irregular (curved) Dirichlet boundary conditions, e.g. as described on p74 of Morton and Meyers Numerical Solution of Partial Differential Equations (2nd edition). (The linear extrapolation version is equivalent to the simpler method of Gibou et al. sciencedirect.com/science/article/pii/S0021999101969773) Also: both linear and quadratic extrapolants give 2nd order accurate solutions, but linear only 1st order gradients. – batty Nov 4 '15 at 1:37