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

Question

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,

Answer 1

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.

Answer 2

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.

Answer 3

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.