Fix up primal vs dual gibbous comparison

demos/gibbous-ice-shelf/Makefile

@@ -1,4 +1,4 @@
-gibbous.pdf volumes.pdf residuals.pdf &: calving.h5 residuals.json make_plots.py
+gibbous.pdf volumes.pdf counts.pdf &: calving.h5 primal-counts.json dual-counts.json make_plots.py
 	python make_plots.py --calving-freq 24.0
 
 calving.h5: steady-state-fine.h5 gibbous.py
@@ -7,8 +7,8 @@ calving.h5: steady-state-fine.h5 gibbous.py
 steady-state-fine.h5: steady-state-coarse.h5 gibbous.py
 	python gibbous.py --resolution 2e3 --final-time 400.0 --num-steps 400 --input $< --output $@
 
-residuals.json: steady-state-coarse.h5 primal_dual_calving.py
-	python primal_dual_calving.py
+%-counts.json: steady-state-coarse.h5 primal_dual_calving.py
+	python primal_dual_calving.py --form $* --hmin 1e-3 --final-time 60 --num-steps 60 --calving-frequency 24 --output $@
 
 steady-state-coarse.h5: gibbous.py
 	python gibbous.py --resolution 5e3 --final-time 400.0 --num-steps 200 --output $@

demos/gibbous-ice-shelf/make_plots.py

@@ -127,11 +127,17 @@
 
 # Make a plot showing the number of Newton iterations required during the
 # calving phase of the experiment
-with open("residuals.json", "r") as residuals_file:
-    residuals = json.load(residuals_file)
+with open("primal-counts.json", "r") as input_file:
+    primal_counts = json.load(input_file)
+
+with open("dual-counts.json", "r") as input_file:
+    dual_counts = json.load(dual_counts)
 
 fig, ax = plt.subplots(figsize=(6.4, 3.2))
-ax.bar(list(range(len(residuals))), [len(r) for r in residuals])
+indices = np.array(list(range(len(primal_counts))))
+width = 0.125
+ax.bar(indices - 2 * width, primal_counts, width=width, label="primal")
+ax.bar(indices, dual_counts, width=width, label="dual")
 ax.set_xlabel("Timestep (years)")
 ax.set_ylabel("Iterations")
-fig.savefig("residuals.pdf", bbox_inches="tight")
+fig.savefig("counts.pdf", bbox_inches="tight")

demos/gibbous-ice-shelf/primal_dual_calving.py

@@ -1,13 +1,25 @@
 import json
+import argparse
 import numpy as np
 from numpy import pi as π
 import tqdm
 import firedrake
 from firedrake import inner, derivative
 import irksome
+import icepack
 from icepack2 import model, solvers
 from icepack2.constants import glen_flow_law as n
 
+parser = argparse.ArgumentParser()
+parser.add_argument("--form", choices=["primal", "dual"])
+parser.add_argument("--hmin", type=float, default=1e-3)
+parser.add_argument("--rtol", type=float, default=1e-4)
+parser.add_argument("--final-time", type=float, default=400.0)
+parser.add_argument("--num-steps", type=int, default=400)
+parser.add_argument("--calving-frequency", type=float, default=24.0)
+parser.add_argument("--output")
+args = parser.parse_args()
+
 # Load the input data
 with firedrake.CheckpointFile("steady-state-coarse.h5", "r") as chk:
     mesh = chk.load_mesh()
@@ -20,45 +32,87 @@
 
 Q = h.function_space()
 V = u.function_space()
-Σ = M.function_space()
-
-Z = V * Σ
-z = firedrake.Function(Z)
-z.sub(0).assign(u)
-z.sub(1).assign(M);
 
 # Set up the Lagrangian
 ε_c = firedrake.Constant(0.01)
 τ_c = firedrake.Constant(0.1)
+h_min = firedrake.Constant(args.hmin)
 
-u, M = firedrake.split(z)
-fields = {
-    "velocity": u,
-    "membrane_stress": M,
-    "thickness": h,
-}
-
-fns = [model.viscous_power, model.ice_shelf_momentum_balance]
-
-rheology = {
-    "flow_law_exponent": n,
-    "flow_law_coefficient": ε_c / τ_c ** n,
+params = {
+    "snes_type": "newtonls",
+    "snes_rtol": args.rtol,
+    "snes_linesearch_type": "nleqerr",
+    "snes_monitor": None,
 }
 
-L = sum(fn(**fields, **rheology) for fn in fns)
-F = derivative(L, z)
-
-# Set up a regularized Lagrangian
-h_min = firedrake.Constant(1e-3)
-rfields = {
-    "velocity": u,
-    "membrane_stress": M,
-    "thickness": firedrake.max_value(h_min, h),
-}
+# Set up the diagnostic solver and boundary conditions and do an initial solve
+if args.form == "primal":
+    from icepack.constants import ice_density as ρ_I, water_density as ρ_W, gravity as g
+    from firedrake import div
+    def gravity(**kwargs):
+        u = kwargs["velocity"]
+        h = kwargs["thickness"]
+        ρ = ρ_I * (1 - ρ_I / ρ_W)
+        return 0.5 * ρ * g * h**2 * div(u)
+
+    def terminus(**kwargs):
+        return firedrake.Constant(0.0)
+
+    _model = icepack.models.IceShelf(gravity=gravity, terminus=terminus)
+
+    u = u.copy(deepcopy=True)
+    A = firedrake.Constant(ε_c / τ_c ** n)
+    opts = {
+        "diagnostic_solver_type": "petsc",
+        "diagnostic_solver_parameters": params,
+    }
+    outer_solver = icepack.solvers.FlowSolver(_model, dirichlet_ids=[1], **opts)
+    outer_solver.diagnostic_solve(velocity=u, thickness=h, fluidity=A)
+    h = outer_solver.fields["thickness"]
+    diagnostic_solver = outer_solver._diagnostic_solver._solver
+
+elif args.form == "dual":
+    Σ = M.function_space()
+    Z = V * Σ
+    z = firedrake.Function(Z)
+    z.sub(0).assign(u)
+    z.sub(1).assign(M);
+
+    u, M = firedrake.split(z)
+    fields = {
+        "velocity": u,
+        "membrane_stress": M,
+        "thickness": h,
+    }
+
+    fns = [model.viscous_power, model.ice_shelf_momentum_balance]
+
+    rheology = {
+        "flow_law_exponent": n,
+        "flow_law_coefficient": ε_c / τ_c ** n,
+    }
+
+    L = sum(fn(**fields, **rheology) for fn in fns)
+    F = derivative(L, z)
+
+    # Set up a regularized Lagrangian
+    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 solver and boundary conditions and do an initial solve
+    bc = firedrake.DirichletBC(Z.sub(0), u_0, (1,))
+    diagnostic_problem = firedrake.NonlinearVariationalProblem(F, z, J=J_r, bcs=bc)
+    params = {"solver_parameters": params}
+    diagnostic_solver = firedrake.NonlinearVariationalSolver(diagnostic_problem, **params)
+    diagnostic_solver.solve()
 
-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 mass balance equation
 prognostic_problem = model.mass_balance(
@@ -69,10 +123,7 @@
     test_function=firedrake.TestFunction(Q),
 )
 
-final_time = 400.0
-num_steps = 400
-
-dt = firedrake.Constant(final_time / num_steps)
+dt = firedrake.Constant(args.final_time / args.num_steps)
 t = firedrake.Constant(0.0)
 method = irksome.BackwardEuler()
 prognostic_params = {
@@ -86,18 +137,6 @@
     prognostic_problem, method, t, dt, h, **prognostic_params
 )
 
-# Set up the diagnostic solver and boundary conditions and do an initial solve
-bc = firedrake.DirichletBC(Z.sub(0), u_0, (1,))
-diagnostic_problem = firedrake.NonlinearVariationalProblem(F, z, J=J_r, bcs=bc)
-params = {
-    "solver_parameters": {
-        "snes_type": "newtonls",
-        "snes_rtol": 1e-4,
-        "snes_linesearch_type": "nleqerr",
-        "snes_monitor": None,
-    },
-}
-diagnostic_solver = firedrake.NonlinearVariationalSolver(diagnostic_problem, **params)
 
 # Create the calving mask -- this describes how we'll remove ice
 radius = firedrake.Constant(60e3)
@@ -106,20 +145,22 @@
 mask = firedrake.conditional(inner(x - y, x - y) < radius**2, 0.0, 1.0)
 
 # Run the simulation
-calving_frequency = 24.0
 time_since_calving = 0.0
 
-for step in range(num_steps):
+num_newton_iterations = []
+for step in range(args.num_steps):
     prognostic_solver.advance()
 
-    if time_since_calving > calving_frequency:
+    if time_since_calving > args.calving_frequency:
         h.interpolate(mask * h)
         time_since_calving = 0.0
     time_since_calving += float(dt)
-    h.interpolate(firedrake.max_value(0, h))
+    expr = firedrake.max_value(0 if args.form == "dual" else h_min, h)
+    h.interpolate(expr)
     diagnostic_solver.solve()
+    num_newton_iterations.append(diagnostic_solver.snes.getIterationNumber())
 
 
 # Save the results to disk
-with open("residuals.json", "w") as residuals_file:
-    json.dump(residuals, residuals_file)
+with open(args.output, "w") as output_file:
+    json.dump(num_newton_iterations, output_file)