JE37: Recovery schemes for variable coefficient diffusion

Note

The Maxima code used in this note is linked here. It requires use of Gkeyll CAS repo and Maxima setup as described in the documentation.

Introduction

Consider computing the weak-form of the term

\[g \equiv \frac{\partial f}{\partial t} = \frac{\partial}{\partial x}\left( \kappa \frac{\partial f}{\partial x} \right)\]

Here \(\kappa = \kappa(x)\) is the position dependent diffusion coefficient. This term appears, for example, in diffusion equation and Poisson equation with variable coefficient. We will at first assume \(\kappa\) is continuous and smooth. Of course, this does not mean that it’s DG representation is continuous.

Multiply \(g\) by a weight \(w\) and integrate by parts once over a single cell \(I_j = [x_{j+1/2}, x_{j-1/2}]\) to get the weak-form:

\[\int_{I_j} w g \thinspace dx = \left. \left(w \kappa \frac{\partial f}{\partial x} \right) \right\rvert_{j-1/2}^{j+1/2} - \int_{I_j} \frac{dw}{dx} \kappa \frac{\partial f}{\partial x} \thinspace dx.\]

From this weak-form we see that we have to compute interface values of \(f\), \(\kappa\) and \(f_x\). To do this we will recover \(\kappa\) across the interfaces \(x_{j\pm 1/2}\). We wil also recover \(f\) across interfaces \(x_{j\pm 1/2}\) for use in the surface terms We then have one of two choices for \(f\) that appears in the volume term, leading to two possible schemes:

1. Scheme 1: Use the DG representation of \(f\) in the volume term.

2. Scheme 2: Recover \(f\) across three cells \([x_{j-1}, x_j, x_{j+1}]\) (enforcing continuity of both \(f\) and \(f_x\) at \(x_{j\pm 1/2}\)) and use this in the volume term.

Note that in each of the above cases we will use the DG representation of \(\kappa\) for the volume term.

Before we use the schemes in solving the diffusion equation, we will see how each of the schemes behave for computing the second derivative, that is, the projection of \(g\) on basis functions. For the tests presented below we will choose

\[\begin{split}\kappa(x) &= 1 + \exp(-10 x^2) \\ f(x) &= 2 + \sin\left(2\pi\frac{x-0.5}{4} \right)\end{split}\]

on the domain \(-2\le x \le 2\).

Computing second derivatives: Scheme 1

Scheme 1, as described above, uses the DG exapansion of \(f\) in the volume term. It turns out that this scheme leads to an inconsistent discretization of the second-order derivative as the mean slope is mispredicted. See below for \(p=1\) case. This figure shows that though the cell-centers are converging correctly, the slopes are incorrect, and, in fact, do not converge as the grid is refined.

Recovery solution for \(p = 1\), \(N_x = 48\) (black) for \(g(x)\) (red) compared to the projection of the exact solution (sky blue). This figure shows that Scheme 1 is inconsistent, that is, the mean slope of \(g\) is mispredicted.

Higher-order DG recovery also shows the same behavior with Scheme 1: the higher (that cell-averages) coefficients are all systematically mispredicted.

Note that this misprediction is not a bug: one can compute the exact errors in the mean slopes and higher-moments by a careful analysis of the stencil for Scheme 1. Further, this issue also persists even with a constant diffusion coefficient.

Recovery solution for \(p = 2\), \(N_x = 10\) (black) for \(g(x)\) (red) compared to the projection of the exact solution (sky blue). This figure shows that Scheme 1 is inconsistent, that is, the mean slope and quadratic moment of \(g\) are mispredicted.

Computing second derivatives: Scheme 2

Convergence for Scheme 2 for \(g(x)\), \(p=0\)

\(N_x\)	\(l_2\)-error	Order
12	1.420
24	\(4.913\times 10^{-1}\)	1.53
48	\(1.407\times 10^{-1}\)	1.80

Recovery solution for \(p = 0\), \(N_x = 48\) (black) for \(g(x)\) (red) compared to the projection of the exact solution (sky blue).

For \(p > 0\) in the table below we also show the convergence of the cell-averages.

Convergence for Scheme 2 for \(g(x)\), \(p=1\)

\(N_x\)	\(l_2\)-error	Order	\(l_2\)-error (\(\overline{g}\))	Order (\(\overline{g}\))
12	\(2.367 \times 10^{-1}\)		\(1.883 \times 10^{-1}\)
24	\(3.352 \times 10^{-2}\)	2.82	\(1.639 \times 10^{-2}\)	3.52
48	\(5.137 \times 10^{-3}\)	2.71	\(1.227\times 10^{-3}\)	3.73

Recovery solution for \(p = 1\), \(N_x = 24\) (black) for \(g(x)\) (red) compared to the projection of the exact solution (sky blue). The exact and recovered solutions are visually indistinguishable.

Convergence for Scheme 2 for \(g(x)\), \(p=2\)

\(N_x\)	\(l_2\)-error	Order	\(l_2\)-error (\(\overline{g}\))	Order (\(\overline{g}\))
8	\(1.127 \times 10^{-1}\)		\(3.903 \times 10^{-2}\)
16	\(7.552 \times 10^{-3}\)	3.90	\(3.224 \times 10^{-3}\)	3.52
32	\(1.783 \times 10^{-4}\)	5.4	\(6.189\times 10^{-5}\)	5.7

Recovery solution for \(p = 2\), \(N_x = 10\) (black) for \(g(x)\) (red) compared to the projection of the exact solution (sky blue). The exact and recovered solutions are visually indistinguishable.

Conclusion

It seems that to get the proper approximation to \(g(x)\) need to use Scheme 2. However, it is likely that both Scheme 1 and Scheme 2 will work for solving parabolic equations, as the initial errors in slopes (and higher-moments) of the DG representation will be damped rapidly. However, this remains to be seen.