Getting mesh movement working for SWE #92
Description
Here I'm writing a plan for everything that needs to happen in SWE for mesh movement.
For now we assume that the new mesh coordinate is given to us via a field, X1, and that the old mesh
coordinate is X0. This is known at the beginning of the timestep.
v stores the vector field mapping from X0 to X1. If we are on the sphere, then this needs to satisfy
X1 = exp(dt*v)X0.
Perhaps for now we should set v, and use it to determine X1 (after all, in the Monge-Ampere framework, we solve for phi where dt*v = grad phi).
We'll also need a v1 such that
X0 = exp(-dt*v1)X1. This is not equal to -v, because the vector must be rotated into the tangent plane at the new base point.
Here we'll define the state variables U even though they are referred to as x in the code. Might be time to rename all of these U given that we are talking about coordinate a lot.
In the nonlinear loop we iteratively solve for Unp1, which is the fields defined on the new mesh.
For each nonlinear loop:
(1) Advect all the quantities according to the following splitting method:
(a) On mesh X0, solve q_t + L_{u_n - v}q = 0 with q=q_n at time t=t_n, (L is the relevant transport operator for quantity q) until t_n+dt/2 for all transported quantities q.
(b) If we are solving continuity or circulation 1-form advection (the equation we solve for u), then compute q' on the new mesh by solving
<q', gamma'> = <q, gamma>
where gamma' is the test function obtained by composing gamma with the transformation from X0 to X1.
(c) On mesh X1, solve q_t + L_{u_{n+1} - v1}q = 0 with q=q' at time t=t_n + dt/2, until t_{n+1} for all transported quantities q.
(2) compute the forcing term on new mesh X1
(3) solve the linear system formulated on the new mesh and update Unp1
I think that the easiest way to manage this is to have three meshes in state: mesh, old_mesh, new_mesh. At the beginning of the timestep, we assign old_mesh.assign(new_mesh), and compute new_mesh using v. Then, inside the nonlinear loop, we assign mesh.assign(old_mesh) until stage 1b, then assign mesh.assign(new_mesh) until the end of the step.
Finally we update X0.assign(X1) and continue.
Due to all of these linear systems changing meshes, we need to make sure constant_Jacobian=False in all solvers.