More on extrapolation

demos/larsen/extrapolate.py

@@ -11,7 +11,7 @@
 parser.add_argument("--input", default="larsen-initial.h5")
 parser.add_argument("--outline", default="larsen.geojson")
 parser.add_argument("--degree", type=int, default=1)
-parser.add_argument("--output")
+parser.add_argument("--output", default="larsen-extrapolated.h5")
 args = parser.parse_args()
 
 # Read in the input data
@@ -20,8 +20,8 @@
     u_input = chk.load_function(input_mesh, name="velocity")
     θ_input = chk.load_function(input_mesh, name="log_fluidity")
 
-    Q = θ_input.function_space()
-    μ = firedrake.interpolate(Constant(1), Q)
+    Δ = firedrake.FunctionSpace(input_mesh, "DG", 0)
+    μ = firedrake.interpolate(Constant(1), Δ)
 
 # Create the mesh and some function spaces
 with open(args.outline, "r") as outline_file:
@@ -36,42 +36,33 @@
 
 mesh = firedrake.Mesh("larsen.msh")
 
-import matplotlib.pyplot as plt
-fig, axes = plt.subplots()
-axes.set_aspect("equal")
-firedrake.triplot(mesh, axes=axes)
-axes.legend()
-plt.show()
-
 Q = firedrake.FunctionSpace(mesh, "CG", args.degree)
 V = firedrake.VectorFunctionSpace(mesh, "CG", args.degree)
 
 # Project the mask, log-fluidity, and velocity onto the larger mesh. The
 # regions with no data will be extrapolated by zero.
-μ = firedrake.project(μ, Q)
-μ.interpolate(firedrake.max_value(0, firedrake.min_value(1, μ)))
+Δ = firedrake.FunctionSpace(mesh, "DG", 0)
+μ = firedrake.project(μ, Δ)
 
 Eθ = firedrake.project(θ_input, Q)
 Eu = firedrake.project(u_input, V)
 
 θ = Eθ.copy(deepcopy=True)
 u = Eu.copy(deepcopy=True)
 
-bc = firedrake.DirichletBC(V, Eu, [4, 6, 7, 8, 9])
-
 α = Constant(8e3)
+
+bc = firedrake.DirichletBC(V, Eu, [4, 7, 8, 9])
 J = 0.5 * (μ * inner(u - Eu, u - Eu) + α**2 * inner(grad(u), grad(u))) * dx
 F = firedrake.derivative(J, u)
 firedrake.solve(F == 0, u, bc)
 
-area = assemble(μ * dx)
-print(np.sqrt(assemble(μ * inner(u - Eu, u - Eu) * dx) / area))
-
-import matplotlib.pyplot as plt
-fig, axes = plt.subplots()
-axes.set_aspect("equal")
-colors = firedrake.tripcolor(u, axes=axes)
-fig.colorbar(colors)
-plt.show()
-
+bc = firedrake.DirichletBC(Q, Eθ, [4, 7, 8, 9])
+J = 0.5 * (μ * inner(θ - Eθ, θ - Eθ) + α**2 * inner(grad(θ), grad(θ))) * dx
+F = firedrake.derivative(J, θ)
+firedrake.solve(F == 0, θ, bc)
 
+with firedrake.CheckpointFile(args.output, "w") as chk:
+    chk.save_mesh(mesh)
+    chk.save_function(u, name="velocity")
+    chk.save_function(θ, name="log_fluidity")