Add mass balance equation

demos/larsen/simulate.py

@@ -3,7 +3,7 @@
 import geojson
 import xarray
 import firedrake
-from firedrake import assemble, Constant, inner, grad, dx
+from firedrake import assemble, Constant, inner, grad, dx, ds
 import icepack
 from icepack.constants import glen_flow_law as n
 from dualform import ice_shelf
@@ -51,7 +51,7 @@
 μ_cutoff = Constant(0.05)
 h.interpolate(firedrake.conditional(μ < μ_cutoff, 0, h))
 
-# Set up the model and solver
+# Set up the momentum balance equation and solver
 T = Constant(260)
 A0 = icepack.rate_factor(T)
 
@@ -103,6 +103,8 @@
 solver_params = {
     "solver_parameters": {
         "snes_monitor": None,
+        "snes_max_it": 10,
+        "snes_convergence_test": "skip",
         "snes_type": "newtonls",
         "ksp_type": "gmres",
         "pc_type": "lu",
@@ -115,9 +117,37 @@
 u_solver = firedrake.NonlinearVariationalSolver(u_problem, **solver_params)
 u_solver.solve()
 
+# Set up the mass balance equation
+h_n = h.copy(deepcopy=True)
+h0 = h.copy(deepcopy=True)
+φ = firedrake.TestFunction(Q)
+final_time = 4.0
+num_steps = 48
+dt = Constant(final_time / 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
+h_problem = firedrake.NonlinearVariationalProblem(G, h)
+h_solver = firedrake.NonlinearVariationalSolver(h_problem)
+
+for step in range(num_steps):
+    h_solver.solve()
+    h.interpolate(firedrake.max_value(0, h))
+    h_n.assign(h)
+    u_solver.solve()
+
 import matplotlib.pyplot as plt
 fig, axes = plt.subplots()
 axes.set_aspect("equal")
-colors = firedrake.tripcolor(z.sub(0), axes=axes)
+colors = firedrake.tripcolor(h, axes=axes)
+fig.colorbar(colors)
+plt.show()
+
+fig, axes = plt.subplots()
+axes.set_aspect("equal")
+δh = firedrake.interpolate(h - h0, Q)
+colors = firedrake.tripcolor(δh, vmin=-50, vmax=+50, cmap="RdBu_r", axes=axes)
 fig.colorbar(colors)
 plt.show()