◆ assemble_matrix()

void pism::stressbalance::SSAFD::assemble_matrix	(	const Inputs &	inputs,
		bool	include_basal_shear,
		Mat	A
	)

protectedvirtual

Assemble the left-hand side matrix for the KSP-based, Picard iteration, and finite difference implementation of the SSA equations.

Recall the SSA equations are

\begin{align*} - 2 \left[\nu H \left(2 u_x + v_y\right)\right]_x - \left[\nu H \left(u_y + v_x\right)\right]_y - \tau_{(b)1} &= - \rho g H h_x, \\ - \left[\nu H \left(u_y + v_x\right)\right]_x - 2 \left[\nu H \left(u_x + 2 v_y\right)\right]_y - \tau_{(b)2} &= - \rho g H h_y, \end{align*}

where \(u\) is the \(x\)-component of the velocity and \(v\) is the \(y\)-component of the velocity.

The coefficient \(\nu\) is the vertically-averaged effective viscosity. (The product \(\nu H\) is computed by compute_nuH_staggered().) The Picard iteration idea is that, to solve the nonlinear equations in which the effective viscosity depends on the velocity, we freeze the effective viscosity using its value at the current estimate of the velocity and we solve the linear equations which come from this viscosity. In abstract symbols, the Picard iteration replaces the above nonlinear SSA equations by a sequence of linear problems

\[ A(U^{(k)}) U^{(k+1)} = b \]

where \(A(U)\) is the matrix from discretizing the SSA equations supposing the viscosity is a function of known velocities \(U\), and where \(U^{(k)}\) denotes the \(k\)th iterate in the process. The current method assembles \(A(U)\).

For ice shelves \(\tau_{(b)i} = 0\) [MacAyealetal]. For ice streams with a basal till modelled as a plastic material, \(\tau_{(b)i} = - \tau_c u_i/|\mathbf{u}|\) where \(\mathbf{u} = (u,v)\), \(|\mathbf{u}| = \left(u^2 + v^2\right)^{1/2}\), where \(\tau_c(t,x,y)\) is the yield stress of the till [SchoofStream]. More generally, ice streams can be modeled with a pseudo-plastic basal till; see IceModel::initBasalTillModel() and IceModel::updateYieldStressUsingBasalWater() and reference [BKAJS]. The pseudo-plastic till model includes all power law sliding relations and the linearly-viscous model for sliding, \(\tau_{(b)i} = - \beta u_i\) where \(\beta\) is the basal drag (friction) parameter [MacAyeal]. In any case, PISM assumes that the basal shear stress can be factored this way, even if the coefficient depends on the velocity, \(\beta(u,v)\). Such factoring is possible even in the case of (regularized) plastic till. This scalar coefficient \(\beta\) is what is returned by IceBasalResistancePlasticLaw::drag().

Note that the basal shear stress appears on the left side of the linear system we actually solve. We believe this is crucial, because of its effect on the spectrum of the linear approximations of each stage. The effect on spectrum is clearest in the linearly-viscous till case but there seems to be an analogous effect in the plastic till case.

This method assembles the matrix for the left side of the above SSA equations. The numerical method is finite difference. Suppose we use difference notation \(\delta_{+x}f^{i,j} = f^{i+1,j}-f^{i,j}\), \(\delta_{-x}f^{i,j} = f^{i,j}-f^{i-1,j}\), and \(\Delta_{x}f^{i,j} = f^{i+1,j}-f^{i-1,j}\), and corresponding notation for \(y\) differences, and that we write \(N = \nu H\) then the first of the two "concrete" SSA equations above has this discretization:

\begin{align*} - &2 \frac{N^{i+\frac{1}{2},j}}{\Delta x} \left[2\frac{\delta_{+x}u^{i,j}}{\Delta x} + \frac{\Delta_{y} v^{i+1,j} + \Delta_{y} v^{i,j}}{4 \Delta y}\right] + 2 \frac{N^{i-\frac{1}{2},j}}{\Delta x} \left[2\frac{\delta_{-x}u^{i,j}}{\Delta x} + \frac{\Delta_y v^{i,j} + \Delta_y v^{i-1,j}}{4 \Delta y}\right] \\ &\qquad- \frac{N^{i,j+\frac{1}{2}}}{\Delta y} \left[\frac{\delta_{+y} u^{i,j}}{\Delta y} + \frac{\Delta_x v^{i,j+1} + \Delta_x v^{i,j}}{4 \Delta x}\right] + \frac{N^{i,j-\frac{1}{2}}}{\Delta y} \left[\frac{\delta_{-y}u^{i,j}}{\Delta y} + \frac{\Delta_x v^{i,j} + \Delta_x v^{i,j-1}}{4 \Delta x}\right] - \tau_{(b)1}^{i,j} = - \rho g H^{i,j} \frac{\Delta_x h^{i,j}}{2\Delta x}. \end{align*}

As a picture, see Figure ssastencil.

ssastencil: Stencil for our finite difference discretization of the first of the two scalar SSA equations. Triangles show staggered grid points where N = nu * H is evaluated. Circles and squares show where u and v are approximated, respectively.

It follows immediately that the matrix we assemble in the current method has 13 nonzeros entries per row because, for this first SSA equation, there are 5 grid values of \(u\) and 8 grid values of \(v\) used in this scheme. For the second equation we also have 13 nonzeros per row.

FIXME: document use of DAGetMatrix and MatStencil and MatSetValuesStencil

Definition at line 491 of file SSAFD.cc.

Referenced by picard_manager().