Change exact solution in momentum balance test

Author: danshapero

test/momentum_balance_test.py

@@ -37,6 +37,7 @@
 def test_convergence_rate_grounded(degree, form):
     Lx, Ly = Constant(20e3), Constant(20e3)
     h0, dh = Constant(500.0), Constant(100.0)
+    s0, ds = Constant(150.0), Constant(90.0)
     u_inflow = Constant(100.0)
 
     n = firedrake.Constant(glen_flow_law)
@@ -46,28 +47,27 @@ def test_convergence_rate_grounded(degree, form):
     ε_c = Constant(0.01)
     A = ε_c / τ_c**n
 
-    height_above_flotation = Constant(10.0)
-    d = Constant(-ρ_I / ρ_W * (h0 - dh) + height_above_flotation)
-    ρ = Constant(ρ_I - ρ_W * d**2 / (h0 - dh) ** 2)
+    h_L = h0 - dh
+    s_L = s0 - ds
+    β = dh / ds * (ρ_I * h_L**2 - ρ_W * (s_L - h_L)**2) / (ρ_I * h_L**2)
 
     # We'll arbitrarily pick this to be the velocity, then we'll find a
     # friction coefficient and surface elevation that makes this velocity
     # an exact solution of the shelfy stream equations.
     def exact_u(x):
-        Z = A * (ρ * g * h0 / 4) ** n
-        q = 1 - (1 - (dh / h0) * (x / Lx)) ** (n + 1)
-        du = Z * q * Lx * (h0 / dh) / (n + 1)
+        ρ = β * ρ_I * ds / dh
+        h = h0 - dh * x / Lx
+        P = ρ * g * h / 4
+        dP = ρ * g * dh / 4
+        P0 = ρ * g * h0 / 4
+        du = Lx * A * (P0 ** (n + 1) - P ** (n + 1)) / ((n + 1) * dP)
         return u_inflow + du
 
 
-    # With this choice of friction coefficient, we can take the surface
-    # elevation to be a linear function of the horizontal coordinate and the
-    # velocity will be an exact solution of the shelfy stream equations.
-    β = Constant(0.5)
-    α = Constant(β * ρ / ρ_I * dh / Lx)
-
     def friction(x):
-        return α * (ρ_I * g * (h0 - dh * x / Lx)) * exact_u(x) ** (-1 / m)
+        h = h0 - dh * x / Lx
+        return (1 - β) * (ρ_I * g * h) * ds / Lx * exact_u(x) ** (-1 / m)
+
 
     errors, mesh_sizes = [], []
     k_min, k_max, num_steps = 5 - degree, 8 - degree, 9
@@ -88,8 +88,7 @@ def friction(x):
         u_exact = Function(V).interpolate(as_vector((exact_u(x), 0)))
 
         h = Function(Q).interpolate(h0 - dh * x / Lx)
-        ds = (1 + β) * ρ / ρ_I * dh
-        s = Function(Q).interpolate(d + h0 - dh + ds * (1 - x / Lx))
+        s = Function(Q).interpolate(s0 - ds * x / Lx)
 
         # TODO: adjust the yield stress so that this has a more sensible value
         C = Function(Q).interpolate(friction(x))