First-Order $x-y$ Stokes MMS Test Cases

Component Being Tested

Albany Land-Ice (ALI) model, First-Order Stokes stress-velocity solver \cite{Dukowicz:2010}. This model is included in the land-ice components of the E3SM, and is implemented within the Albany code base \cite{Albany:2016}, available on github at the following URL: https://github.com/SNLComputation/Albany.

Description

It is commontly accepted the land ice behaves like a power-law viscous, incompressible fluid in a low Reynolds number flow, described by the nonlinear Stokes flow equations for glaciers and ice sheets \cite{Dukowicz:2010, Schoof:2010}. Here, we consider the First-Order (FO) approximation to these equations, also referred to as the Blatter–Pattyn model \cite{Pattyn:2003, Blatter:1995}. This model is derived from the full Stokes model under the assumption of a small geometric aspect ratio, and the assumption that the normal vectors to the ice sheet’s upper and lower surfaces are nearly vertical.

3D First-Order Stokes Equations

The First-Order (FO) Stokes equations \cite{Pattyn:2003, Blatter:1995} are given by:

\begin{equation}\label{eq:1} \left\{\begin{array}{rcl} -\nabla \cdot (2 \mu \dot{{\mathbf{ \epsilon}}}_1) &=& -\rho g \frac{\partial s}{\partial x}, \\ -\nabla \cdot (2 \mu \dot{{\mathbf{\epsilon}}}_2) &=& -\rho g \frac{\partial s}{\partial y} \\ \end{array}\right. \end{equation}

Here, in the most general (3D) case,

\begin{equation}\label{eq:2} \dot{{\mathbf{\epsilon}}}_1^T = \left(\begin{array}{ccc} 2\dot{\epsilon}_{xx} + \dot{\epsilon}_{yy}, &\dot{\epsilon}_{xy},& \dot{\epsilon}_{xz}\end{array}\right) \end{equation}\begin{equation}\label{eq:3} \dot{{\mathbf{\epsilon}}}_2^T = \left( \begin{array}{ccc}\dot{\epsilon}_{xy}, & \dot{\epsilon}_{xx} + 2\dot{\epsilon}_{yy}, &\dot{\epsilon}_{yz} \end{array}\right) \end{equation}\begin{equation}\label{eq:4} \dot{\epsilon}_{ij} = \frac{1}{2}\left( \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i}\right) \end{equation}

for $i,j \in \{ x,y,z\}$. The viscosity $\mu$ is given by Glen's law:

\begin{equation}\label{eq:5} \mu = \frac{1}{2} A^{-\frac{1}{n}} \left(\dot{\epsilon}_{xx}^2 + \dot{\epsilon}_{yy}^2 + \dot{\epsilon}_{xx} \dot{\epsilon}_{yy} + \dot{\epsilon}_{xy}^2 + \dot{\epsilon}_{xz}^2 + \dot{\epsilon}_{yz}^2 \right)^{\left( \frac{1}{2n}-\frac{1}{2}\right)} \end{equation}

Here, $A$ is the flow rate factor, and $n$ is the Glen's law exponent ($n=3$ for ice sheets).

2D $x-y$ First-Order Stokes Equations

The MMS problems given in this document are for a simplified version of the first order Stokes equations, namely (\ref{eq:1}) in two spatial dimensions (2D) and on a rectangular geometry aligned with the $x$- and $y$-axes: $\Omega = (0,1) \times (0,1)$ [6]. For this geometry, $s(x,y) = const$, so the gravity-based source terms in (\ref{eq:1}) vanish. Some additional sources $f_1$ and $f_2$ are added to the equations, however. Toward this effect, the governing equations for the MMS problems presented here are the following:

\begin{equation}\label{eq:6} \left\{\begin{array}{rcl} -\nabla \cdot (2 \mu \dot{{\mathbf{ \epsilon}}}_1) +f_1&=& 0, \\ -\nabla \cdot (2 \mu \dot{{\mathbf{\epsilon}}}_2) + f_2&=& 0\\ \end{array}\right. \end{equation}

where

\begin{equation}\label{eq:7} \dot{{\mathbf{\epsilon}}}_1^T = \left(\begin{array}{cc} 2\dot{\epsilon}_{xx} + \dot{\epsilon}_{yy}, &\dot{\epsilon}_{xy}\end{array}\right) \end{equation}\begin{equation}\label{eq:8} \dot{{\mathbf{\epsilon}}}_2^T = \left( \begin{array}{cc}\dot{\epsilon}_{xy}, & \dot{\epsilon}_{xx} + 2\dot{\epsilon}_{yy} \end{array}\right) \end{equation}

and

\begin{equation}\label{eq:9} \mu = \frac{1}{2} A^{-\frac{1}{n}} \left(\dot{\epsilon}_{xx}^2 + \dot{\epsilon}_{yy}^2 + \dot{\epsilon}_{xx} \dot{\epsilon}_{yy} + \dot{\epsilon}_{xy}^2 \right)^{\left( \frac{1}{2n}-\frac{1}{2}\right)} \end{equation}

with $\dot{\epsilon}_{ij}$ given by (\ref{eq:4}). Written out, (\ref{eq:6}) has the form:

\begin{equation}\label{eq:10} \left\{\begin{array}{rcl} -\frac{\partial}{\partial x} \left( 4 \mu \frac{\partial u}{\partial x} + 2\mu \frac{\partial v}{\partial y}\right) - \frac{\partial}{\partial y} \left( \mu \frac{\partial u}{\partial y} + \mu \frac{\partial v}{\partial x}\right) + f_1&=& 0, \\ -\frac{\partial}{\partial x} \left( \mu \frac{\partial u}{\partial x} + \mu \frac{\partial v}{\partial y}\right) - \frac{\partial}{\partial y} \left(2\mu \frac{\partial u}{\partial x} + 4\mu \frac{\partial v}{\partial y} \right) +f_2 &=& 0\\ \end{array}\right. \end{equation}

SinCos2D MMS Test Case

The first MMS test case derived will be referred to as the Sincos2D test case. It is an Dirichlet/homogeneous Neumann BVP on the unit square: $\Omega = (0,1)\times (0,1)$. The force terms in (\ref{eq:6}) and boundary conditions are derived such that the exact solution to the problem is given by:

\begin{equation}\label{eq:11} \begin{array}{rcl} u &= & \sin(2 \pi x + \phi) \cos(2\pi y + \psi) + 3\pi x \\ v &=& -\cos(2\pi x + \phi) \sin(2\pi y + \psi) - 3 \pi y \end{array} \end{equation}

The parameters $\phi, \psi \in [0, 2\pi)$ in (\ref{eq:11}) are phase shifts that can be used to generate a family of solutions. In order to obtain this exact solution, the source terms in (\ref{eq:6}) are given by:

\begin{equation}\label{eq:12} f_1 = -8\mu \pi^2 \sin(2\pi x + \phi)\cos(2\pi y + \psi) + A^{-\frac{1}{n}} \left(\frac{1}{n} - 1\right)\tilde{\mu}^{\frac{1}{n} - 2}\left(\frac{\partial \tilde{\mu}}{\partial x}(2\epsilon_{xx} + \epsilon_{yy}) + \frac{\partial \tilde{\mu}}{\partial y}\epsilon_{xy}\right) \end{equation}\begin{equation}\label{eq:13} f_2 = 8\mu\pi^2\cos(2\pi x + \phi)\sin(2\pi y + \psi) + A^{-\frac{1}{n}} \left(\frac{1}{n} - 1\right)\tilde{\mu}^{\frac{1}{n} - 2}\left(\frac{\partial \tilde{\mu}}{\partial x} \epsilon_{xy} + \frac{\partial \tilde{\mu}}{\partial y}(\epsilon_{xx} + 2\epsilon_{yy})\right) \end{equation}

where

\begin{equation}\label{eq:14} \tilde{\mu} = 2\pi\cos(2\pi x + \phi)\cos(2\pi y + \psi) +3\pi x \end{equation}\begin{equation}\label{eq:15} \frac{\partial \tilde{\mu} }{\partial x} = -4\pi^2\sin(2\pi x + \phi)\cos(2\pi y + \psi) \end{equation}\begin{equation}\label{eq:16} \frac{\partial \tilde{\mu}}{\partial y} = -4\pi^2\cos(2\pi x + \phi)\sin(2\pi y + \psi) \end{equation}\begin{equation}\label{eq:17} \epsilon_{xx} = 2\pi\cos(2\pi x + \phi)\cos(2\pi y + \psi) + 3\pi \end{equation}\begin{equation}\label{eq:18} \epsilon_{yy} = -2\pi\cos(2\pi x + \phi)\cos(2\pi y + \psi) - 3\pi \end{equation}\begin{equation}\label{eq:19} \epsilon_{xy} = 0 \end{equation}\begin{equation}\label{eq:20} \mu = \frac{1}{2}A^{-\frac{1}{n}} \tilde{\mu}^{\frac{1}{n} - 1} \end{equation}

The following piece of code specifies $f_1$ and $f_2$ (taken from the Albany code base):

       double r = 3.0*pi;
       MeshScalarT x2pi = 2.0*pi*coordVec(cell,qp,0);
       MeshScalarT y2pi = 2.0*pi*coordVec(cell,qp,1);
       MeshScalarT muargt = 2.0*pi*cos(x2pi + xphase)*cos(y2pi + yphase) + r;
       MeshScalarT muqp = 0.5*pow(A, -1.0/n)*pow(muargt, 1.0/n - 1.0);
       MeshScalarT dmuargtdx = -4.0*pi*pi*sin(x2pi + xphase)*cos(y2pi + yphase);
       MeshScalarT dmuargtdy = -4.0*pi*pi*cos(x2pi + xphase)*sin(y2pi + yphase);
       MeshScalarT exx = 2.0*pi*cos(x2pi + xphase)*cos(y2pi + yphase) + r;
       MeshScalarT eyy = -2.0*pi*cos(x2pi + xphase)*cos(y2pi + yphase) - r;
       MeshScalarT exy = 0.0;
       f[0] = 2.0*muqp*(-4.0*pi*pi*sin(x2pi + xphase)*cos(y2pi + yphase))
            + 2.0*0.5*pow(A, -1.0/n)*(1.0/n - 1.0)*pow(muargt, 1.0/n
            - 2.0)*(dmuargtdx*(2.0*exx + eyy) + dmuargtdy*exy);
       f[1] = 2.0*muqp*(4.0*pi*pi*cos(x2pi + xphase)*sin(y2pi + yphase))
            + 2.0*0.5*pow(A, -1.0/n)*(1.0/n - 1.0)*pow(muargt, 1.0/n
            - 2.0)*(dmuargtdx*exy + dmuargtdy*(exx + 2.0*eyy));

To finish the problem specification, we still need to give the boundary conditions. They are:

\begin{equation}\label{eq:21} \begin{array}{rcll} u&=& 0, & \text{at } x=0 \\ u &=& 3\pi, & \text{at } x = 1\\ 2\mu \dot{{\bf{ \epsilon}}}_2 \cdot \mathbf{n} &=& 0, & \text{at } y=0, y=1 \\ 2\mu \dot{{\bf{\epsilon}}}_1 \cdot \mathbf{n} &=& 0, & \text{at } x=0, x=1\\ v &=& 0, & \text{at } y = 0 \\ v &=& -3\pi, & \text{at } y = 1 \end{array} \end{equation}

Note that the Sincos2Dtest case can be turned into a constant viscosity test case by setting $n = 1$.

Variable Definitions

The input file for the SinCos2D test case is named input_fo_sincos2DT.yaml and can be found in the tests/small/LandIce/FO_MMS directory of the Albany code. Here, the following parameter values are utilized:

$$ \begin{array}{l|l} \text{Variable} & \text{Value} \\ \hline n & 3 \\ A & 1.0 \\ \phi & 0.0 \\ \psi & 0.0 \\ \end{array} $$

Discussion

The exact solution to the SinCos2D problem is plotted below.

In [ ]:
from mpl_toolkits.mplot3d.axes3d import Axes3D
import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
import numpy as np


figU = plt.figure()
axU = figU.gca(projection='3d')


figV = plt.figure()
axV = figV.gca(projection='3d')

phi = 0.0
psi = 0.0

# Make data.
X = np.arange(0.0, 1.0, 0.01)
Y = np.arange(0.0, 1.0, 0.01)
X, Y = np.meshgrid(X, Y)
U = np.sin(2*np.pi*X + phi)*np.cos(2*np.pi*Y + psi) + 3*np.pi*X
V = -np.cos(2*np.pi*X + phi)*np.sin(2*np.pi*Y + psi) - 3*np.pi*Y

# Plot the surface.
surf = axU.plot_surface(X, Y, U, cmap=cm.coolwarm,
                       linewidth=0, antialiased=False)

figU.suptitle('u-velocity', fontsize=16)
axU.set_xlabel('x')
axU.set_ylabel('y')

surf = axV.plot_surface(X, Y, V, cmap=cm.coolwarm,
                       linewidth=0, antialiased=False)

figV.suptitle('v-velocity', fontsize=16)
axV.set_xlabel('x')
axV.set_ylabel('y')

The figure below shows the expected and observed convergence rates for several different finite elements.

In [ ]:
from IPython.display import Image
Image(filename='fo_sincos2D_convrate_plot_slides.png')

CosExp2D MMS Test Case

The second MMS test case derived will be referred to as the CosExp2D test case, and is detailed in [6]. The problem specification includes Robin boundary conditions resembling the basal sliding boundary condition. For this problem, the force terms in (\ref{eq:6}) and boundary conditions are derived such that the exact solution to the problem is given by:

\begin{equation}\label{eq:22} \begin{array}{rcl} u &= & e^x \sin(2\pi y) \\ v &=& e^x \cos(2 \pi y) \end{array} \end{equation}

In order to obtain this exact solution, the source terms in (\ref{eq:6}) are given by:

\begin{equation}\label{eq:23} f_1 = 2\mu\exp(x)\sin(2\pi y)\left[2 - 3\pi - 2\pi^2\right] + A^{-\frac{1}{n}}\left(\frac{1}{n}-1\right)\tilde{\mu}^{ \frac{1}{n}-2}\left(\frac{\partial \tilde{\mu}}{\partial x}(2\epsilon_{xx} + \epsilon_{yy}) + \frac{\partial \tilde{\mu}}{\partial y}\epsilon_{xy}\right) \end{equation}\begin{equation}\label{eq:24} f_2 = 2\mu\exp(x)\cos(2\pi y)\left[3\pi + \frac{1}{2} - 8\pi^2\right] + A^{-\frac{1}{n}}\left(\frac{1}{n}-1\right)\tilde{\mu}^{\frac{1}{n}-2} \left(\frac{\partial \tilde{\mu}}{\partial x}\epsilon_{xy} + \frac{\partial \tilde{\mu}}{\partial y}(\epsilon_{xx} + 2\epsilon_{yy})\right) \end{equation}

where

\begin{equation}\label{eq:25} \tilde{\mu} = \exp(x) \sqrt{(1 + 4\pi^2 - 2\pi)\sin^2(2\pi y) + \frac{1}{4}(2\pi+1)^2\cos^2(2\pi y)} \end{equation}\begin{equation}\label{eq:26} \frac{\partial \tilde{\mu} }{\partial x} = \tilde{\mu} \end{equation}\begin{equation}\label{eq:27} \frac{\partial \tilde{\mu}}{\partial y} = \frac{3}{2}\frac{\pi(1+4\pi^2-4\pi)\cos(2\pi y)\sin(2\pi y)\exp(x)}{\sqrt{(1 + 4\pi^2 - 2\pi)\sin^2(2\pi y) + \frac{1}{4}(2\pi+1)^2\cos^2(2\pi y)}} \end{equation}\begin{equation}\label{eq:28} \epsilon_{xx} = \exp(x)\sin(2\pi y) \end{equation}\begin{equation}\label{eq:29} \epsilon_{yy} = -2\pi\exp(x)\sin(2\pi y) \end{equation}\begin{equation}\label{eq:30} \epsilon_{xy} = \frac{1}{2}(2\pi+1)\exp(x)\cos(2\pi y) \end{equation}\begin{equation}\label{eq:31} \mu = \frac{1}{2}A^{-\frac{1}{n}} \tilde{\mu}^{\frac{1}{n} - 1} \end{equation}

The following piece of code specifies $f_1$ and $f_2$ (taken from the Albany code base):

       double a = 1.0;
       MeshScalarT x2pi = 2.0*pi*coordVec(cell,qp,0);
       MeshScalarT x = coordVec(cell,qp,0);
       MeshScalarT y2pi = 2.0*pi*coordVec(cell,qp,1);
       MeshScalarT muargt = (a*a + 4.0*pi*pi - 2.0*pi*a)*sin(y2pi)*sin(y2pi)
                          + 1.0/4.0*(2.0*pi+a)*(2.0*pi+a)*cos(y2pi)*cos(y2pi);
       muargt = sqrt(muargt)*exp(a*x);
       MeshScalarT muqp = 1.0/2.0*pow(A, -1.0/n)*pow(muargt, 1.0/n - 1.0);
       MeshScalarT dmuargtdx = a*muargt;
       MeshScalarT dmuargtdy = 3.0/2.0*pi*(a*a+4.0*pi*pi
                             -4.0*pi*a)*cos(y2pi)*sin(y2pi)*exp(a*x)/sqrt((a*a
                             + 4.0*pi*pi - 2.0*pi*a)*sin(y2pi)*sin(y2pi)
                             + 1.0/4.0*(2.0*pi+a)*(2.0*pi+a)*cos(y2pi)*cos(y2pi));
       MeshScalarT exx = a*exp(a*x)*sin(y2pi);
       MeshScalarT eyy = -2.0*pi*exp(a*x)*sin(y2pi);
       MeshScalarT exy = 1.0/2.0*(2.0*pi+a)*exp(a*x)*cos(y2pi);
       f[0] = 2.0*muqp*(2.0*a*a*exp(a*x)*sin(y2pi) - 3.0*pi*a*exp(a*x)*sin(y2pi)
            - 2.0*pi*pi*exp(a*x)*sin(y2pi))
            + 2.0*0.5*pow(A, -1.0/n)*(1.0/n-1.0)*pow(muargt, 1.0/n-2.0)*(dmuargtdx*(2.0*exx
            + eyy) + dmuargtdy*exy);
       f[1] = 2.0*muqp*(3.0*a*pi*exp(a*x)*cos(y2pi) + 1.0/2.0*a*a*exp(a*x)*cos(y2pi)
            - 8.0*pi*pi*exp(a*x)*cos(y2pi)) + 2.0*0.5*pow(A, -1.0/n)*(1.0/n-1.0)*pow(muargt,
            1.0/n-2.0)*(dmuargtdx*exy + dmuargtdy*(exx + 2.0*eyy));

To finish the problem specification, we still need to give the boundary conditions. They are:

\begin{equation}\label{eq:32} \begin{array}{rcll} u &=& 0, & \text{at } y=0, 1 \\ v &=& 0, & \text{at } (x,y) = (0,0) \\ 2\mu \dot{{\mathbf{ \epsilon}}}_1 \cdot \mathbf{n} &=& 4(\pi-1)\mu u, & \text{at } x=0 \\ 2\mu \dot{{\mathbf{ \epsilon}}}_1 \cdot \mathbf{n} &=& -4(\pi -1) \mu u, & \text{at } x=1 \\ 2\mu \dot{{\mathbf{ \epsilon}}}_2 \cdot \mathbf{n} &=& -(2\pi + 1) v, & \text{at } x=0\\ 2\mu \dot{{\mathbf{ \epsilon}}}_2 \cdot \mathbf{n} &=& (2\pi + 1) v, & \text{at } x=1\\ 2\mu \dot{{\mathbf{ \epsilon}}}_2 \cdot \mathbf{n} &=& 0, & \text{at } y=0, 1 \end{array} \end{equation}

The boundary condition at the single point $(x,y) = (0,0)$ for $v$ is specified to make $v$ unique (otherwise, $v$ could be non-unique up to a constant).

Note that the CosExp2D test case can be turned into a constant viscosity test case by setting $n = 1$.

Variable Definitions

The input file for the CosExp2D test case is named input_fo_cosexp_basal_all_glensLawT.yaml and can be found in the tests/small/LandIce/FO_MMS directory of the Albany code. This is the same problem detailed in Section 4.1 of [6]. The following parameter values are utilized:

$$ \nonumber \begin{array}{l|l} \text{Variable} & \text{Value} \\ \hline n & 3 \\ A & 1.0 \\ \end{array} $$

Discussion

The exact solution to the CosExp2D problem is plotted below.

In [2]:
from mpl_toolkits.mplot3d.axes3d import Axes3D
import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
import numpy as np


figU = plt.figure()
axU = figU.gca(projection='3d')


figV = plt.figure()
axV = figV.gca(projection='3d')

# Make data.
X = np.arange(0.0, 1.0, 0.01)
Y = np.arange(0.0, 1.0, 0.01)
X, Y = np.meshgrid(X, Y)
U = np.exp(X)*np.sin(2*np.pi*Y)
V = np.exp(X)*np.cos(2*np.pi*Y)

# Plot the surface.
surf = axU.plot_surface(X, Y, U, cmap=cm.coolwarm,
                       linewidth=0, antialiased=False)

figU.suptitle('u-velocity', fontsize=16)
axU.set_xlabel('x')
axU.set_ylabel('y')

surf = axV.plot_surface(X, Y, V, cmap=cm.coolwarm,
                       linewidth=0, antialiased=False)

figV.suptitle('v-velocity', fontsize=16)
axV.set_xlabel('x')
axV.set_ylabel('y')
Out[2]:
Text(0.5, 0, 'y')

The figure below shows the expected and observed convergence rates for several different finite elements (from [6]).

In [1]:
from IPython.display import Image
Image(filename='fo_cosexp_basal_convrate_plot.png')
Out[1]:
In [ ]:
In [ ]:

References

[1] Dukowicz J., Price S. and Lipscomb W., ``Consistent approximations and boundary conditions for ice-sheet dynamics from a principle of least action'', Journal of Glaciology, vol. 56, number , pp. 480-496, 2010.

[2] Salinger A., Bartlett R., Bradley A. et al., ``Albany: Using Agile Components to Develop a Flexible, Generic Multiphysics Analysis Code'', Int. J. Multiscale Comput. Engng, vol. 14, number , pp. 415-438, 2016.

[3] Schoof C. and Hindmarsh R., ``Thin-Film Flows withWall Slip: An Asymptotic Analysis of Higher Order Glacier Flow Models'', Q.J. Mech. Appl. Math, vol. 63, number , pp. 73-114, 2010.

[4] Pattyn F., ``A new three-dimensional higher-order thermomechanical ice-sheet model: basic sensitivity, ice stream development, and ice flow across subglacial lakes'', J. Geophys. Res., vol. 108, number , pp. 1-15, 2003.

[5] Blatter H., ``Velocity and stress fields in grounded glaciers: a simple algorithm for including deviatoric stress gradients'', J. Glaciology, vol. 41, number , pp. 333-344, 1995.