Solution of the Blatter-Pattyn model

The final form of the equations we'd like to solve is:

\begin{align} & x:\quad 4\frac{\partial }{\partial x}\left( \eta \frac{\partial u}{\partial x} \right)+\frac{\partial }{\partial y}\left( \eta \frac{\partial u}{\partial y} \right)+\frac{\partial }{\partial z}\left( \eta \frac{\partial u}{\partial z} \right)=-2\frac{\partial }{\partial x}\left( \eta \frac{\partial v}{\partial y} \right)-\frac{\partial }{\partial y}\left( \eta \frac{\partial v}{\partial x} \right)+\rho g\frac{\partial s}{\partial x} \\ & y:\quad 4\frac{\partial }{\partial y}\left( \eta \frac{\partial v}{\partial y} \right)+\frac{\partial }{\partial x}\left( \eta \frac{\partial v}{\partial x} \right)+\frac{\partial }{\partial z}\left( \eta \frac{\partial v}{\partial z} \right)=-2\frac{\partial }{\partial y}\left( \eta \frac{\partial u}{\partial x} \right)-\frac{\partial }{\partial x}\left( \eta \frac{\partial u}{\partial y} \right)+\rho g\frac{\partial s}{\partial y} \\ \end{align}

Again, note that for the x equation we've moved all the terms containing gradients in v to the right-hand side (RHS).

We've set it up this way in order to solve the equations using an operator splitting approach; for the x equation, we treat v as known (where we take the values of v from the previous iteration) and solve for u, and vice versa when we solve they y equation for v. The "splitting" refers to the fact that we are breaking the multi-dimensional divergence operation into multiple steps. Rather than solving one big matrix equation for u and v simultaneously we solve two smaller matrix equations in sequence with one of the unknowns treated as a known "source" term.

... unsure on how much more detail to go into on this ...

- symmetry in the equations means that the same subroutines can be used for the discretization of both the u and v equations. For example,