Revert changes to mass balance, use discontinuous velocity

The changes to the mass balance function were making the tests fail
because we don't actually want the factor of `h` inside the `max_value`
that defines the upwind flux. We can get some of the effect using a
discontinuous velocity field but all this does is change the thickness
downstream to its upstream value * (u / (u + u_c)) which is not zero.

Author: danshapero

notebooks/frontal-ablation.ipynb

@@ -26,8 +26,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "nx, ny = 32, 32\n",
-    "mesh = firedrake.UnitSquareMesh(nx, ny, diagonal=\"crossed\")\n",
+    "nx, ny = 64, 64\n",
+    "mesh = firedrake.UnitSquareMesh(nx, ny, quadrilateral=True)\n",
     "δ = 1.0 / nx"
    ]
   },
@@ -39,7 +39,7 @@
    "outputs": [],
    "source": [
     "degree = 0\n",
-    "element = firedrake.FiniteElement(\"DG\", \"triangle\", degree)\n",
+    "element = firedrake.FiniteElement(\"DQ\", \"quadrilateral\", degree)\n",
     "Q = firedrake.FunctionSpace(mesh, element)"
    ]
   },
@@ -90,12 +90,11 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "u_γ = Constant(-1.0)\n",
-    "x_γ = Constant(0.5)\n",
-    "h_γ = Constant(0.5)\n",
+    "U_c = Constant(-1.0)\n",
+    "x_c = Constant(0.5)\n",
     "\n",
-    "f_γ = h_γ * firedrake.as_vector(\n",
-    "    (firedrake.conditional(x[0] >= x_γ, u_γ, 0), 0)\n",
+    "u_c = firedrake.as_vector(\n",
+    "    (firedrake.conditional(x[0] >= x_c, U_c, 0), 0)\n",
     ")"
    ]
   },
@@ -109,10 +108,9 @@
     "ϕ = firedrake.TestFunction(Q)\n",
     "F = model.mass_balance(\n",
     "    thickness=h,\n",
-    "    velocity=u,\n",
+    "    velocity=u - u_c,\n",
     "    accumulation=a + m,\n",
     "    thickness_inflow=h_in,\n",
-    "    frontal_ablation=f_γ,\n",
     "    test_function=ϕ,\n",
     ")"
    ]
@@ -135,10 +133,9 @@
     "\n",
     "params = {\n",
     "    \"solver_parameters\": {\n",
+    "        \"snes_monitor\": \":frontal-ablation.log\",\n",
+    "        \"snes_atol\": 1e-12,\n",
     "        \"snes_type\": \"vinewtonrsls\",\n",
-    "        \"ksp_type\": \"gmres\",\n",
-    "        \"pc_type\": \"lu\",\n",
-    "        \"pc_factor_mat_solver_type\": \"mumps\",\n",
     "    },\n",
     "    \"stage_type\": \"value\",\n",
     "    \"basis_type\": \"Bernstein\",\n",
@@ -157,7 +154,7 @@
    "source": [
     "hs = [h.copy(deepcopy=True)]\n",
     "\n",
-    "final_time = 1.0\n",
+    "final_time = 2.0\n",
     "num_steps = int(final_time / float(dt))\n",
     "for step in trange(num_steps):\n",
     "    solver.advance()\n",
@@ -203,6 +200,28 @@
    "source": [
     "HTML(animation.to_html5_video())"
    ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "e05e0c73-52bc-4c80-9cfc-0411dc53b568",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "hs[-1].at((0.75, 0.5))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "602b42eb-5285-4838-87af-ee1ed9f4f030",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from firedrake import dx\n",
+    "expr = firedrake.conditional(x[0] <= x_c, h_in, h_in / 2)\n",
+    "firedrake.assemble(abs(hs[-1] - expr) * dx) / firedrake.assemble(abs(hs[-1]) * dx)"
+   ]
   }
  ],
  "metadata": {

src/icepack2/model/__init__.py

@@ -11,18 +11,16 @@ def mass_balance(**kwargs):
     field_names = ("thickness", "velocity", "accumulation", "test_function")
     h, u, a, φ = map(kwargs.get, field_names)
     h_in = kwargs.get("thickness_inflow", Constant(0.0))
-    f_c = kwargs.get("frontal_ablation", Constant((0.0, 0.0)))
 
-    cell_balance = (Dt(h) * φ - inner(h * u - f_c, grad(φ)) - a * φ) * dx
+    cell_balance = (Dt(h) * φ - inner(h * u, grad(φ)) - a * φ) * dx
 
     mesh = ufl.domain.extract_unique_domain(h)
     ν = FacetNormal(mesh)
-    F = max_value(0, inner(h * u - f_c, ν))
+    f = h * max_value(0, inner(u, ν))
 
-    outflow = F * φ * ds
+    outflow = f * φ * ds
     inflow = h_in * min_value(0, inner(u, ν)) * φ * ds
-    boundary_balance = inflow + outflow
 
-    facet_balance = jump(F) * jump(φ) * dS
+    facet_balance = jump(f) * jump(φ) * dS
 
-    return cell_balance + facet_balance + boundary_balance
+    return cell_balance + facet_balance + inflow + outflow

test/mass_balance_test.py

@@ -2,7 +2,7 @@
 from numpy import pi as π
 import pytest
 import firedrake
-from firedrake import dx, inner, max_value, Constant
+from firedrake import dx, inner, max_value, Constant, conditional
 import irksome
 from icepack2 import model
 
@@ -49,19 +49,12 @@ def test_convergence_rate(degree):
         dt = firedrake.Constant(δt)
         final_time = 2 * π
 
-        lower, upper = firedrake.Function(Q), firedrake.Function(Q)
-        upper.assign(+np.inf)
-        bounds = ("stage", lower, upper)
         params = {
             "solver_parameters": {
-                "snes_type": "vinewtonrsls",
+                "snes_type": "ksponly",
                 "ksp_type": "gmres",
-                "pc_type": "lu",
-                "pc_factor_mat_solver_type": "mumps",
+                "pc_type": "bjacobi",
             },
-            "stage_type": "value",
-            "basis_type": "Bernstein",
-            "bounds": bounds,
         }
         solver = irksome.TimeStepper(problem, tableau, t, dt, h, **params)
 
@@ -80,3 +73,81 @@ def test_convergence_rate(degree):
     slope, intercept = np.polyfit(log_mesh_sizes, log_errors, 1)
     print(f"degree {degree}: log(error) ~= {slope:g} * log(dx) {intercept:+g}")
     assert slope > degree - 0.55
+
+
+@pytest.mark.parametrize("degree", [0, 1])
+def test_frontal_ablation(degree):
+    errors, mesh_sizes = [], []
+    k_min, k_max, num_steps = 5 - degree, 8 - degree, 9
+    for nx in np.logspace(k_min, k_max, num_steps, base=2, dtype=int):
+        mesh = firedrake.UnitSquareMesh(nx, nx, quadrilateral=True)
+        x = firedrake.SpatialCoordinate(mesh)
+
+        element = firedrake.FiniteElement("DQ", "quadrilateral", degree)
+        # Try Bernstein?
+
+        Q = firedrake.FunctionSpace(mesh, element)
+
+        u_max = Constant(1.0)
+        u = Constant((u_max, 0.0))
+
+        h_in = Constant(1.0)
+        x = firedrake.SpatialCoordinate(mesh)
+        L = Constant(0.25)
+        expr = conditional(x[0] < L, h_in, 0.0)
+        h_0 = firedrake.Function(Q).project(expr)
+        h = h_0.copy(deepcopy=True)
+
+        a = Constant(0.0)
+
+        U_c = Constant(0.5)
+        x_c = Constant(0.5)
+        u_c = firedrake.as_vector((conditional(x[0] >= x_c, U_c, 0), 0))
+
+        h_exact = conditional(x[0] <= x_c, h_in, h_in * u_max / (u_max + U_c))
+
+        ϕ = firedrake.TestFunction(Q)
+        problem = model.mass_balance(
+            thickness=h,
+            velocity=u + u_c,
+            accumulation=a,
+            thickness_inflow=h_in,
+            test_function=ϕ,
+        )
+
+        δx = mesh.cell_sizes.dat.data_ro.max()
+        δt = 0.5 * δx / u_max
+
+        t = Constant(0.0)
+        dt = Constant(δt)
+        final_time = 2.0
+
+        lower, upper = firedrake.Function(Q), firedrake.Function(Q)
+        upper.assign(np.inf)
+        params = {
+            "solver_parameters": {
+                "snes_type": "vinewtonrsls",
+                "snes_atol": 1e-12,
+            },
+            "stage_type": "value",
+            "basis_type": "Bernstein",
+            "bounds": ("stage", lower, upper),
+        }
+
+        tableau = irksome.BackwardEuler()
+        solver = irksome.TimeStepper(problem, tableau, t, dt, h, **params)
+
+        while float(t) < final_time:
+            if float(t) + float(dt) > final_time:
+                dt.assign(final_time - float(t))
+            solver.advance()
+            t.assign(t + dt)
+
+        mesh_sizes.append(δx)
+        errors.append(firedrake.assemble(abs(h - h_exact) * dx))
+        print(".", end="", flush=True)
+
+    log_mesh_sizes = np.log2(np.array(mesh_sizes))
+    log_errors = np.log2(np.array(errors))
+    slope, intercept = np.polyfit(log_mesh_sizes, log_errors, 1)
+    print(f"degree {degree}: log(error) ~= {slope:g} * log(dx) {intercept:+g}")