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

(→The Jacobian Free Approach) |
|||

(62 intermediate revisions by 2 users not shown) | |||

Line 1: | Line 1: | ||

+ | ===Governing Equations=== | ||

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

Line 7: | Line 8: | ||

\end{align}</math> | \end{align}</math> | ||

+ | ===Coordinate Transform=== | ||

− | + | For ice sheet modeling, it is convenient to recast the governing equations using a dimensionless, stretched vertical coordinate (often called a [http://atmo.tamu.edu/class/metr452/models/2001/vertres.html ''sigma coordinates'']). The stretched vertical coordinate is defined as: | |

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

− | + | This means that at the surface of the ice sheet <big><math>\sigma = 0</math></big>, and at the base <big><math>\sigma = 1</math></big> regardless of the ice thickness. As a result of this transformation, a coordinate <big><math>(x,y,z)</math></big> is mapped to <big><math>(x',y',\sigma)</math></big>. 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>. We can simplify 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. 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 <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> | ||

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

Line 24: | Line 66: | ||

− | + | 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). | |

+ | |||

+ | A similar transform is applied to each of the terms in the governing equations given above. At any point within the grid, the grid spacing, coordinate transform, and viscosity information associated with the unknown velocity components (''U'' and ''V'') is made discrete using [[Finite differencing | finite differences]]. This information ultimately equates to coefficients on the unknown velocities, allowing the governing equations over the entire grid (with appropriate discretizations for boundary conditions) to be recast as a system of ''n'' equations in ''n'' unknowns. In turn, this system is solved using standard linear algebraic methods for large, [[http://en.wikipedia.org/wiki/Sparse_matrix sparse]] systems of equations. | ||

+ | |||

+ | ===Operating Splitting=== | ||

+ | |||

+ | In the governing equations given above, note that for the ''x'' equation we have moved all terms containing gradients in ''v'' to the right-hand side (RHS) (and vice-versa for the ''y'' equation). | ||

+ | |||

+ | This allows us 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, as discussed further below) and solve for ''u'', and vice versa when we solve the ''y'' equation for ''v''. The "splitting" refers to the fact that we are breaking the multi-dimensional divergence operation into two 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. This procedure was probably more common and important years ago when it was desirable to keep the matrix equations as small as possible for memory management issues. On today's machines, with fewer memory limitations (in particular when dealing with codes designed to run on parallel, distributed memory architectures) this splitting is not necessary and may even lead to some undesirable numerical side effects (i.e. a slow-down in the convergence of iterations used to treat nonlinearity in the governing equations). | ||

+ | |||

+ | A general matrix form of the "split" 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 | ||

+ | |||

+ | |||

+ | <math>\begin{matrix} | ||

+ | \left[ \begin{matrix} | ||

+ | \mathbf{A}_{\mathbf{uu}} & \mathbf{0} \\ | ||

+ | \mathbf{0} & \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{A}_{\mathbf{uv}}\mathbf{v} \\ | ||

+ | \mathbf{b}_{\mathbf{v}}-\mathbf{A}_{\mathbf{vu}}\mathbf{u} \\ | ||

+ | \end{matrix} \right] \\ | ||

+ | \\ | ||

+ | \mathbf{A}_{\mathbf{uu}}\mathbf{u}=\mathbf{b}_{\mathbf{u}}-\mathbf{A}_{\mathbf{uv}}\mathbf{v},\quad \quad \mathbf{A}_{\mathbf{vv}}\mathbf{v}=\mathbf{b}_{\mathbf{v}}-\mathbf{A}_{\mathbf{vu}}\mathbf{u} \\ | ||

+ | \end{matrix}</math> | ||

+ | |||

+ | |||

+ | 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 single subscripts '''u''' and '''v''' are attached to the geometric source terms for the ''x'' and ''y'' components of velocity, 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 [http://en.wikipedia.org/wiki/Fixed_point_iteration "fixed-point iteration"]. A general fixed point iteration for a vector of unknowns ''u'' can be written as | ||

+ | |||

+ | |||

+ | |||

+ | <math>u^{k}=\mathbf{B}\left( u^{k-1} \right)</math>, | ||

+ | |||

+ | |||

+ | where ''k'' is the index for the ''u'' being solved for and '''B''' is a matrix operation performed on the components of ''u'' obtained at the previous iteration, ''k-1''. The "fixed point" occurs when the values of ''u'' at ''k'' and ''k-1'' are equal to within some given tolerance (at which point the iteration process is halted). CISM has options for implementing both "Picard" and "Newton"-based fixed-point iterations. For the Picard iteration (standard in CISM), the matrix coefficients with a velocity dependence are simply based on the velocities obtained at the previous iteration. In most cases, this equates to using velocities obtained from the previous iteration to calculate the strain rate components that go into the calculation of the effective viscosity, <big><math>\eta</math></big>. | ||

+ | |||

+ | ===Final Matrix Form=== | ||

+ | |||

+ | When accounting for both the operator splitting and the Picard iteration on the effective viscosity, the final form of the matrix equations solved in CISM becomes | ||

+ | |||

+ | |||

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

+ | |||

+ | |||

+ | 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 (again, here the lagging is primarily with respect to the effective viscosity, the value of which is calculated using velocity gradients obtained at the end of the previous iteration). The final form of the matrix equations given above represents a linear system; for the solution at any particular nonlinear iteration ''k'', all of the coefficients on the unknown velocity components ''u'' and ''v'' are held "frozen" during the solution of the linear system. This linear system can be solved using any practical method. For large, sparse systems, some variant on the iterative conjugate gradient method (e.g. BiCG, GMRES) is generally the most efficient. In this case the linear system is not solved exactly but is solved to within some small tolerance of the "true" solution. | ||

+ | |||

+ | ===Newton-based Methods for Solutions of the Non-linear System=== | ||

+ | |||

+ | Without any operator splitting, the generic matrix form of the equations to be solved can be written as | ||

+ | |||

+ | |||

+ | <math>\mathbf{A}(\mathbf{u})\mathbf{u}=\mathbf{b}</math>. | ||

+ | |||

+ | |||

+ | The linearized form of the equations to be solved using the Picard solution can be written as | ||

+ | |||

+ | |||

+ | <math>\mathbf{u}^{k}=\mathbf{A}(\mathbf{u}^{k-1})^{-1}\mathbf{b}</math>. | ||

+ | |||

+ | |||

+ | The full nonlinear system to be solved can be written as | ||

+ | |||

+ | |||

+ | <math>\mathbf{F}(\mathbf{u})=\mathbf{A}(\mathbf{u})\mathbf{u}-\mathbf{b}</math> | ||

+ | |||

+ | |||

+ | with the solution for the uknown vector ''u'' given by | ||

+ | |||

+ | |||

+ | <math>\mathbf{F}(\mathbf{u})=0</math>. | ||

+ | |||

+ | |||

+ | |||

+ | A Newton-based solution for this system of equations, based on a first-order Taylor series expansion about the solution for ''u'' at iteration ''k-1'', can be written as | ||

+ | |||

+ | |||

+ | <math>\mathbf{F}(\mathbf{u}^{k})=\mathbf{F}(\mathbf{u}^{k-1})+\mathbf{{F}'}(\mathbf{u}^{k-1})\delta \mathbf{u}^{k-1}</math>, | ||

+ | |||

+ | |||

+ | where | ||

+ | |||

+ | |||

+ | <math>\mathbf{{F}'}(\mathbf{u}^{k-1})=\mathbf{J}^{k-1}</math> | ||

+ | |||

+ | |||

+ | is the system Jacobian with individual components give by | ||

+ | |||

+ | |||

+ | |||

+ | <math>J_{ij}=\frac{\partial F_{i}(\mathbf{u}^{k-1})}{\partial u_{j}}</math> | ||

+ | |||

+ | |||

+ | and | ||

+ | |||

+ | |||

+ | <math>\delta \mathbf{u}^{k-1}=\mathbf{u}^{k}-\mathbf{u}^{k-1}</math> | ||

+ | |||

+ | |||

+ | is the Newton update to be solved for. One method for doing so is by solving | ||

+ | |||

+ | |||

+ | <math>\delta \mathbf{u}^{k-1}=-\left( \mathbf{J}^{k-1} \right)^{-1}\mathbf{F}(\mathbf{u}^{k-1})</math>. | ||

+ | |||

+ | |||

+ | The advantage of Newton-based methods is that, with a good initial guess for the solution, convergence rates are very often quadratic (e.g. the residual decreases quadratically, so that at iteration ''k'' one has a residual of 0.1, at iteration ''k+1'', a residual of 0.01, and at iteration ''k+2'' a residual of 0.0001), whereas Picard-based iterations are much slower to converge. A figure comparing rates of convergence for Picard versus Newton on a CISM tests case is shown below. The Newton method is based on the work of Lemieux et. al (submitted to ''JCP''), which is discussed further below. | ||

+ | |||

+ | [[Image:picard_vs_newton.jpg|Rates of convergence for the nonlinear iteration in CISM.]] | ||

+ | |||

+ | ===The Jacobian Free Approach=== | ||

+ | |||

+ | In practice, the model Jacobian may either be too difficult or to expensive too form. A "Jacobian Free Newton-Krylov" (JFNK) approach has recently been implemented in CISM (Leimieux et al., submitted to ''JCP''), largely following methods discussed in [http://www.sciencedirect.com/science?_ob=ArticleURL&_udi=B6WHY-49KSNJ6-4&_user=2493154&_coverDate=01%2F20%2F2004&_alid=1678413922&_rdoc=1&_fmt=high&_orig=search&_origin=search&_zone=rslt_list_item&_cdi=6863&_sort=r&_st=13&_docanchor=&view=c&_ct=1&_acct=C000057551&_version=1&_urlVersion=0&_userid=2493154&md5=887f929a8adbac698ecf277ac0c99cef&searchtype=a Knoll and Keyes (2004)]. The crux of the method comes from noting that, when solving the last equation above using a Krylov method (e.g. Conjugate Gradients, GMRES, etc.) the solution for the Newton update is taken from a combination of Krylov vectors that span the subspace | ||

+ | |||

+ | |||

+ | <math>span\left\{ \mathbf{r}_{0},\mathbf{Jr}_{0},\mathbf{J}^{2}\mathbf{r}_{0},...,\mathbf{J}^{n-1}\mathbf{r}_{0} \right\}=span\left\{ \mathbf{r}_{0},\mathbf{Jv}_{1},\mathbf{Jv}_{2},...,\mathbf{Jv}_{n-1} \right\}</math>. | ||

+ | |||

+ | |||

+ | This implies that, when using a Krylov method, one only ever needs to calculate matrix vector products of the form <big><math>\mathbf{Jv}</math></big> when building up the subspace that approximates the solution vector <math>\delta \mathbf{u}</math>. | ||

+ | |||

+ | |||

+ | Following Knoll and Keyes (2004), note that the necessary matrix vector products can be approximated through nonlinear function evaluations and a perturbation as | ||

+ | |||

+ | |||

+ | <math>\mathbf{Jv}\approx \frac{\mathbf{F}\left( \mathbf{u}+\varepsilon \mathbf{v} \right)-\mathbf{F}\left( \mathbf{u} \right)}{\varepsilon }</math>. | ||

+ | |||

+ | |||

+ | It is not immediately obvious why this approximation is valid. To verify this, take a few steps back and consider a nonlinear system of equations of two variables, ''u1'' and ''u2''. The right-hand side of the above equation can be expanded as | ||

+ | |||

+ | |||

+ | <math>\frac{\mathbf{F}\left( \mathbf{u}+\varepsilon \mathbf{v} \right)-\mathbf{F}\left( \mathbf{u} \right)}{\varepsilon }=\left[ \begin{matrix} | ||

+ | \frac{F_{1}\left( u_{1}+\varepsilon v_{1},u_{2}+\varepsilon v_{2} \right)-F_{1}(u_{1},u_{2})}{\varepsilon } \\ | ||

+ | \frac{F_{2}\left( u_{1}+\varepsilon v_{1},u_{2}+\varepsilon v_{2} \right)-F_{2}(u_{1},u_{2})}{\varepsilon } \\ | ||

+ | \end{matrix} \right]</math>. | ||

+ | |||

+ | |||

+ | A first-order Taylor series expansion approximation to this is given by | ||

+ | |||

+ | |||

+ | <math>\frac{\mathbf{F}\left( \mathbf{u}+\varepsilon \mathbf{v} \right)-\mathbf{F}\left( \mathbf{u} \right)}{\varepsilon }\approx \left[ \begin{matrix} | ||

+ | \frac{F_{1}\left( u_{1},u_{2} \right)+\varepsilon v_{1}\frac{\partial F_{1}}{\partial u_{1}}+\varepsilon v_{2}\frac{\partial F_{1}}{\partial u_{2}}-F_{1}(u_{1},u_{2})}{\varepsilon } \\ | ||

+ | \frac{F_{2}\left( u_{1},u_{2} \right)+\varepsilon v_{1}\frac{\partial F_{2}}{\partial u_{1}}+\varepsilon v_{2}\frac{\partial F_{2}}{\partial u_{2}}-F_{2}(u_{1},u_{2})}{\varepsilon } \\ | ||

+ | \end{matrix} \right]</math>, | ||

+ | |||

+ | |||

+ | which collapses to | ||

+ | |||

+ | |||

+ | <math>\frac{\mathbf{F}\left( \mathbf{u}+\varepsilon \mathbf{v} \right)-\mathbf{F}\left( \mathbf{u} \right)}{\varepsilon }\approx\left[ \begin{matrix} | ||

+ | v_{1}\frac{\partial F_{1}}{\partial u_{1}}+v_{2}\frac{\partial F_{1}}{\partial u_{2}} \\ | ||

+ | v_{1}\frac{\partial F_{2}}{\partial u_{1}}+v_{2}\frac{\partial F_{2}}{\partial u_{2}} \\ | ||

+ | \end{matrix} \right]</math>. | ||

+ | |||

+ | |||

+ | Finally, note that the right-hand side of the above equation is equal to | ||

+ | |||

+ | |||

+ | <math>\mathbf{Jv}\approx\left[ \begin{matrix} | ||

+ | v_{1}\frac{\partial F_{1}}{\partial u_{1}}+v_{2}\frac{\partial F_{1}}{\partial u_{2}} \\ | ||

+ | v_{1}\frac{\partial F_{2}}{\partial u_{1}}+v_{2}\frac{\partial F_{2}}{\partial u_{2}} \\ | ||

+ | \end{matrix} \right]</math>, | ||

+ | |||

+ | |||

+ | with the Jacobian matrix given by | ||

+ | |||

+ | |||

+ | |||

+ | <math>\mathbf{J}=\left[ \begin{matrix} | ||

+ | \frac{\partial F_{1}}{\partial u_{1}} & \frac{\partial F_{1}}{\partial u_{2}} \\ | ||

+ | \frac{\partial F_{2}}{\partial u_{1}} & \frac{\partial F_{2}}{\partial u_{2}} \\ | ||

+ | \end{matrix} \right]</math>. | ||

+ | |||

+ | |||

+ | |||

+ | This matrix vector product is what needs to be calculated repeatedly while building up the Krylov subspace vectors that combine to approximate the Newton update vector <math>\delta \mathbf{u}</math>. The important point is that at no point in this process does one need to calculate the entire Jacobian matrix. Another important point is that the accuracy of the approximation to the Jacobian is proportional to the small perturbation term, <math>\varepsilon</math>. |

## Latest revision as of 16:10, 14 March 2011

## Contents |

### Governing Equations

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

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

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

Similarly for and . We can simplify this by assuming that

and

.

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

Rescaling parameters , , , , and are defined. For the *x* derivative case (the *y* derivative case is analogous) we have

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

A similar transform is applied to each of the terms in the governing equations given above. At any point within the grid, the grid spacing, coordinate transform, and viscosity information associated with the unknown velocity components (*U* and *V*) is made discrete using finite differences. This information ultimately equates to coefficients on the unknown velocities, allowing the governing equations over the entire grid (with appropriate discretizations for boundary conditions) to be recast as a system of *n* equations in *n* unknowns. In turn, this system is solved using standard linear algebraic methods for large, [sparse] systems of equations.

### Operating Splitting

In the governing equations given above, note that for the *x* equation we have moved all terms containing gradients in *v* to the right-hand side (RHS) (and vice-versa for the *y* equation).

This allows us 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, as discussed further below) and solve for *u*, and vice versa when we solve the *y* equation for *v*. The "splitting" refers to the fact that we are breaking the multi-dimensional divergence operation into two 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. This procedure was probably more common and important years ago when it was desirable to keep the matrix equations as small as possible for memory management issues. On today's machines, with fewer memory limitations (in particular when dealing with codes designed to run on parallel, distributed memory architectures) this splitting is not necessary and may even lead to some undesirable numerical side effects (i.e. a slow-down in the convergence of iterations used to treat nonlinearity in the governing equations).

A general matrix form of the "split" 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

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 single subscripts **u** and **v** are attached to the geometric source terms for the *x* and *y* components of velocity, 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". A general fixed point iteration for a vector of unknowns *u* can be written as

,

where *k* is the index for the *u* being solved for and **B** is a matrix operation performed on the components of *u* obtained at the previous iteration, *k-1*. The "fixed point" occurs when the values of *u* at *k* and *k-1* are equal to within some given tolerance (at which point the iteration process is halted). CISM has options for implementing both "Picard" and "Newton"-based fixed-point iterations. For the Picard iteration (standard in CISM), the matrix coefficients with a velocity dependence are simply based on the velocities obtained at the previous iteration. In most cases, this equates to using velocities obtained from the previous iteration to calculate the strain rate components that go into the calculation of the effective viscosity, .

### Final Matrix Form

When accounting for both the operator splitting and the Picard iteration on the effective viscosity, the final form of the matrix equations solved in CISM becomes

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 (again, here the lagging is primarily with respect to the effective viscosity, the value of which is calculated using velocity gradients obtained at the end of the previous iteration). The final form of the matrix equations given above represents a linear system; for the solution at any particular nonlinear iteration *k*, all of the coefficients on the unknown velocity components *u* and *v* are held "frozen" during the solution of the linear system. This linear system can be solved using any practical method. For large, sparse systems, some variant on the iterative conjugate gradient method (e.g. BiCG, GMRES) is generally the most efficient. In this case the linear system is not solved exactly but is solved to within some small tolerance of the "true" solution.

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

Without any operator splitting, the generic matrix form of the equations to be solved can be written as

.

The linearized form of the equations to be solved using the Picard solution can be written as

.

The full nonlinear system to be solved can be written as

with the solution for the uknown vector *u* given by

.

A Newton-based solution for this system of equations, based on a first-order Taylor series expansion about the solution for *u* at iteration *k-1*, can be written as

,

where

is the system Jacobian with individual components give by

and

is the Newton update to be solved for. One method for doing so is by solving

.

The advantage of Newton-based methods is that, with a good initial guess for the solution, convergence rates are very often quadratic (e.g. the residual decreases quadratically, so that at iteration *k* one has a residual of 0.1, at iteration *k+1*, a residual of 0.01, and at iteration *k+2* a residual of 0.0001), whereas Picard-based iterations are much slower to converge. A figure comparing rates of convergence for Picard versus Newton on a CISM tests case is shown below. The Newton method is based on the work of Lemieux et. al (submitted to *JCP*), which is discussed further below.

### The Jacobian Free Approach

In practice, the model Jacobian may either be too difficult or to expensive too form. A "Jacobian Free Newton-Krylov" (JFNK) approach has recently been implemented in CISM (Leimieux et al., submitted to *JCP*), largely following methods discussed in Knoll and Keyes (2004). The crux of the method comes from noting that, when solving the last equation above using a Krylov method (e.g. Conjugate Gradients, GMRES, etc.) the solution for the Newton update is taken from a combination of Krylov vectors that span the subspace

.

This implies that, when using a Krylov method, one only ever needs to calculate matrix vector products of the form when building up the subspace that approximates the solution vector .

Following Knoll and Keyes (2004), note that the necessary matrix vector products can be approximated through nonlinear function evaluations and a perturbation as

.

It is not immediately obvious why this approximation is valid. To verify this, take a few steps back and consider a nonlinear system of equations of two variables, *u1* and *u2*. The right-hand side of the above equation can be expanded as

.

A first-order Taylor series expansion approximation to this is given by

,

which collapses to

.

Finally, note that the right-hand side of the above equation is equal to

,

with the Jacobian matrix given by

.

This matrix vector product is what needs to be calculated repeatedly while building up the Krylov subspace vectors that combine to approximate the Newton update vector . The important point is that at no point in this process does one need to calculate the entire Jacobian matrix. Another important point is that the accuracy of the approximation to the Jacobian is proportional to the small perturbation term, .