Smoothing the diffusion coefficient to improve convergence

Question

I have been reading a book by Thomee and he considers the case of $u_t=(au_x)_x$, for the case of $a$ possibly being discontinuous. Then he says that the problems with convergence might occur, and thus to fix that one can smooth the coefficient $a$ and restore convergence rate.

Please correct me if I am wrong but discontinuous coefficients affect the smoothness of the solution of a pde at the first place, and therefore affects the local truncation error of a numerical method. How does smoothing of the coefficient affects convergence if the convergence depends on local truncation error, which depends on the properties of the function itself and its derivatives by the original stated problem which we can't control, because it is an input to the problem?

Answer 1

I don't know what exactly Thomee had in mind, but if you can't get solution at all for a discontinuous coefficient and a particular discretization scheme, you can smooth the coefficient on a lengthscale proportional to the mesh size $h$ and obtain a problem that's solvable. Then you make the smoothing distance go to zero as you refine the mesh and you obtain a sequence of numerical solutions that approximate a sequence of exact solutions of perturbed problems that converge to the solution of the original problem. That's a common trick.

Answer 2

Strong form discretizations are extremely cumbersome for the case of discontinuous coefficients. It is better to use weak form methods that do not require differentiating the solution. Furthermore while the product $a u_x$ is continuous even for discontinuous coefficients, $u_x$ alone is not. Therefore even evaluating one derivative of the solution is something for a numerical method to be very careful with. The standard (robust) approaches are $H(\textrm{div})$ mixed finite elements, mimetic finite differences, and primal finite volume/difference methods with a choice of quadrature making them equivalent to certain mixed methods.

I started my reply with the note above because a numerical discretization must be stable and conserve the correct properties for it to be worth trying to solve the resulting algebraic system. Make sure the spatial discretization is robust for the chosen purpose before proceeding further. Local conservation, for example, is critical for any problem exhibiting non-smoothness. Also confirm that your quadrature respects the non-smoothness.

When the coefficient has discontinuities, but when the jumps conform to the grid, you can use discrete spaces that have a discontinuous gradient at the interface, thus allowing accurate approximation without elaborate tricks. This is most natural in the finite element context, which partly explains its popularity for structural mechanics.

For the case of the coefficient $a$ being rough at a subgrid scale, there are numerous homogenization techniques. For problems with scale separation in 1D, the simple approach of replacing the rough coefficient with its harmonic mean on each element is optimal. For multiple dimensions or for problems that do not have separated scales, more sophisticated techniques are needed.

In all cases, simple smoothing of the coefficients degrades accuracy significantly, and you will have no hope of accurate large-scale solutions until the mesh resolves the fine scale structure and the scheme will generally not be convergent in the energy norm.