Add variational form of ice shelf momentum balance

Author: danshapero

src/icepack2/model/variational.py

@@ -78,3 +78,10 @@ def momentum_balance(**kwargs):
     return cell_balance + facet_balance
 
 
+def ice_shelf_momentum_balance(**kwargs):
+    field_names = ("membrane_stress", "thickness", "test_function")
+    M, h, v = map(kwargs.get, field_names)
+    ε = sym(grad(v))
+
+    ρ = ρ_I * (1 - ρ_I / ρ_W)
+    return (-h * inner(M, ε) + 0.5 * ρ * g * h**2 * div(v)) * dx