Poisson equation with discontinuous variable coefficient

Question

Trying to solve numerically the 2D Poisson equation with variable diffusion coefficient $K$, that in general can be discontinuous.

$$ \frac{\partial}{\partial x} ( K(x,y)\frac{\partial C}{\partial x}) + \frac{\partial}{\partial y} ( K(x,y)\frac{\partial C}{\partial y}) =a$$

where $a$ is a constant. Would this second order central difference approximation for the first LHS term (and a similar one for the second one) be the right choice? $$ \frac1h\Bigl(K(x+h,h)\frac{c(x+2h,y)-c(x+h,y)}{h}-K(x,y)\frac{c(x+h,y)-c(x,y)}{h}\Bigr). $$ I assumed the same grid size $h$ in $x$ and $y$. But where would you take the $K$ coefficient? E.g. for $K(x,y)\frac{\partial C}{\partial x}$, discretized using the nodes in $(x+h,y)$ and $(x,y)$, would you take $K(x,y)$ or $K(x+h,y)$? Seems to me that it would be wise to take it at the location where $C$ is higher to mimic some sort of upwind approach for the diffusive flux, but not sure if it makes any sense in an elliptical PDE. Does it even matter? E.g.

Answer 1

First, some comments on discretization schemes for elliptic problems in general. There are some improvements you can make on the difference scheme that you've derived. The idealized problem is to solve a linear system with a self-adjoint, positive-definite linear operator. Your discretization has the right order of accuracy but fails to be self-adjoint. As Wolfgang alluded to in the comments, the usual approach is to evaluate the diffusion coefficient at the midpoint. Let's just look at a 1D problem to simplify things; you can think of this as evaluating the fluxes at cell midpoints: $$F(x + h/2) = K\cdot\frac{dC}{dx}(x + h / 2) \approx K(x + h / 2)\frac{C(x + h) - C(x)}{h}$$ and then evaluating the divergence at the midpoint of midpoints, i.e. back on the nodes: $$\nabla\cdot F(x) \approx \frac{F(x + h / 2) - F(x - h/2)}{h}.$$ For 2D problems you can play the same game in the $y$ direction as well. You can find more information about this in $\S$2.15 of Leveque's book.

But you don't just have any old elliptic problem; you've specified that the diffusion coefficient is discontinuous, and in this case a few things break down. First, it's rare that a PDE is actually a first principle. They really start life as a conservation law:

$$\forall \omega \subset \Omega, \quad \int_{\partial\omega}K\nabla C\cdot\nu\,ds = \int_\omega a\,dx$$

where $\nu$ is the unit outward normal vector to $\partial\omega$, or a variational form:

$$\forall v \in H^1(\Omega), \quad \int_\Omega K\nabla C\cdot\nabla v\,dx = \int_\Omega a\cdot v\, dx.$$

To derive a PDE, you need to apply the divergence theorem one way or the other, which assumes everything is differentiable. PDEs are convenient lies that we tell ourselves in order to write down analytical solutions in nice geometries with constant coefficients. For this reason, I don't recommend using the manipulations in Željko_JL's answer -- those only make sense when $K$ is differentiable.

Likewise, a lot of the convergence theory for finite difference methods relies on the same assumptions of continuity and differentiability. Now it's possible that a finite difference method will still converge, but it's typical to lose half an order of convergence. For example, the difference scheme I wrote down above has relative $L^2$-norm error that goes like $\mathscr{O}(h)$ for nice problems, but instead converges like $\mathscr{O}(h^{1/2})$ when the diffusion coefficient is discontinuous. This is a pretty dramatic loss of efficiency.

To work around this loss of accuracy, you have two options. First, there are some hacks that make the finite difference method yield higher-order convergence for problems with discontinuous coefficients. You can find this in Li and Ito. Second, you can try to align the grid exactly with the discontinuities in $K$. If these discontinuities are geometrically simple you might still be able to use the finite difference method. If they aren't, you may have an easier time using an unstructured mesh and switching to using the finite element method.

The classic thing to test when the diffusion coefficient is discontinuous is that the derivative has a jump at the interface:

$$K_1\nabla C_1\cdot \nu_1 + K_2\nabla C_2\cdot\nu_2 = 0$$

which is again something you can derive from the conservation law or the variational form. To stress test your solver, in order of difficulty, try a problem where (1) the jumps in $K$ are aligned with the coordinate axes and they fall on the nodes of your finite difference grid, (2) the jumps are still aligned with the coordinate axes but they fall off the grid, (3) the jumps fall on straight interfaces that are not aligned with the coordinate axes, and (4) the jumps are on curved interfaces.

Answer 2

I am not sure if what you use is the central difference. For example, in this post a discretization scheme is provided for the equation of the type $$\nabla(\varepsilon(x,y)\nabla\varphi)=f(x,y).$$ Your case demands that $f=a$, and $\varepsilon=K$, where both $\varepsilon$ and $K$ are functions of $x$ and $y$. In said post central finite differences are used for both second and first derivatives in the form (adapted to your 2D notation, with general assumption that $h$ differs for $x$ and $y$): $$\nabla(K\nabla\varphi) = \nabla K\nabla\varphi+K\nabla^2\varphi\to\frac{(K_{i+1,j}-K_{i-1,j})(\varphi_{i+1,j}-\varphi_{i-1,j})}{4h_x^2}+\frac{(K_{i,j+1}-K_{i,j-1})(\varphi_{i,j+1}-\varphi_{i,j-1})}{4h_y^2}+K_{i,j}\frac{\varphi_{i+1,j}-2\varphi_{i,j}+\varphi_{i-1,j}}{4h_x^2}+K_{i,j}\frac{\varphi_{i,j+1}-2\varphi_{i,j}+\varphi_{i,j-1}}{4h_y^2}.$$ Then for the term $K\frac{\partial^2\varphi}{\partial x^2}$ one can take $K_{i,j}$, and for the term $\frac{\partial K}{\partial x}\frac{\partial\varphi}{\partial x}$ one should use both $K_{i-1, j}$ and $K_{i+1, j}$ which will multiply the discretized $\frac{\partial\varphi}{\partial x}$. Similar approach is then done in $y$ direction.

Hope this helps.