Test if manufactured solution is a solution

Author: danshapero

test/composite_rheology_test.py

@@ -1,4 +1,5 @@
 import numpy as np
+import sympy
 import firedrake
 from firedrake import (
     as_vector,
@@ -122,3 +123,48 @@ def perturb_u(x, y):
     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.9
+
+
+def test_manufactured_solution():
+    L, ρ_I, ρ_W, g = sympy.symbols("L ρ_I ρ_W g", real=True, positive=True)
+    A_1, A_3 = sympy.symbols("A_1 A_3", real=True, positive=True)
+    K_1, K_3 = sympy.symbols("K_1 K_3", real=True, positive=True)
+    n, m = sympy.symbols("n m", integer=True, positive=True)
+
+    def strain_rate(M, A_1, A_3):
+        return 2 * A_1 * M + 2 * A_3 * M ** 3
+
+    def sliding_velocity(τ, K_1, K_3):
+        return K_1 * τ + K_3 * τ**3
+
+    def stress_balance(x, h, s, M, τ):
+        return diff(h * M, x) + τ - ρ_I * g * h * diff(s, x)
+
+    def boundary_condition(x, u, h, s, M):
+        d = (s - h).subs(x, L)
+        τ_c = (ρ_I * g * h**2 - ρ_W * g * d**2) / 2
+        return (h * M - τ_c).subs(x, L)
+
+    x = sympy.symbols("x", real=True)
+    h_0, δh = sympy.symbols("h_0 δh", real=True, positive=True)
+    h = h_0 - δh * x / L
+
+    # How the hell did I pick this again
+    h_f = sympy.symbols("h_f", real=True, positive=True)
+    d = -ρ_I / ρ_W * h.subs(x, L) + h_f
+    ρ = (ρ_I - ρ_W * d**2 / h**2).subs(x, L)
+
+    u_0 = sympy.symbols("u_0", real=True, positive=True)
+    P = ρ * g * h / 4
+    P_0 = ρ * g * h_0 / 4
+    δP = ρ * g * δh / 4
+    du_1 = L * A * (P_0 ** 2 - P ** 2) / (2 * δP)
+    du_3 = L * A * (P_0 ** (n + 1) - P ** (n + 1)) / ((n + 1) * δP)
+    u = u_0 + du_1 + du_3
+
+    β = 1 / 2
+    α = β * ρ / ρ_I * δh / L
+    # Miracle occurs??
+
+    ds = (1 + β) * ρ / ρ_I * δh
+    s = d + h.subs(x, L) + ds * (1 - x / L)