Difference between revisions of "Solution of the Blatter-Pattyn model"

From Interactive System for Ice sheet Simulation
Jump to: navigation, search
(Coordinate Transform)
Line 20: Line 20:
 
</math>
 
</math>
  
Similarly for <math>\frac{\partial f}{\partial y}</math> and <math>\frac{\partial f}{\partial z}</math>. Pattyn simplifies this by assuming that  
+
Similarly for <math>\frac{\partial f}{\partial y}</math> and <math>\frac{\partial f}{\partial z}</math>. We can simplify this by assuming that  
  
 
<math>\frac{\partial x'}{\partial x}, \frac{\partial y'}{\partial y} = 1</math>  
 
<math>\frac{\partial x'}{\partial x}, \frac{\partial y'}{\partial y} = 1</math>  
Line 28: Line 28:
 
<math>\frac{\partial x'}{\partial y}, \frac{\partial x'}{\partial z}, \frac{\partial y'}{\partial x}, \frac{\partial y'}{\partial z} = 0</math>.   
 
<math>\frac{\partial x'}{\partial y}, \frac{\partial x'}{\partial z}, \frac{\partial y'}{\partial x}, \frac{\partial y'}{\partial z} = 0</math>.   
  
This assumption is valid if the bed and surface gradients are not too large (Pattyn 2003). This simplifies the above to:
+
This assumption is valid if the bed and surface gradients are not too large. This simplifies the above to:
  
 
<math>\frac{\partial f}{\partial x} = \frac{\partial f}{\partial x'} + \frac{\partial f}{\partial \sigma}\frac{\partial \sigma}{\partial x}</math>
 
<math>\frac{\partial f}{\partial x} = \frac{\partial f}{\partial x'} + \frac{\partial f}{\partial \sigma}\frac{\partial \sigma}{\partial x}</math>
Line 36: Line 36:
 
<math>\frac{\partial f}{\partial z} = \frac{\partial f}{\partial \sigma}\frac{\partial \sigma}{\partial z}</math>
 
<math>\frac{\partial f}{\partial z} = \frac{\partial f}{\partial \sigma}\frac{\partial \sigma}{\partial z}</math>
  
Rescaling parameters <math>a_x</math>, <math>a_y</math>, <math>b_x</math>, <math>b_y</math>, and <math>c_{xy}</math> are defined. Presenting the x derivative case, as the y derivative case is analogous,
+
Rescaling parameters <big><math>a_x</math></big>, <big><math>a_y</math></big>, <big><math>b_x</math></big>, <big><math>b_y</math></big>, and <big><math>c_{xy}</math></big> are defined. For the ''x'' derivative case (the ''y'' derivative case is analogous) we have
  
 
<math>a_x = \frac{1}{H}(\frac{\partial s}{\partial x'} - \sigma \frac{\partial H}{\partial x'})</math>
 
<math>a_x = \frac{1}{H}(\frac{\partial s}{\partial x'} - \sigma \frac{\partial H}{\partial x'})</math>
Line 50: Line 50:
 
</math>
 
</math>
  
Using these, expressions for the x derivatives become:
+
Using these, expressions for the ''x'' derivatives become:
  
 
<math>
 
<math>
Line 66: Line 66:
  
  
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 ...
+
One term has become five terms and each one of those is pretty ugly looking on its own. Luckily, there is 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 (to either the ''U'' or the ''V'' discretization) or by passing the appropriate arguments (by passing either the grid spacing in the ''X'' direction or the ''Y'' direction, where appropriate).
  
 
===Operating Splitting===
 
===Operating Splitting===

Revision as of 11:37, 14 March 2011

Contents

Governing Equations

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}

Coordinate Transform

For ice sheet modeling, it is convenient to recast the governing equations using a dimensionless, stretched vertical coordinate (often called a sigma coordinates). The stretched vertical 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. As a result of this transformation, a coordinate (x,y,z) is mapped to (x',y',\sigma). 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}. We can simplify this by assuming that

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

and

\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. 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. For the x derivative case (the y derivative case is analogous) we have

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 ugly looking on its own. Luckily, there is 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 (to either the U or the V discretization) or by passing the appropriate arguments (by passing either the grid spacing in the X direction or the Y direction, where appropriate).

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


\begin{matrix}
  \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} \\ 
\end{matrix}


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

Go to Blatter-Pattyn model.

Go to Blatter-Pattyn Boundary Conditions.

Go to CISM governing equations and numerical solution.

Return to main coarse page