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

Line 11: | Line 11: | ||

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. | 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. | ||

+ | |||

+ | ===Sigma Coordinate=== | ||

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 | 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 <math>\zeta</math> by Pattyn. I will refer to it as <math>\sigma</math> for consistency with CISM's existing notation. The rescaled coordinate is defined as: | ||

+ | |||

+ | <math>\sigma = \frac{(s - z)}{H}</math> | ||

+ | |||

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

+ | |||

+ | <math> | ||

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

+ | </math> | ||

+ | |||

+ | Similarly for <math>\frac{\partial f}{\partial y}</math> and <math>\frac{\partial f}{\partial z}</math>. Pattyn simplifies this by assuming that | ||

+ | |||

+ | <math>\frac{\partial x'}{\partial x}, \frac{\partial y'}{\partial y} = 1</math> | ||

+ | |||

+ | and | ||

+ | |||

+ | <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: | ||

+ | |||

+ | <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 y} = \frac{\partial f}{\partial y'} + \frac{\partial f}{\partial \sigma}\frac{\partial \sigma}{\partial y}</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, | ||

+ | |||

+ | <math>a_x = \frac{1}{H}(\frac{\partial s}{\partial x'} - \sigma \frac{\partial H}{\partial x'})</math> | ||

+ | |||

+ | <math> | ||

+ | 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'}) | ||

+ | </math> | ||

+ | |||

+ | <math> | ||

+ | 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} | ||

+ | </math> | ||

+ | |||

+ | Using these, expressions for the x derivatives become: | ||

+ | |||

+ | <math> | ||

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

+ | </math> | ||

## Revision as of 08:54, 13 August 2009

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

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.

### Sigma Coordinate

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 by Pattyn. I will refer to it as for consistency with CISM's existing notation. The rescaled coordinate is defined as:

This means that at the surface of the ice sheet , and at the base regardless of the ice thickness (Pattyn 2003). As a result of this transformation, a coordinate is mapped to (Pattyn 2003). This means that function derivatives must be re-written (using as an example) as:

Similarly for and . Pattyn simplifies this by assuming that

and

.

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

Rescaling parameters , , , , and are defined. Presenting the x derivative case, as the y derivative case is analogous,

Using these, expressions for the x derivatives become:

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

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 ...