Solution of the Blatter-Pattyn model

From Interactive System for Ice sheet Simulation
Revision as of 01:01, 12 March 2011 by Sprice (Talk | contribs)

Jump to: navigation, search


Governing Equations

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

  & 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} \\ 

Coordinate Transform

As with the 0-order model, we need to change from Cartesian to sigma coordinates. The first normal-stress term first term on the left-hand side becomes

Pattyn introduces a dimensionless vertical coordinate so that variations in ice thickness do not complicate the numerical treatment considerably (Pattyn 2003). This dimensionless vertical coordinate is referred to as \zeta by Pattyn. I will refer to it as \sigma for consistency with CISM's existing notation. The rescaled coordinate is defined as:

\sigma = \frac{(s - z)}{H}

This means that at the surface of the ice sheet \sigma = 0, and at the base \sigma = 1 regardless of the ice thickness (Pattyn 2003). As a result of this transformation, a coordinate (x,y,z) is mapped to (x',y',\sigma) (Pattyn 2003). This means that function derivatives must be re-written (using \frac{\partial f}{\partial x} as an example) as:

\frac{\partial f}{\partial x} = \frac{\partial f}{\partial x'} \frac{\partial x'}{\partial x} + \frac{\partial f}{\partial y'} \frac{\partial y'}{\partial x} + \frac{\partial f}{\partial \sigma} \frac{\partial \sigma}{\partial x}

Similarly for \frac{\partial f}{\partial y} and \frac{\partial f}{\partial z}. Pattyn simplifies this by assuming that

\frac{\partial x'}{\partial x}, \frac{\partial y'}{\partial y} = 1


\frac{\partial x'}{\partial y}, \frac{\partial x'}{\partial z}, \frac{\partial y'}{\partial x}, \frac{\partial y'}{\partial z} = 0.

This assumption is valid if the bed and surface gradients are not too large (Pattyn 2003). This simplifies the above to:

\frac{\partial f}{\partial x} = \frac{\partial f}{\partial x'} + \frac{\partial f}{\partial \sigma}\frac{\partial \sigma}{\partial x}

\frac{\partial f}{\partial y} = \frac{\partial f}{\partial y'} + \frac{\partial f}{\partial \sigma}\frac{\partial \sigma}{\partial y}

\frac{\partial f}{\partial z} = \frac{\partial f}{\partial \sigma}\frac{\partial \sigma}{\partial z}

Rescaling parameters a_x, a_y, b_x, b_y, and c_{xy} are defined. Presenting the x derivative case, as the y derivative case is analogous,

a_x = \frac{1}{H}(\frac{\partial s}{\partial x'} - \sigma \frac{\partial H}{\partial x'})

b_x = \frac{\partial a_x}{\partial x'} + a_x \frac{\partial a_x}{\partial \sigma} 
    = \frac{1}{H} (\frac{\partial^2 s}{\partial x'^2} - \sigma \frac{\partial^2 H}{\partial x'^2} - 2a_x \frac{\partial H}{\partial x'})

c_{xy} = \frac{\partial a_y}{\partial x'} + a_x \frac{\partial a_y}{\partial \sigma} 
       = \frac{\partial a_x}{\partial y'} + a_y \frac{\partial a_x}{\partial \sigma}

Using these, expressions for the x derivatives become:

\frac{\partial f}{\partial x} = \frac{\partial f}{\partial x'} + a_x \frac{\partial f}{\partial \sigma}

\frac{\partial }{\partial x}\left( \eta \frac{\partial u}{\partial x} \right)=\frac{\partial }{\partial \hat{x}}\left( \eta \frac{\partial u}{\partial \hat{x}} \right)+\frac{\partial \sigma }{\partial \hat{x}}\frac{\partial }{\partial \sigma }\left( \eta \frac{\partial u}{\partial \hat{x}} \right)+\frac{\partial \sigma }{\partial \hat{x}}\frac{\partial }{\partial \hat{x}}\left( \eta \frac{\partial u}{\partial \sigma } \right)+\left( \frac{\partial \sigma }{\partial \hat{x}} \right)^{2}\frac{\partial }{\partial \sigma }\left( \eta \frac{\partial u}{\partial \sigma } \right)+\left( \frac{\partial _{{}}^{2}\sigma }{\partial \hat{x}_{{}}^{2}} \right)\eta \frac{\partial u}{\partial \sigma }

where hatted values refer to the coordinate directions in sigma coordinates. Similarly, the first cross-stress term on the RHS is given by

\frac{\partial }{\partial x}\left( \eta \frac{\partial v}{\partial y} \right)=\underset{{}}{\mathop{\frac{\partial }{\partial \hat{x}}\left( \eta \frac{\partial v}{\partial \hat{y}} \right)}}\,+\underset{{}}{\mathop \frac{\partial \sigma }{\partial \hat{x}}\frac{\partial }{\partial \sigma }\left( \eta \frac{\partial v}{\partial \hat{y}} \right)}\,+\underset{{}}{\mathop \frac{\partial \sigma }{\partial \hat{y}}\frac{\partial }{\partial \hat{x}}\left( \eta \frac{\partial v}{\partial \sigma } \right)}\,+\underset{{}}{\mathop \frac{\partial \sigma }{\partial \hat{x}}\frac{\partial \sigma }{\partial \hat{y}}\frac{\partial }{\partial \sigma }\left( \eta \frac{\partial v}{\partial \sigma } \right)}\,+\underset{{}}{\mathop \frac{\partial _{{}}^{2}\sigma }{\partial \hat{x}\partial \hat{y}}\eta \frac{\partial v}{\partial \sigma }}\,

One term has become five terms and each one of those is pretty scary looking on its own. Luckily, there is also a lot of symmetry here. Notice that if we wanted to design subroutines to discretize the terms on the RHS, we could re-use a lot of them by either applying them to the correct velocity component (either to the u OR the v discretization) or by passing the appropriate arguments (pass either the grid spacing in the x direction OR the y direction, where appropriate). There is still a lot of work here ...

Operating Splitting

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.

General Matrix Form

A general matrix form of the equations, where coefficients on the u and v velocity components (i.e. viscosity, grid spacing, scalars) are contained in the block matrices A, is given by

\left[ \begin{matrix}
   \mathbf{A}_{\mathbf{uu}} & \mathbf{A}_{\mathbf{uv}}  \\
   \mathbf{A}_{\mathbf{vu}} & \mathbf{A}_{\mathbf{vv}}  \\
\end{matrix} \right]\left[ \begin{matrix}
   \mathbf{u}  \\
   \mathbf{v}  \\
\end{matrix} \right]=\left[ \begin{matrix}
   \mathbf{b}_{\mathbf{u}}  \\
   \mathbf{b}_{\mathbf{v}}  \\
\end{matrix} \right]

where the uu subscript denotes block matrices containing coefficients for gradients on u in the equation for the x component of velocity (i.e. u). The subscript uv denotes block matrices containing coefficients for gradients on v in the equation for the x component of velocity (and similarly for the vv and vu subscripts). On the right-hand side, the subscripts denote the geometric source terms for the x and y components of velocity (subscripts u and v, respectively).

Solution of the Non-linear System Through a Fixed Point Iteration

The non-linearity in the equations - the fact that the coefficients on the velocity components (the viscosity) are dependent on the velocity (or more specifically, the velocity gradients) - is handled through a fixed-point iteration.

Go to Strang for this ..

Final Matrix Form

The final matrix form of the equations, accounting for the Picard iteration on the viscosity, is given by

  \left[ \begin{matrix}
   \mathbf{A}_{\mathbf{uu}}^{k-1} & \mathbf{0}  \\
   \mathbf{0} & \mathbf{A}_{\mathbf{vv}}^{k-1}  \\
\end{matrix} \right]\left[ \begin{matrix}
   \mathbf{u}^{k}  \\
   \mathbf{v}^{k}  \\
\end{matrix} \right]=\left[ \begin{matrix}
   \mathbf{b}_{\mathbf{u}}-\mathbf{A}_{\mathbf{uv}}^{k-1}\mathbf{v}^{k-1}  \\
   \mathbf{b}_{\mathbf{v}}-\mathbf{A}_{\mathbf{vu}}^{k-1}\mathbf{u}^{k-1}  \\
\end{matrix} \right] \\ 
  \mathbf{A}_{\mathbf{uu}}^{k-1}\mathbf{u}=\mathbf{b}_{\mathbf{u}}-\mathbf{A}_{\mathbf{uv}}^{k-1}\mathbf{v}^{k-1},\quad \quad \mathbf{A}_{\mathbf{vv}}^{k-1}\mathbf{v}=\mathbf{b}_{\mathbf{v}}-\mathbf{A}_{\mathbf{vu}}^{k-1}\mathbf{u}^{k-1} \\ 

where the index k denotes an unknown value being solved for during the current non-linear iteration and the index k-1 denotes a lagged value taken from solution at the end of the previous non-linear iteration.

Solution of the Linear System

Newton-based Methods for Solutions of the Non-linear System