Numerical solution of DIVA / SSA

The DIVA and SSA momentum balances share the same 2D depth-integrated form: a nonlinear elliptic problem for the depth-averaged horizontal velocity $\bar{\mathbf u} = (\bar u, \bar v)$, with coefficients $(\mu, \beta_\mathrm{eff})$ that depend on the solution. Yelmo discretises this problem on the Arakawa C-grid (the ac-grid) and linearises it by Picard iteration, freezing $(\mu, \beta_\mathrm{eff}, H)$ at the previous iterate, solving a linear system for the new $\bar{\mathbf u}$, relaxing, and repeating until convergence.

Two solvers are provided for the linear step. Both produce the same solution at interior cells (the two systems are related by a sign and a cell-area factor), but they differ in the structure of the assembled matrix and in how boundary conditions are imposed. The choice is made at runtime via the parameter ydyn.ssa_solver = "residual" | "energy", which the DIVA solver dispatches on.

Discretisation conventions

Yelmo uses the Arakawa C / staggered AC grid:

• aa-nodes (cell centres): scalars — $H$, $s$, $\bar\mu$, $\bar\mu\,H$, $\beta$ (before staggering).
• acx-nodes (right face of each aa-cell): $\bar u$, $\tau_{d,x}$, $\beta_{\mathrm{eff},x}$.
• acy-nodes (top face): $\bar v$, $\tau_{d,y}$, $\beta_{\mathrm{eff},y}$.
• ab-nodes (cell corners): cross-coupling viscosity $(\bar\mu H)^{\mathrm{ab}}$, obtained by averaging $\bar\mu\,H$ from its four neighbouring aa-cells (stagger_visc_aa_ab).

The two unknowns per cell ($\bar u$ at the right face and $\bar v$ at the top face) are interleaved into a single state vector $\mathbf x = [\,\bar u_{1}, \bar v_{1}, \bar u_{2}, \bar v_{2}, \dots\,]$, so the assembled linear system has dimension $2 N_\mathrm{cells}$. The matrix is stored in CSR format and solved with LIS (Library of Iterative Solvers for Linear Systems); the iterative method and preconditioner are configured at runtime via ydyn.ssa_lis_opt.

Per-row solver masks (ssa_mask_acx, ssa_mask_acy) classify each ac-node as one of:

• 0: Dirichlet zero velocity (frozen bed, edge of domain, etc.);
• -1: Dirichlet prescribed velocity (e.g. observed input);
• 1: solve (interior);
• 3: lateral / calving-front Neumann row driven by the depth-integrated lateral stress $\tau_{l,\mathrm{int}}$.

The two assemblers share the same argument list so they are drop-in interchangeable from the Picard loop in calc_velocity_diva.

Solver A — residual form (Yelmo v1)

The residual assembler (solver_ssa_ac.f90) builds a non-symmetric matrix $A_\mathrm{res}$ and right-hand side $\mathbf b_\mathrm{res}$ directly from the strong form of the SSA PDE. Writing $N \equiv \bar\mu\,H$ for brevity, the interior $\bar u$ row is a 5-point stencil in $\bar u$ plus cross-coupling terms in $\bar v$ on the neighbouring acy-faces, balanced against the basal drag and the driving stress. Schematically (cell $(i,j)$, with $N^{\mathrm{aa}}$ on aa-nodes and $N^{\mathrm{ab}}$ on ab-nodes):

$$\sum_{(i',j')}\alpha^{x}_{i',j'}\,\bar u_{i',j'} \;+\; \sum_{(i',j')}\gamma^{xy}_{i',j'}\,\bar v_{i',j'} \;-\; \beta_\mathrm{eff}\,\bar u_{i,j} \;=\; \tau_{d,x}(i,j).$$

The $\bar u$ stencil coefficients carry the membrane stretching and shearing,

$$\alpha^{x}_{i-1,j} \;=\; \tfrac{4}{\Delta x^{2}}\,N^{\mathrm{aa}}_{i,j},$$

$$\alpha^{x}_{i+1,j} \;=\; \tfrac{4}{\Delta x^{2}}\,N^{\mathrm{aa}}_{i+1,j},$$

$$\alpha^{x}_{i,j-1} \;=\; \tfrac{1}{\Delta y^{2}}\,N^{\mathrm{ab}}_{i,j-1},$$

$$\alpha^{x}_{i,j+1} \;=\; \tfrac{1}{\Delta y^{2}}\,N^{\mathrm{ab}}_{i,j},$$

with the centre coefficient $\alpha^{x}_{i,j}$ being minus the sum of its four neighbours. The $\gamma^{xy}$ coefficients couple the row to four neighbouring $\bar v$-faces via the $\partial \bar v/\partial x$ and $\partial \bar v/\partial y$ cross-terms in the membrane stress. The $\bar v$ row has the analogous structure.

The right-hand side is the driving stress itself (no cell-area factor). The system $A_\mathrm{res}\,\mathbf x = \mathbf b_\mathrm{res}$ is non-symmetric, and is solved with a Krylov method such as BiCGStab plus an algebraic preconditioner.

• Lateral / calving-front rows substitute the membrane-stress balance with the prescribed depth-integrated lateral stress $\tau_{l,\mathrm{int}}$ in a row-specific stencil.
• Dirichlet rows replace the equation by $\bar u_{i,j} = u^*$ with the corresponding column kept in place (non-symmetric).
• Free-slip, no-slip and periodic conditions are applied per side via the boundary-code helper get_neighbor_indices_bc_codes.

This is the formulation inherited from Yelmo v1 and is the legacy default.

Solver B — energy form (new)

The energy assembler (solver_ssa_ac_energy.f90) builds the Hessian of a discrete energy functional and solves $K\,\mathbf x = \mathbf b$ for the velocity that minimises that energy. With $(\mu, \beta, H)$ frozen during each Picard step the energy is quadratic, so $K = \frac{\partial^{2} W}{\partial \mathbf x^{\,2}}$ is symmetric positive (semi-)definite and the linear step can use a symmetric Krylov method — CG with an AMG preconditioner — in place of BiCGStab.

Energy density

The continuum energy density underlying the SSA momentum balance is the sum of a membrane (deformation) term, a basal-drag term, and a gravitational potential-energy term:

$$\begin{aligned} W \;=\; \bar\mu\,H\,&\biggl( 2\!\left(\frac{\partial \bar u}{\partial x}\right)^{\!2} + 2\!\left(\frac{\partial \bar v}{\partial y}\right)^{\!2} + 2\,\frac{\partial \bar u}{\partial x}\,\frac{\partial \bar v}{\partial y} + \tfrac{1}{2}\!\left(\frac{\partial \bar u}{\partial y} + \frac{\partial \bar v}{\partial x}\right)^{\!2} \,\biggr) \\[4pt] &+\; \tfrac{1}{2}\,\beta\,(\bar u^{\,2} + \bar v^{\,2}) \;+\; \rho_i\,g\,H\,\!\left(\bar u\,\frac{\partial s}{\partial x} + \bar v\,\frac{\partial s}{\partial y}\right). \end{aligned}$$

The first line is $2\,\bar\mu\,H\,\dot{\bar\varepsilon}_{ij}\,\dot{\bar\varepsilon}_{ij}$ written out for the depth-averaged horizontal strain rates. Stationarity of the integral $\mathcal W = \int W \,\mathrm dx\,\mathrm dy$ with respect to $(\bar u, \bar v)$ reproduces exactly the SSA / DIVA strong form, so any critical point of $\mathcal W$ is a solution of the momentum balance.

Discrete assembly

Yelmo's discrete energy is the cell-by-cell evaluation of $W$ on the C-grid, with the derivatives $\frac{\partial \bar u}{\partial x}, \frac{\partial \bar v}{\partial y}, \frac{\partial \bar u}{\partial y}, \frac{\partial \bar v}{\partial x}$ expressed as the natural finite differences between adjacent ac-nodes. The Hessian $K$ then has the same stencil graph as the residual matrix but is symmetric in $(\bar u, \bar v)$. At inner cells the two formulations are related by an exact algebraic identity (documented in the header of the energy assembler):

$$K_\mathrm{inner} \;=\; -\,A_\mathrm{res, inner}\cdot \Delta x\,\Delta y, \qquad \mathbf b_\mathrm{inner} \;=\; -\,\boldsymbol\tau_d \cdot \Delta x\,\Delta y.$$

So the energy formulation is, at interior cells, the residual formulation rescaled by the cell area and a sign — the physical solution at interior cells is identical to machine precision. Where the two solvers differ is at boundaries:

• Lateral / calving-front BC: the front stress enters the energy as a boundary-work term $\pm\,\tau_{l,\mathrm{int}}\,\Delta y$ on the RHS, with the sign determined by the outward normal. This is the variational form of the Neumann condition and is symmetric by construction.
• Dirichlet rows: prescribed values are imposed by static condensation — the prescribed column is multiplied by the known velocity and moved to the RHS — instead of by row replacement. This preserves the symmetry of $K$ and so allows CG / AMG to be used without spoiling the SPD structure.

The viscosity staggering aa $\to$ ab is the same routine (stagger_visc_aa_ab) used by the residual assembler.

Why bother?

Two practical advantages flow from the SPD structure:

1. Linear solver choice: CG with AMG is typically faster and more robust than BiCGStab + ILU for large, well-conditioned SPD systems, and converges monotonically in the energy norm.
2. Physical interpretability and discrete consistency: every term in $K$ and $\mathbf b$ corresponds to a contribution to a discrete energy. Boundary conditions that are natural for the continuum functional (e.g. Neumann front stress) become natural for the discrete one. This makes it easier to add new physics — e.g. alternative friction laws, additional body forces — in a way that provably preserves the variational structure.

The Picard loop, the viscosity update, the F-integral closure and the basal-stress and 3D velocity diagnostics are identical between the two solvers: only the linear-system assembly differs.