Add figure for number of Newton iterations

demos/gibbous-ice-shelf/Makefile

@@ -1,4 +1,4 @@
-gibbous.pdf volumes.pdf &: calving.h5 make_plots.py
+gibbous.pdf volumes.pdf residuals.pdf &: calving.h5 residuals.json make_plots.py
 	python make_plots.py --calving-freq 24.0
 
 calving.h5: steady-state-fine.h5 gibbous.py
@@ -7,5 +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
+
 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

@@ -1,4 +1,5 @@
 import argparse
+import json
 import numpy as np
 import matplotlib.pyplot as plt
 from matplotlib_scalebar.scalebar import ScaleBar
@@ -122,3 +123,15 @@
 ax.annotate("fine\nspin-up", xy=(t1, tlevel), xytext=((t1 + t2) / 2, tlevel), **kwargs)
 ax.annotate("calving\ncycle", xy=(t2, tlevel), xytext=((t2 + t3) / 2, tlevel), **kwargs)
 fig.savefig("volumes.pdf", bbox_inches="tight")
+
+
+# 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)
+
+fig, ax = plt.subplots(figsize=(6.4, 3.2))
+ax.bar(list(range(len(residuals))), [len(r) for r in residuals])
+ax.set_xlabel("Timestep (years)")
+ax.set_ylabel("Iterations")
+fig.savefig("residuals.pdf", bbox_inches="tight")

demos/gibbous-ice-shelf/primal_dual_calving.py

@@ -0,0 +1,118 @@
+import json
+import numpy as np
+from numpy import pi as π
+import tqdm
+import firedrake
+from firedrake import inner, derivative
+import irksome
+from icepack2 import model, solvers
+from icepack2.constants import glen_flow_law as n
+
+# Load the input data
+with firedrake.CheckpointFile("steady-state-coarse.h5", "r") as chk:
+    mesh = chk.load_mesh()
+    idx = len(chk.h5pyfile["timesteps"]) - 1
+    h = chk.load_function(mesh, "thickness", idx=idx)
+    h0 = h.copy(deepcopy=True)
+    u = chk.load_function(mesh, "velocity", idx=idx)
+    M = chk.load_function(mesh, "membrane_stress", idx=idx)
+
+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)
+
+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
+h_min = firedrake.Constant(0.1)
+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)
+H_r = derivative(F_r, z)
+
+# Set up the mass balance equation
+prognostic_problem = model.mass_balance(
+    thickness=h,
+    velocity=u,
+    accumulation=firedrake.Constant(0.0),
+    thickness_inflow=h0,
+    test_function=firedrake.TestFunction(Q),
+)
+
+final_time = 400.0
+num_steps = 400
+
+dt = firedrake.Constant(final_time / num_steps)
+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 the diagnostic solver and boundary conditions and do an initial solve
+bc = firedrake.DirichletBC(Z.sub(0), firedrake.Constant((0, 0)), (1,))
+problem = solvers.ConstrainedOptimizationProblem(L, z, H=H_r, bcs=bc)
+diagnostic_solver = solvers.NewtonSolver(problem, tolerance=1e-4)
+
+residuals = [list(diagnostic_solver.solve())]
+
+# Create the calving mask -- this describes how we'll remove ice
+radius = firedrake.Constant(60e3)
+x = firedrake.SpatialCoordinate(mesh)
+y = firedrake.Constant((0.0, radius))
+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 tqdm.trange(num_steps):
+    prognostic_solver.advance()
+
+    if time_since_calving > calving_frequency:
+        h.interpolate(mask * h)
+        time_since_calving = 0.0
+    time_since_calving += float(dt)
+    h.interpolate(firedrake.max_value(0, h))
+
+    residuals.append(list(diagnostic_solver.solve()))
+
+# Save the results to disk
+with open("residuals.json", "w") as residuals_file:
+    json.dump(residuals, residuals_file)

dual-problems.tex

@@ -583,21 +583,29 @@ \subsection{Gibbous ice shelf} \label{sec:gibbous-ice-shelf}
 These results were consistent across different mesh resolutions and when run several times on multiple machines.
 Since the dual problem has 4x as many unknowns, the added cost that we found experimentally is less than what we would expect if we naively assumed that cost is proportional to the number of degrees of freedom.
 
-In the calving phase of the experiment, our solver for the dual problem still performed gracefully in ice-free areas.
+In the calving phase of the experiment, our solver for the dual problem still worked in ice-free areas.
 This feature offers the possibility of implementing physically-based calving models in a simple way.
 Figure \ref{fig:gibbous-calving-volumes} shows the evolution of the volume of ice in the shelf over the two spin-up phases and the calving phase.
 The 24-year recurrence interval is not enough time for the ice to advance back to the original edge of the computational domain.
 Using a longer interval would allow the calving terminus to advance back to its original position.
 Figure \ref{fig:gibbous}(d) shows the thickness of the ice shelf immediately after the calving event (e) and (f) show the magnitude of the change in speed and stress respectively.
 In particular, the stress field shows a discontinuity at the new calving front as expected.
 
-The number of nonlinear solver iterations to recompute the ice velocity after the calving event is much greater than after a normal timestep.
+\begin{figure}[h]
+    \begin{center}
+        \includegraphics[width=0.75\linewidth]{demos/gibbous-ice-shelf/residuals.pdf}
+    \end{center}
+    \caption{The number of Newton iterations to compute the ice velocity at each step of the calving phase of the experiment.
+    Calving occurs every 24 years.}
+    \label{fig:gibbous-residuals}
+\end{figure}
+
+Figure \ref{fig:gibbous-residuals} shows the number of Newton iterations needed to compute the ice velocity at each timestep.
+The number of iterations after a calving event is much greater than after a normal timestep.
 This behavior is to be expected if we run the solver to convergence because it introduces a type of shock into the system.
 As the system relaxes back, the number of iterations decreases again.
-Moreover, we had to do some manual adjustment of the convergence tolerances.
 Several different strategies can alleviate the need for manual adjustment (see the appendix).
-\textcolor{red}{Give a plot showing the number of Newton solves.
-Do a comparison to the primal form with e.g. a small minimum thickness.}
+\textcolor{red}{Do a comparison to the primal form with e.g. a small minimum thickness.}
 
 \subsection{Larsen C Ice Shelf}