Manufactured solution for $-\operatorname{div}(a(x) \nabla{u}) = f$ when $\alpha(x)$ is discontinuous

Question

I'm studying the dealii tutorial number 4,5 and I understand the workflow. I've also been able to find the EOC by using manufactured solution where $f$ is a smooth r.h.s. and $\alpha(x)$ smooth too.

Now I'm interested in checking the EOC when $\alpha(x)$ is discontinuous (link to Doxygen docs of dealii)

Let's say I have $\Omega = [-1,1]\times[-1,1]$ and I assume $u(x,y)=\sin(\pi x ) \sin(\pi y)$ to be the solution of

$$- \nabla \cdot (\alpha(x) \nabla{u}) = f$$ with homogeneous Dirichlet b.cs and $\alpha(x)$ defined as in the link. The forcing term I obtain is:

$$f(x,y)= \begin{cases} 20 \cdot 2\pi^2 \sin(\pi x)\cos(\pi y) \\ 2 \pi^2 \sin(\pi x)\cos(\pi y) \end{cases}$$

Looking at the EOC, it seems that my numerical solution is not converging to the true solution $u(x,y)$ (I obtain an EOC which decreases from 0.23 to 0.02...) So I'm wondering if this approach is not okay.

• How can I construct a manufactured solution when the $\alpha(x)$ is not continuous?

---

EDIT:

Using that $\alpha(x)$, I could set $\nabla{u}=[1,20]^T$ on $B_{1/2}(0,0)$ and $\nabla{u}=[20,40]^T$ outside $B_{1/2}(0,0)$.

This implies that $\alpha(x,y) \nabla{u}$ is continuous, actually it's constant, hence I obtain $f=0$. Also, $$u(x,y)= x+20y$$ in $B_{1/2}(0,0)$, and $$u(x,y)=20x+40y$$ outside $B_{1/2}(0,0)$. With this argument, I always have a jump in $\nabla{u}$ due to $\alpha(x,y)$.

However, I'm still facing issues in showing the convergence, and I think this may be due to the fact that I have $f=0$ as forcing term.

---

EDIT^2:

Following the comment from @MaximUmansky, let $\alpha \nabla{u} = \vec{f}$ such that $$-\operatorname{div}(\vec{f}) = f_x + f_y= \sin(x)+\cos(y)$$

so $\vec{f}=[\cos(x), - \sin(y)]^T$

With $\alpha$ as above, we need $\nabla{u} = [\frac{\cos(x)}{20}, -\frac{\sin(y)}{20}]$ in $B_{1/2}(0,0)$

and $\nabla{u} = [\cos(x),-\sin(y)]$ outside $B_{1/2}(0,0)$.

So I have, in $B_{1/2}(0,0)$, $$u(x,y)=\frac{1}{20}(\sin(x)+\cos(y))$$

and $$u(x,y)=\sin(x) +\cos(y)$$ outside $B_{1/2}(0,0)$.

Answer 1

@cfdlab's answer gives a general way to construct solutions for discontinuous coefficients, but here is a slightly more theoretical perspective to it as well. All of this can be understood in 1d, so imagine all of the functions below to be functions of just one argument $x$ for a moment.

First, you can't just choose any $u$ and $\alpha$ to obtain a function $f(x)$. That's because in general, if $u$ is continuous and $\alpha$ is discontinuous, then $\alpha \nabla u = \alpha\frac{du}{dx}$ is a discontinuous function and if you take the negative divergence of it (we're in 1d, so the divergence is again just $\frac{d}{dx}$), you end up with a function $f$ that is not, in fact, a function but a distribution: It has delta functions in all of those places where $\alpha \nabla u$ is discontinuous. One can, in theory, solve the PDE you are interested in with such right hand sides $f$, but not easily with finite elements and for sure not if you use quadrature, because the way you integrate the right hand side terms doesn't see these delta functions.

Secondly, the Poisson equation models some physical behavior. In general, this would be a diffusive process of some sort where $j=-\alpha\nabla u$ is a flux and the divergence of the flux equals the source terms $f$. Fluxes are generally thought to be continuous functions, though not necessarily smooth: they can have kinks, and so their derivatives can be discontinuous. That means that if $\alpha$ is discontinuous, then $\nabla u$ must be discontinuous as well, in such a way that $-\alpha\nabla u$ is continuous again. What this implies is that you can't just choose $\alpha$ and $u$ independently. At its heart, that is why @cfdlab in their answer above starts by prescribing $\alpha$ and the flux (rather than $\alpha$ and $u$) and then deriving what $u$ needs to be.

Answer 2

If you want to use the coefficient of the form $$ \alpha = \begin{cases} \alpha_1 & r < 0.5 \\ \alpha_2 & r > 0.5 \end{cases} $$ it is useful to work in polar coordinates. You build two pieces of the solution $$ u_1(r,\theta), \qquad r < 0.5 $$ and $$ u_2(r,\theta), \qquad r > 0.5 $$ Then at the interface you ensure solution and flux continuity $$ u_1(0.5, \theta) = u_2(0.5, \theta), \qquad \alpha_1 \partial_r u_1(0.5,\theta) = \alpha_2 \partial_r u_2(0.5,\theta) $$ Maybe you can try to build a solution of the form $$ u_i(r,\theta) = f_i(r) g(\theta), \qquad i=1,2 $$ You can first try to build a purely radial solution, which should be easy to do, it will be just like a 1d problem.

Note that the discontinuity is on a circle and to observe good convergence you must make a mesh that aligns with this circle.

If you know fenics, you can see below example where the mesh is made in gmsh https://github.com/cpraveen/fembook/tree/master/fenics/step-10