Put bed topo in separate file

demos/mismip/bed_topo.py

@@ -0,0 +1,24 @@
+from firedrake import exp, max_value, Constant
+
+
+def mismip_bed_topography(x, y, Lx, Ly):
+    x_c = Constant(300e3)
+    X = x / x_c
+
+    B_0 = Constant(-150)
+    B_2 = Constant(-728.8)
+    B_4 = Constant(343.91)
+    B_6 = Constant(-50.57)
+    B_x = B_0 + B_2 * X**2 + B_4 * X**4 + B_6 * X**6
+
+    f_c = Constant(4e3)
+    d_c = Constant(500)
+    w_c = Constant(24e3)
+
+    B_y = d_c * (
+        1 / (1 + exp(-2 * (y - Ly / 2 - w_c) / f_c)) +
+        1 / (1 + exp(+2 * (y - Ly / 2 + w_c) / f_c))
+    )
+
+    z_deep = Constant(-720)
+    return max_value(B_x + B_y, z_deep)

demos/mismip/mismip.py

@@ -4,8 +4,6 @@
 import numpy as np
 import firedrake
 from firedrake import (
-    sqrt,
-    exp,
     min_value,
     max_value,
     as_vector,
@@ -25,6 +23,7 @@
     weertman_sliding_law as m,
 )
 from icepack2 import model
+import bed_topo
 
 
 # Get the command-line options.
@@ -57,27 +56,8 @@
 
 # Set up the basal elevation.
 x, y = firedrake.SpatialCoordinate(mesh)
-
-x_c = Constant(300e3)
-X = x / x_c
-
-B_0 = Constant(-150)
-B_2 = Constant(-728.8)
-B_4 = Constant(343.91)
-B_6 = Constant(-50.57)
-B_x = B_0 + B_2 * X**2 + B_4 * X**4 + B_6 * X**6
-
-f_c = Constant(4e3)
-d_c = Constant(500)
-w_c = Constant(24e3)
-
-B_y = d_c * (
-    1 / (1 + exp(-2 * (y - Ly / 2 - w_c) / f_c)) +
-    1 / (1 + exp(+2 * (y - Ly / 2 + w_c) / f_c))
-)
-
-z_deep = Constant(-720)
-b = interpolate(max_value(B_x + B_y, z_deep), S)
+expr = bed_topo.mismip_bed_topography(x, y, Lx, Ly)
+b = interpolate(expr, S)
 
 # Some physical constants; `A = ε_c / τ_c ** n`, `C = τ_c / u_c ** (1 / m)`.
 # Rather than work with the fluidity `A` and friction `C` directly, we use