Iterative Block Matrix Splitting for Multiphysics Simulation

Question

I have a problem of the form

$$\left[\begin{array}{cc} -(\lambda+2\mu)\frac{d^2}{dx^2} & \alpha\frac{d}{dx} \\ \frac{\alpha}{\Delta t}\frac{d}{dx} & \frac{c_0}{\Delta t}I-\frac{\kappa}{\mu_f}\frac{d^2}{dx^2} \end{array}\right] \left[ \begin{array}{c} u\\ p\end{array}\right] =\left[\begin{array}{c} f\\g \end{array}\right]$$

When I apply a finite difference scheme to this system, I obtain a block matrix system of the form

$$\left[ \begin{array}{cc} m_{11} & m_{12}\\ m_{21} & m_{22}\end{array}\right]\left[ \begin{array}{c} \vec{u}\\ \vec{p}\end{array}\right] =\left[\begin{array}{c} f\\g \end{array}\right]$$

where $m_{ij}$ are all sparse symmetric matrices and $u_i$, $p_i$ are the discrete values of u & p at the discrete points $x_i$ in the domain.

I was thinking of solving this system of equations by a fixed point iteration by splitting the block matrices into the form

$$\left[ \begin{array}{cc} m_{11} & m_{12}\\ 0 & m_{22}\end{array}\right]\left[ \begin{array}{c} \vec{u}\\ \vec{p}\end{array}\right]^{k+1} +\left[ \begin{array}{cc} 0 & 0\\ m_{21} & 0\end{array}\right] \left[ \begin{array}{c} \vec{u}\\ \vec{p}\end{array}\right]^k =\left[\begin{array}{c} f\\g \end{array}\right]$$

I recall from basic numerical analysis that if $A=M+N$ and we iterate $Mx^{k+1}=-Nx^k+b$, then the method is guaranteed to converge from any starting vector $x^0$ if the eigenvalues of $M^{-1}N$ are $<1$ in absolute value. Is there an easy way to verify whether this is true for this particular matrix splitting? Can I apply an analogous proof as that of standard gauss-seidel or jacobi methods?

Note: I know that this system can be solved by other (possibly faster) means in the 1D case, but I'm particularly interested in this form because it allows me to decouple the multiphysics simulation and I can implement fast solvers for each equation in the more complicated multidimensional case.

Answer 1

When you write down the product

$$M^{-1}N = \begin{pmatrix} m_{11}^{-1} & -m_{11}^{-1} m_{12} m_{22}^{-1} \\ 0 & m_{22}^{-1} \end{pmatrix} \begin{pmatrix} 0 & 0 \\ m_{21} & 0 \end{pmatrix} = \begin{pmatrix} -m_{11}^{-1} m_{12} m_{22}^{-1} m_{21} & 0 \\ m_{22}^{-1} m_{21} & 0 \end{pmatrix}$$

it is clear that the eigenvalues of $M^{-1} N$ depend on $\Delta t$. Note that you can consider the generalized eigenvalue problem $m_{12} m_{22}^{-1} m_{21} x = \lambda m_{11} x$. The stability of this iteration depends on $\Delta t$ and the relative sizes of the two terms in $m_{22}$.

If you want a robust solver based on splitting, you should look at approximate commutator preconditioners. The idea there is to approximate the Schur complement $m_{11} - m_{12} m_{22}^{-1} m_{21}$ or $m_{22} - m_{21} m_{11}^{-1} m_{12}$ using something that is feasible to apply. These are usually combined with a Krylov method to catch stray eigenvalues arising because the approximate commutator arguments are inexact, especially at boundaries.