Use Rainier demo as unit test

Author: danshapero

test/dome_test.py

@@ -0,0 +1,217 @@
+import numpy as np
+from numpy import pi as π
+import firedrake
+from firedrake import inner, grad, dx, ds, exp, min_value, max_value, Constant
+import irksome
+from icepack2.model import mass_balance
+from icepack2.model.variational import momentum_balance, flow_law, friction_law
+from icepack2.constants import gravity, ice_density, glen_flow_law
+
+
+def run_simulation(refinement_level: int):
+    radius = Constant(12e3)
+    mesh = firedrake.UnitDiskMesh(refinement_level)
+    mesh.coordinates.dat.data[:] *= float(radius)
+
+    # Make a bunch of finite elements and function spaces
+    cg1 = firedrake.FiniteElement("CG", "triangle", 1)
+    dg1 = firedrake.FiniteElement("DG", "triangle", 1)
+    dg0 = firedrake.FiniteElement("DG", "triangle", 0)
+    S = firedrake.FunctionSpace(mesh, cg1)
+    Q = firedrake.FunctionSpace(mesh, dg1)
+    V = firedrake.VectorFunctionSpace(mesh, cg1)
+    Σ = firedrake.TensorFunctionSpace(mesh, dg0, symmetry=True)
+    T = firedrake.VectorFunctionSpace(mesh, cg1)
+    Z = V * Σ * T
+    W = V * Σ * T * Q
+
+    x = firedrake.SpatialCoordinate(mesh)
+
+    # Make the bed topography
+    B = Constant(4e3)
+    r_b = Constant(150e3 / (2 * π))
+    expr = B * exp(-inner(x, x) / r_b**2)
+    b = firedrake.Function(S).interpolate(expr)
+
+    # Make the mass balance field
+    z_measured = Constant(1600.0)
+    a_measured = Constant(-0.917 * 8.7)
+    a_top = Constant(0.7)
+    z_top = Constant(4e3)
+    δa_δz = (a_top - a_measured) / (z_top - z_measured)
+    a_max = Constant(0.7)
+
+    def smb(z):
+        return min_value(a_max, a_measured + δa_δz * (z - z_measured))
+
+    # Make the initial thickness
+    r_h = Constant(5e3)
+    H = Constant(100.0)
+    expr = H * firedrake.max_value(0, 1 - inner(x, x) / r_h**2)
+    h = firedrake.Function(Q).interpolate(expr)
+    h_0 = h.copy(deepcopy=True)
+
+    s = firedrake.Function(Q).interpolate(b + h)
+    a = firedrake.Function(Q).interpolate(smb(s))
+
+    # Fluidity of ice in yr⁻¹ MPa⁻³ at 0C
+    A = Constant(158.0)
+
+    # Make an initial guess for the velocity using SIA
+    u = firedrake.Function(V)
+    v = firedrake.TestFunction(V)
+
+    ρ_I = Constant(ice_density)
+    g = Constant(gravity)
+
+    n = Constant(glen_flow_law)
+
+    P = ρ_I * g * h
+    S_n = inner(grad(s), grad(s))**((n - 1) / 2)
+    u_shear = -2 * A * P ** n / (n + 2) * h * S_n * grad(s)
+    F = inner(u - u_shear, v) * dx
+
+    degree = 1
+    qdegree = max(8, degree ** glen_flow_law)
+    pparams = {"form_compiler_parameters": {"quadrature_degree": qdegree}}
+    firedrake.solve(F == 0, u, **pparams)
+
+    # Compute the initial velocity using the dual form of SSA
+    z = firedrake.Function(Z)
+    z.sub(0).assign(u);
+
+    τ_c = Constant(0.1)
+    ε_c = Constant(A * τ_c ** n)
+
+    K = h * A / (n + 2)
+    U_c = Constant(100.0)
+    u_c = K * τ_c ** n + U_c
+
+    glen_rheology = {
+        "flow_law_exponent": n,
+        "flow_law_coefficient": ε_c / τ_c ** n,
+        "sliding_exponent": n,
+        "sliding_coefficient": u_c / τ_c ** n,
+    }
+
+    α = firedrake.Constant(1e-4)
+    linear_rheology = {
+        "flow_law_exponent": 1,
+        "flow_law_coefficient": ε_c / τ_c,
+        "sliding_exponent": 1,
+        "sliding_coefficient": u_c / τ_c,
+    }
+
+    u, M, τ = firedrake.split(z)
+    fields = {
+        "velocity": u,
+        "membrane_stress": M,
+        "basal_stress": τ,
+        "thickness": h,
+        "surface": s,
+    }
+
+    sparams = {
+        "solver_parameters": {
+            "snes_monitor": ":rainier-output-initial.log",
+            "snes_type": "newtonls",
+            "snes_max_it": 200,
+            "snes_linesearch_type": "nleqerr",
+            "ksp_type": "gmres",
+            "pc_type": "lu",
+            "pc_factor_mat_solver_type": "mumps",
+        },
+    }
+
+    v, N, σ = firedrake.TestFunctions(Z)
+    F = (
+        momentum_balance(**fields, test_function=v)
+        + firedrake.replace(flow_law(**fields, **glen_rheology, test_function=N), {h: H})
+        + α * firedrake.replace(flow_law(**fields, **linear_rheology, test_function=N), {h: H})
+        + friction_law(**fields, **glen_rheology, test_function=σ)
+        + α * friction_law(**fields, **linear_rheology, test_function=σ)
+    )
+
+    momentum_problem = firedrake.NonlinearVariationalProblem(F, z, **pparams)
+    momentum_solver = firedrake.NonlinearVariationalSolver(momentum_problem, **sparams)
+
+    num_continuation_steps = 5
+    for exponent in np.linspace(1.0, 3.0, num_continuation_steps):
+        n.assign(exponent)
+        momentum_solver.solve()
+
+    # Time-dependent solve
+    w = firedrake.Function(W)
+    w.sub(0).assign(z.sub(0))
+    w.sub(1).assign(z.sub(1))
+    w.sub(2).assign(z.sub(2))
+    w.sub(3).assign(h);
+
+    u, M, τ, h = firedrake.split(w)
+    v, N, σ, η = firedrake.TestFunctions(W)
+
+    fields = {
+        "velocity": u,
+        "membrane_stress": M,
+        "basal_stress": τ,
+        "thickness": h,
+        "surface": b + h,
+    }
+
+    F_momentum = (
+        momentum_balance(**fields, test_function=v)
+        + firedrake.replace(flow_law(**fields, **glen_rheology, test_function=N), {h: H})
+        + α * firedrake.replace(flow_law(**fields, **linear_rheology, test_function=N), {h: H})
+        + friction_law(**fields, **glen_rheology, test_function=σ)
+        + α * friction_law(**fields, **linear_rheology, test_function=σ)
+    )
+    F_mass = mass_balance(thickness=h, velocity=u, accumulation=a, test_function=η)
+    F = F_momentum + F_mass
+
+    tableau = irksome.BackwardEuler()
+    t = Constant(0.0)
+    dt = Constant(1.0 / 6)
+
+    lower = firedrake.Function(W)
+    upper = firedrake.Function(W)
+    lower.assign(-np.inf)
+    upper.assign(+np.inf)
+    lower.subfunctions[3].assign(0.0)
+    bounds = ("stage", lower, upper)
+
+    bparams = {
+        "solver_parameters": {
+            "snes_monitor": ":rainier-output-vi.log",
+            "snes_type": "vinewtonrsls",
+            "snes_max_it": 200,
+            "ksp_type": "gmres",
+            "pc_type": "lu",
+            "pc_factor_mat_solver_type": "mumps",
+        },
+        "stage_type": "value",
+        "basis_type": "Bernstein",
+        "bounds": bounds,
+    }
+
+    solver = irksome.TimeStepper(F, tableau, t, dt, w, **bparams, **pparams)
+
+    us = [w.subfunctions[0].copy(deepcopy=True)]
+    hs = [w.subfunctions[3].copy(deepcopy=True)]
+
+    final_time = 300.0
+    num_steps = int(final_time / float(dt))
+    for step in num_steps:
+        solver.advance()
+        h = w.subfunctions[3]
+        a.interpolate(smb(b + h))
+
+        us.append(w.subfunctions[0].copy(deepcopy=True))
+        hs.append(w.subfunctions[3].copy(deepcopy=True))
+
+    return hs, us
+
+
+def test_dome_problem():
+    hs, us = run_simulation(4)
+    volumes = [firedrake.assemble(h * dx) for h in hs]
+    assert volumes[0] <= volumes[-1]