Use icepack2 library in demos

demos/convergence-tests/ice_shelf.py

@@ -3,8 +3,13 @@
 import numpy as np
 import firedrake
 from firedrake import Constant
-from icepack.constants import ice_density, water_density, gravity, glen_flow_law
-from dualform import ice_shelf
+from icepack2.constants import (
+    ice_density as ρ_I,
+    water_density as ρ_W,
+    gravity as g,
+    glen_flow_law as n,
+)
+from icepack2 import model
 
 
 parser = argparse.ArgumentParser()
@@ -15,15 +20,9 @@
 parser.add_argument("--output", default="results.json")
 args = parser.parse_args()
 
-# Physical constants
-# TODO: Copy the values from icepack
-ρ_I = Constant(ice_density)
-ρ_W = Constant(water_density)
+# Physical constants at about -14C, with an applied stress of 100 kPa, the
+# glacier will experience a strain rate of 10 (m / yr) / km at -14C.
 ρ = ρ_I * (1 - ρ_I / ρ_W)
-g = Constant(gravity)
-
-# At about -14C, with an applied stress of 100 kPa, the glacier will experience
-# a strain rate of 10 (m / yr) / km at -14C.
 τ_c = Constant(0.1)  # MPa
 ε_c = Constant(0.01) # (km / yr) / km
 
@@ -34,7 +33,6 @@ def exact_velocity(x, inflow_velocity, inflow_thickness, thickness_change, lengt
     dh = Constant(thickness_change)
     lx = Constant(length)
 
-    n = Constant(glen_flow_law)
     A = ε_c / τ_c**n
     u_0 = Constant(100.0)
     h_0, dh = Constant(500.0), Constant(100.0)
@@ -79,19 +77,24 @@ def exact_velocity(x, inflow_velocity, inflow_thickness, thickness_change, lengt
 
     z = firedrake.Function(Z)
     u, M = firedrake.split(z)
-    kwargs = {
+    fields = {
         "velocity": u,
         "membrane_stress": M,
         "thickness": h,
-        "viscous_yield_strain": ε_c,
-        "viscous_yield_stress": τ_c,
-        "outflow_ids": outflow_ids,
     }
-    fns = [ice_shelf.viscous_power, ice_shelf.boundary, ice_shelf.constraint]
-    J_l = sum(fn(**kwargs, flow_law_exponent=1) for fn in fns)
+    rheology = {
+        "flow_law_exponent": n,
+        "flow_law_coefficient": ε_c / τ_c**n,
+    }
+    linear_rheology = {
+        "flow_law_exponent": 1,
+        "flow_law_coefficient": ε_c / τ_c,
+    }
+    fns = [model.viscous_power, model.ice_shelf_momentum_balance]
+    J_l = sum(fn(**fields, **linear_rheology) for fn in fns)
     F_l = firedrake.derivative(J_l, z)
 
-    J = sum(fn(**kwargs, flow_law_exponent=glen_flow_law) for fn in fns)
+    J = sum(fn(**fields, **rheology) for fn in fns)
     F = firedrake.derivative(J, z)
 
     inflow_bc = firedrake.DirichletBC(Z.sub(0), u_in, inflow_ids)

demos/convergence-tests/ice_stream.py

@@ -3,15 +3,15 @@
 import tqdm
 import numpy as np
 import firedrake
-from firedrake import interpolate, as_vector, max_value, Constant
-from icepack.constants import (
+from firedrake import interpolate, as_vector, max_value, Constant, derivative
+from icepack2.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
+from icepack2 import model
 
 
 parser = argparse.ArgumentParser()
@@ -54,7 +54,7 @@ def friction(x):
     return α * (ρ_I * g * (h0 - dh * x / Lx)) * exact_u(x) ** (-1 / m)
 
 
-errors = []
+errors, mesh_sizes = [], []
 k_min, k_max, num_steps = args.log_nx_min, args.log_nx_max, args.num_steps
 for nx in np.logspace(k_min, k_max, num_steps, base=2, dtype=int):
     mesh = firedrake.RectangleMesh(nx, nx, float(Lx), float(Ly), diagonal="crossed")
@@ -65,8 +65,7 @@ def friction(x):
     Q = firedrake.FunctionSpace(mesh, cg)
     V = firedrake.VectorFunctionSpace(mesh, cg)
     Σ = firedrake.TensorFunctionSpace(mesh, dg, symmetry=True)
-    # TODO: investigate using DG for this?
-    T = firedrake.VectorFunctionSpace(mesh, cg)
+    T = firedrake.VectorFunctionSpace(mesh, dg)
     Z = V * Σ * T
     z = firedrake.Function(Z)
     z.sub(0).assign(Constant((u_inflow, 0)))
@@ -81,67 +80,80 @@ def friction(x):
     C = interpolate(friction(x), Q)
     u_c = interpolate((τ_c / C)**m, Q)
 
+    # Create the boundary conditions
     inflow_ids = (1,)
     outflow_ids = (2,)
     side_wall_ids = (3, 4)
 
+    inflow_bc = firedrake.DirichletBC(Z.sub(0), Constant((u_inflow, 0)), inflow_ids)
+    side_wall_bc = firedrake.DirichletBC(Z.sub(0).sub(1), 0, side_wall_ids)
+    bcs = [inflow_bc, side_wall_bc]
+
+    # Get the Lagrangian, material parameters, and input fields
+    fns = [
+        model.viscous_power,
+        model.friction_power,
+        model.calving_terminus,
+        model.momentum_balance,
+    ]
+
+    rheology = {
+        "flow_law_exponent": n,
+        "flow_law_coefficient": ε_c / τ_c ** n,
+        "sliding_exponent": m,
+        "sliding_coefficient": u_c / τ_c ** m,
+    }
+
     u, M, τ = firedrake.split(z)
-    kwargs = {
+    fields = {
         "velocity": u,
         "membrane_stress": M,
         "basal_stress": τ,
         "thickness": h,
         "surface": s,
-        "viscous_yield_strain": ε_c,
-        "viscous_yield_stress": τ_c,
-        "friction_yield_speed": u_c,
-        "friction_yield_stress": τ_c,
-        "outflow_ids": outflow_ids,
     }
 
-    inflow_bc = firedrake.DirichletBC(Z.sub(0), Constant((u_inflow, 0)), inflow_ids)
-    side_wall_bc = firedrake.DirichletBC(Z.sub(0).sub(1), 0, side_wall_ids)
-    bcs = [inflow_bc, side_wall_bc]
+    boundary_ids = {"outflow_ids": outflow_ids}
 
-    fns = [
-        ice_stream.viscous_power,
-        ice_stream.friction_power,
-        ice_stream.boundary,
-        ice_stream.constraint,
-    ]
+    L = sum(fn(**fields, **rheology, **boundary_ids) for fn in fns)
+    F = derivative(L, z)
+
+    # Regularize second derivative of the Lagrangian
+    λ = firedrake.Constant(1e-3)
+    regularized_rheology = {
+        "flow_law_exponent": 1,
+        "flow_law_coefficient": λ * ε_c / τ_c,
+        "sliding_exponent": 1,
+        "sliding_coefficient": λ * u_c / τ_c,
+    }
 
-    J_l = sum(
-        fn(**kwargs, flow_law_exponent=1, sliding_law_exponent=1) for fn in fns
+    K = sum(
+        fn(**fields, **regularized_rheology)
+        for fn in [model.viscous_power, model.friction_power]
     )
-    F_l = firedrake.derivative(J_l, z)
-    firedrake.solve(F_l == 0, z, bcs=bcs)
+    J = derivative(derivative(L + K, z), z)
 
     qdegree = max(8, args.degree ** n)
-    params = {
-        "form_compiler_parameters": {"quadrature_degree": qdegree}
-    }
+    pparams = {"form_compiler_parameters": {"quadrature_degree": qdegree}}
+    problem = firedrake.NonlinearVariationalProblem(F, z, bcs, J=J, **pparams)
+    solver = firedrake.NonlinearVariationalSolver(problem)
 
-    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=bcs, **params)
+    solver.solve()
+    λ.assign(0.0)
+    solver.solve()
 
     u, M, τ = z.subfunctions
     error = firedrake.norm(u - u_exact) / firedrake.norm(u_exact)
     δx = mesh.cell_sizes.dat.data_ro.min()
-    errors.append((δx, error))
+    mesh_sizes.append(δx)
+    errors.append(error)
 
 try:
     with open(args.output, "r") as input_file:
         results = json.load(input_file)
 except FileNotFoundError:
     results = {}
 
-results.update({f"degree-{args.degree}": errors})
+results.update({f"degree-{args.degree}": list(zip(mesh_sizes, errors))})
 with open(args.output, "w") as output_file:
     json.dump(results, output_file)

demos/gibbous-ice-shelf/gibbous.py

@@ -10,11 +10,13 @@
     grad,
     dx,
     ds,
+    derivative,
     NonlinearVariationalProblem,
     NonlinearVariationalSolver,
 )
-from icepack.constants import glen_flow_law
-from dualform import ice_shelf
+import irksome
+from icepack2 import model
+from icepack2.constants import glen_flow_law as n
 
 parser = argparse.ArgumentParser()
 parser.add_argument("--input")
@@ -86,10 +88,11 @@
 
 # Create some function spaces and some fields
 cg = firedrake.FiniteElement("CG", "triangle", 1)
-dg = firedrake.FiniteElement("DG", "triangle", 0)
-Q = firedrake.FunctionSpace(mesh, cg)
+dg0 = firedrake.FiniteElement("DG", "triangle", 0)
+dg1 = firedrake.FiniteElement("DG", "triangle", 1)
+Q = firedrake.FunctionSpace(mesh, dg1)
 V = firedrake.VectorFunctionSpace(mesh, cg)
-Σ = firedrake.TensorFunctionSpace(mesh, dg, symmetry=True)
+Σ = firedrake.TensorFunctionSpace(mesh, dg0, symmetry=True)
 Z = V * Σ
 
 h0 = firedrake.interpolate(h_expr, Q)
@@ -129,40 +132,34 @@
 else:
     z.sub(0).assign(u0)
 
-# Set up the diagnostic problem and compute an initial guess by solving a
-# Picard linearization of the problem
+# Set up the diagnostic problem
 u, M = firedrake.split(z)
 fields = {
     "velocity": u,
     "membrane_stress": M,
     "thickness": h,
 }
 
-h_min = firedrake.Constant(1.0)
-rfields = {
-    "velocity": u,
-    "membrane_stress": M,
-    "thickness": firedrake.max_value(h_min, h),
-}
+fns = [model.viscous_power, model.ice_shelf_momentum_balance]
 
-params = {
-    "viscous_yield_strain": ε_c,
-    "viscous_yield_stress": τ_c,
-    "outflow_ids": (2,),
+rheology = {
+    "flow_law_exponent": n,
+    "flow_law_coefficient": ε_c / τ_c ** n,
 }
 
-fns = [ice_shelf.viscous_power, ice_shelf.boundary, ice_shelf.constraint]
-L_1 = sum(fn(**fields, **params, flow_law_exponent=1) for fn in fns)
-F_1 = firedrake.derivative(L_1, z)
+L = sum(fn(**fields, **rheology) for fn in fns)
+F = derivative(L, z)
 
-L = sum(fn(**fields, **params, flow_law_exponent=glen_flow_law) for fn in fns)
-F = firedrake.derivative(L, z)
+# Make an initial guess by solving a Picard linearization
+linear_rheology = {
+    "flow_law_exponent": 1,
+    "flow_law_coefficient": ε_c / τ_c,
+}
 
-L_r = sum(fn(**rfields, **params, flow_law_exponent=glen_flow_law) for fn in fns)
-F_r = firedrake.derivative(L_r, z)
-J_r = firedrake.derivative(F_r, z)
+L_1 = sum(fn(**fields, **linear_rheology) for fn in fns)
+F_1 = derivative(L_1, z)
 
-qdegree = int(glen_flow_law) + 2
+qdegree = n + 2
 bc = firedrake.DirichletBC(Z.sub(0), u0, (1,))
 problem_params = {
     #"form_compiler_parameters": {"quadrature_degree": qdegree},
@@ -179,7 +176,17 @@
 }
 firedrake.solve(F_1 == 0, z, **problem_params, **solver_params)
 
-# Create a regularized Lagrangian
+# Create an approximate linearization
+h_min = firedrake.Constant(1.0)
+rfields = {
+    "velocity": u,
+    "membrane_stress": M,
+    "thickness": firedrake.max_value(h_min, h),
+}
+
+L_r = sum(fn(**rfields, **rheology) for fn in fns)
+F_r = derivative(L_r, z)
+J_r = derivative(F_r, z)
 
 # Set up the diagnostic problem and solver
 # TODO: possibly use a regularized Jacobian
@@ -188,16 +195,26 @@
 diagnostic_solver.solve()
 
 # Set up the prognostic problem and solver
-h_n = h.copy(deepcopy=True)
-φ = firedrake.TestFunction(Q)
+prognostic_problem = model.mass_balance(
+    thickness=h,
+    velocity=u,
+    accumulation=firedrake.Constant(0.0),
+    thickness_inflow=h0,
+    test_function=firedrake.TestFunction(Q),
+)
 dt = firedrake.Constant(args.final_time / args.num_steps)
-flux_cells = ((h - h_n) / dt * φ - inner(h * u, grad(φ))) * dx
-ν = firedrake.FacetNormal(mesh)
-flux_in = h0 * firedrake.min_value(0, inner(u, ν)) * φ * ds
-flux_out = h * firedrake.max_value(0, inner(u, ν)) * φ * ds
-G = flux_cells + flux_in + flux_out
-prognostic_problem = NonlinearVariationalProblem(G, h)
-prognostic_solver = NonlinearVariationalSolver(prognostic_problem)
+t = firedrake.Constant(0.0)
+method = irksome.BackwardEuler()
+prognostic_params = {
+    "solver_parameters": {
+        "snes_type": "ksponly",
+        "ksp_type": "gmres",
+        "pc_type": "bjacobi",
+    },
+}
+prognostic_solver = irksome.TimeStepper(
+    prognostic_problem, method, t, dt, h, **prognostic_params
+)
 
 # Set up a calving mask
 R = firedrake.Constant(60e3)
@@ -214,7 +231,7 @@
     chk.save_function(h, name="thickness", idx=0)
 
     for step in tqdm.trange(args.num_steps):
-        prognostic_solver.solve()
+        prognostic_solver.advance()
 
         if args.calving_freq != 0.0:
             if time_since_calving > args.calving_freq:
@@ -223,7 +240,6 @@
             time_since_calving += float(dt)
 
         h.interpolate(firedrake.max_value(0, h))
-        h_n.assign(h)
 
         diagnostic_solver.solve()