The effect of decoupling a coupled system of PDEs

Question

I asked a somewhat similar question previously but perhaps it might have been too specific for anyone to really answer. Here is a bit more general of a question that I am struggling with. Consider the following system:

$$ -\nabla\cdot(D_{1}(u_{2})\nabla u_{1}) = \nabla\cdot\mathbf{f}_{1}(u_{2}) $$ $$ \frac{\partial u_{2}}{\partial t} + \nabla\cdot\mathbf{f}_{2}(u_{1},u_{2}) - \nabla\cdot(D_{2}(u_{2})\nabla u_{2}) = 0 $$

assuming a general set of BCs: $$ u_{i} = u_{i,D}, \;\;\;\text{on}\;\Gamma_{D} $$ $$ D_{i}\nabla u_{i}\cdot\mathbf{n} = u_{i,N}, \;\;\;\text{on} \;\Gamma_{N} $$

Using DGFEM for spatial discretizations and backwards Euler for the time derivative. We decouple like this:

1. Solve for $u_{1}^{k+1}$ using $u_{2}^{k}$: $$ -\nabla\cdot(D_{1}(u_{2}^{k})\nabla u_{1}^{k+1}) = \nabla\cdot\mathbf{f}_{1}(u_{2}^{k}) $$

2. Solve for $u_{2}^{k+1}$ using $u_{1}^{k+1}$: $$ \frac{\partial u_{2}}{\partial t} + \nabla\cdot\mathbf{f}_{2}(u_{1}^{k+1},u_{2}^{k},u_{2}^{k+1}) - \nabla\cdot(D_{2}(u_{2}^{k+1})\nabla u_{2}^{k+1}) = 0 $$

Newton's method is used to handle the non-linearity.

In the particular example I'm looking at, I have $$ \mathbf{f}_{2}(u_{1},u_{2}) = u_{2}(-\nabla u_{1} + u_{2}\nabla u_{2}) $$

so

$$ \mathbf{f}_{2}(u_{1}^{k+1},u_{2}^{k},u_{2}^{k+1}) = u_{2}^{k+1}(-\nabla u_{1}^{k+1} + u_{2}^{k}\nabla u_{2}^{k}) $$

The issue is that the evaluation of $\mathbf{f}_{2}$ at $u_{1}^{k+1}$ and $u_{2}^{k}$ is really messing up my solution of $u_{2}^{k+1}$ to the point where Newton's method is not converging quadratically and I'm not observing the correct convergence of the discretization errors. I originally wanted to evaluate $\mathbf{f}_{2}$ at $u_{2}^{k}$ so that the vector part, $-\nabla u_{1}^{k+1} + u_{2}^{k}\nabla u_{2}^{k}$ could be treated as a constant at each Newton iteration. Otherwise, it would have to be added to the second order term, which I wanted to avoid doing.

To test my code, I'm using method of manufactured solutions. As a test, I'm considering the second equation by itself and evaluating $\mathbf{f}_{2} = \mathbf{f}_{2}(u_{1,\text{exact}}^{k+1}, u_{2}^{k+1}, u_{2,\text{exact}}^{k+1}) = u_{2}^{k+1}(-\nabla u_{1,\text{exact}}^{k+1} + u_{2,\text{exact}}^{k+1}\nabla u_{2,\text{exact}}^{k+1})$ using the exact solutions and there's, no problem here. Newton's method converges quadratically, I'm observing correct convergence of discretization error, etc. so I'm pretty sure the error here lies with the decoupling and not the discretization.

Is this suboptimal performance in the first scenario expected given the decoupling? I know the decoupling will limit the time step, but I'm pretty sure I'm taking the time step sufficiently small.

The next step I was going to take was to split up $\mathbf{f}_{2}$ and toss the $u_{2}\nabla u_{2}$ into the second order term. This would have the advantage of making the formulation more implicit since the only explicit part would be the use of $u_{1}^{k+1}$ in the convective term.

I hadn't been able to find much literature on what effect decoupling would have on convergence, so if anyone could point me in the right direction or offer some advice, that would be great.

Answer 1

This is definitely whats known as a splitting scheme and are very popular in nonlinear optics as well as Quantum simulations, an example being the dirac equation but you can also do it on nonlinear PDES'.

I think the issue you're having is that you're solving for $u_1^{k+1}$ first in terms of $u_2^k$, when you should really be solving for $u_1^{k+1}$ in terms of $u_2^{k+1}$, (if you want to use an upwind scheme in (2).

To see this, notice that the first equation is linear in $u_1$, with coefficients and an inhomogeneity only depending on $u_2$. So 'in principle' (it could be brutally tough), you could find you're solution to $u_1$ from the first equation (1), analytically in terms of a Greens function (which also may depend on $u_2$) convoluted with the LHS of one. ie. let $ Lu_1$ be the RHS of (1) and $f=\nabla \cdot f_2(u_2)$ and let $G(x;u_2)$ satisfy $LG=\delta(x-s)$ and the boundary conditions for $u_1$, then

$u_1=\int_{\Omega} G((x-s);u_2)fds$

So you've reduced the coupled set to an integro-differential equation for (2). However, now look what happens if you want to do an upwind scheme for $u_2^{k+1}$, you find that $u_1^{k+1}$ depends on $u_2^{k+1}$.

There is definitely flexibility though in this decoupling or splitting though that could make for a nice scheme in your equations. A nice introduction to splitting in the context you're interested in can be found in Chapter 13 'Splitting and its cousins' of Boyd 'Chebyshev and Fourier Spectral Methods'.

Hope this helps

Answer 2

It's a bit hard to say without knowing precisely what PDE you're solving; that said, a few places to maybe look are

1. For coupled linear systems, Uzawa and augmented Lagrangian methods are examples of general splitting schemes. There are known results concerning convergence of these methods based on the spectra of the blocks and the Schur complement of the coupled block system (see Bramble et al, on solving saddle point problems).
2. There's a wealth of material on different splitting methods for PDEs. Though these are largely formulated for linear problems, there is ample use and theory for their extension to nonlinear equations - for example, at least for the nonlinear PDEs governing types of convection-diffusion and incompressible flow, splitting/projection schemes have a pretty rich history.

One thought is that solutions to the decoupled PDEs may or may not be close to solutions of the full coupled systems. Thus, solving only one equation "exactly" with Newton might not buy you more (or may be worse) than simply taking one single Newton iteration and updating both solutions (i.e. taking a small implicit step, then updating $u_1$ and $u_2$ again, instead of solving Newton for $u_1$ and $u_2$ separately at each implicit step). But then again, I'm not sure about this.