Checkpoint, start on MISMIP test case

demos/mismip/mismip.py

@@ -0,0 +1,134 @@
+import argparse
+import numpy as np
+import firedrake
+from firedrake import (
+    sqrt, exp, min_value, max_value, inner, as_vector, Constant, dx, interpolate
+)
+from icepack.constants import (
+    ice_density as ρ_I,
+    water_density as ρ_W,
+    gravity as g,
+    glen_flow_law as n,
+    weertman_sliding_law as m,
+)
+from dualform import ice_stream
+
+
+lx, ly = 640e3, 80e3
+Lx, Ly = Constant(lx), Constant(ly)
+ny = 20
+nx = int(lx / ly) * ny
+area = lx * ly
+
+mesh = firedrake.RectangleMesh(nx, ny, lx, ly, diagonal="crossed")
+
+# Set up some function spaces. TODO: Make it work for higher degree.
+cg = firedrake.FiniteElement("CG", "triangle", 1)
+dg = firedrake.FiniteElement("DG", "triangle", 0)
+
+Q = firedrake.FunctionSpace(mesh, cg)
+V = firedrake.VectorFunctionSpace(mesh, cg)
+Σ = firedrake.TensorFunctionSpace(mesh, dg, symmetry=True)
+T = firedrake.VectorFunctionSpace(mesh, cg)
+Z = V * Σ * T
+
+# Set up the basal elevation.
+x, y = firedrake.SpatialCoordinate(mesh)
+
+x_c = Constant(300e3)
+X = x / x_c
+
+B_0 = Constant(-150)
+B_2 = Constant(-728.8)
+B_4 = Constant(343.91)
+B_6 = Constant(-50.57)
+B_x = B_0 + B_2 * X**2 + B_4 * X**4 + B_6 * X**6
+
+f_c = Constant(4e3)
+d_c = Constant(500)
+w_c = Constant(24e3)
+
+B_y = d_c * (
+    1 / (1 + exp(-2 * (y - Ly / 2 - w_c) / f_c)) +
+    1 / (1 + exp(+2 * (y - Ly / 2 + w_c) / f_c))
+)
+
+z_deep = Constant(-720)
+b = interpolate(max_value(B_x + B_y, z_deep), Q)
+
+# Some physical constants; `A = ε_c / τ_c ** n`, `C = τ_c / u_c ** (1 / m)`.
+# Rather than work with the fluidity `A` and friction `C` directly, we use
+# these stress, strain rate, and velocity scales so that we can easily rescale
+# the physical constants under changing flow law and sliding exponents.
+ε_c = Constant(0.02)
+τ_c = Constant(0.1)
+u_c = Constant(1000.0)
+
+a = Constant(0.3)
+
+# Set up the boundary conditions.
+inflow_ids = (1,)
+outflow_ids = (2,)
+side_wall_ids = (3, 4)
+
+inflow_bc = firedrake.DirichletBC(Z.sub(0), Constant((0, 0)), inflow_ids)
+side_wall_bc = firedrake.DirichletBC(Z.sub(0).sub(1), 0, side_wall_ids)
+bcs = [inflow_bc, side_wall_bc]
+
+# Set up the solution variables, input data, Lagrangian, and solvers.
+fns = [
+    ice_stream.viscous_power,
+    ice_stream.friction_power,
+    ice_stream.boundary,
+    ice_stream.constraint,
+]
+
+z = firedrake.Function(Z)
+δu = Constant(90)
+z.sub(0).interpolate(as_vector((δu * x / Lx, 0)))
+
+u, M, τ = firedrake.split(z)
+h = interpolate(Constant(100.0), Q)
+s = interpolate(max_value(b + h, (1 - ρ_I / ρ_W) * h), Q)
+h_min = Constant(10.0)
+p_I = ρ_I * g * max_value(h_min, h)
+p_W = -ρ_W * g * min_value(0, s - h)
+f = (1 - p_W / p_I) ** m
+kwargs = {
+    "velocity": u,
+    "membrane_stress": M,
+    "basal_stress": τ,
+    "thickness": h,
+    "surface": s,
+    "floating": f,
+    "viscous_yield_strain": ε_c,
+    "viscous_yield_stress": τ_c,
+    "friction_yield_speed": u_c,
+    "friction_yield_stress": τ_c,
+    "outflow_ids": outflow_ids,
+}
+
+linear_exponents = {"flow_law_exponent": 1, "sliding_law_exponent": 1}
+J_l = sum(fn(**kwargs, **linear_exponents) for fn in fns)
+F_l = firedrake.derivative(J_l, z)
+firedrake.solve(F_l == 0, z, bcs)
+
+num_continuation_steps = 4
+for t in np.linspace(0.0, 1.0, num_continuation_steps):
+    exponents = {
+        "flow_law_exponent": Constant((1 - t) + t * n),
+        "sliding_law_exponent": Constant((1 - t) + t * m),
+    }
+    J = sum(fn(**kwargs, **exponents) for fn in fns)
+    F = firedrake.derivative(J, z)
+    firedrake.solve(F == 0, z, bcs)
+
+import matplotlib.pyplot as plt
+u, M, τ = z.subfunctions
+fig, axes = plt.subplots()
+axes.set_aspect("equal")
+axes.get_xaxis().set_visible(False)
+axes.get_yaxis().set_visible(False)
+colors = firedrake.tripcolor(u, axes=axes)
+fig.colorbar(colors, orientation="horizontal")
+plt.show()