From f24037aedc18dde62a8af3c78539dc64f5b65059 Mon Sep 17 00:00:00 2001
From: Daniel Shapero <shapero.daniel@gmail.com>
Date: Fri, 25 Aug 2023 11:53:47 -0700
Subject: [PATCH] Save results to file

---
 demos/mismip+/mismip.py | 52 +++++++++++++++++++----------------------
 1 file changed, 24 insertions(+), 28 deletions(-)

diff --git a/demos/mismip+/mismip.py b/demos/mismip+/mismip.py
index b4b9735..16d7a56 100644
--- a/demos/mismip+/mismip.py
+++ b/demos/mismip+/mismip.py
@@ -2,7 +2,7 @@
 import tqdm
 import firedrake
 from firedrake import (
-    sqrt, exp, max_value, inner, grad, div, as_vector, Constant, dx, ds, interpolate
+    sqrt, exp, max_value, min_value, inner, grad, div, Constant, dx, ds
 )
 import icepack
 from dualform import ice_stream
@@ -26,7 +26,7 @@
 Q = firedrake.FunctionSpace(mesh, cg1)
 V = firedrake.VectorFunctionSpace(mesh, cg1)
 Σ = firedrake.TensorFunctionSpace(mesh, dg0, symmetry=True)
-T = firedrake.VectorFunctionSpace(mesh, cg1)
+T = firedrake.VectorFunctionSpace(mesh, dg0)
 Z = V * Σ * T
 
 # Set up the bed topography
@@ -51,7 +51,7 @@
 )
 
 z_deep = Constant(-720)
-b = interpolate(max_value(B_x + B_y, z_deep), Q)
+b = firedrake.interpolate(max_value(B_x + B_y, z_deep), Q)
 
 # Choose the yield stress, strain rate, and velocity to give the specified
 # fluidity value of A = 20 MPa^{-3} yr^{-1} and slipperiness of
@@ -61,7 +61,7 @@
 ε_c = Constant(0.02)    # 1 / yr
 u_c = Constant(1000.0)  # m / yr
 
-h_0 = interpolate(Constant(100), Q)
+h_0 = firedrake.interpolate(Constant(100), Q)
 s_0 = icepack.compute_surface(thickness=h_0, bed=b)
 
 # Set up the variational forms for the problem and a linearization
@@ -69,7 +69,7 @@
 z = firedrake.Function(Z)
 u, M, τ = firedrake.split(z)
 
-u_0 = interpolate(as_vector((90 * x / Lx, 0)), V)
+u_0 = firedrake.interpolate(firedrake.as_vector((90 * x / Lx, 0)), V)
 z.sub(0).assign(u_0)
 fields = {
     "velocity": u,
@@ -130,7 +130,7 @@
 bcs = [bc_inflow, bc_sides]
 
 # We're starting from nothing, so do an initial solve for a velocity field from
-# a linearized equation set to get us halfway there
+# a linearized equation set to get us halfway there.
 fparams = {"form_compiler_parameters": {"quadrature_degree": 6}}
 
 lparams = {
@@ -147,7 +147,7 @@
 lsolver.solve()
 
 # Now we'll use a continuation method in the flow law and sliding parameters to
-# get to the real starting velocity field
+# get to the real starting velocity field.
 nparams = {
     "solver_parameters": {
         "snes_type": "newtonls",
@@ -167,44 +167,40 @@
     n.assign(t * icepack.constants.glen_flow_law + (1 - t) * n_0)
     nsolver.solve()
 
-# Form the mass conservation equation and solver
+# Form the mass conservation equation and solver; we use a Lax-Wendroff or 2nd-
+# order Taylor discretization in time. This adds some diffusion along-flow,
+# which increases the order and also helps with stability.
 a = firedrake.Constant(0.3)
-dt = firedrake.Constant(0.125)
+dt = firedrake.Constant(1.0)
 ν = firedrake.FacetNormal(mesh)
 h_k = h.copy(deepcopy=True)
 φ = firedrake.TestFunction(Q)
 G = (
     ((h - h_k) / dt * φ - (inner(h * u, grad(φ)) + a * φ)) * dx +
     0.5 * dt * div(h * u) * inner(u, grad(φ)) * dx +
-    h * firedrake.max_value(0, inner(u, ν)) * φ * ds -
-    0.5 * dt * div(h * u) * firedrake.max_value(0, inner(u, ν)) * φ * ds -
-    0.5 * dt * div(h_0 * u) * firedrake.min_value(0, inner(u, ν)) * φ * ds +
-    h_0 * firedrake.min_value(0, inner(u, ν)) * φ * ds
+    h * max_value(0, inner(u, ν)) * φ * ds -
+    0.5 * dt * div(h * u) * max_value(0, inner(u, ν)) * φ * ds -
+    0.5 * dt * div(h_0 * u) * min_value(0, inner(u, ν)) * φ * ds +
+    h_0 * min_value(0, inner(u, ν)) * φ * ds
 )
 
 hproblem = firedrake.NonlinearVariationalProblem(G, h)
 hsolver = firedrake.NonlinearVariationalSolver(hproblem, **lparams)
 
 # Run the simulation
-final_time = 120
+final_time = 10000
 num_steps = int(final_time / float(dt))
 for step in tqdm.trange(num_steps):
     hsolver.solve()
-    h.interpolate(firedrake.max_value(0, h))
+    h.interpolate(max_value(0, h))
     h_k.assign(h)
 
     nsolver.solve()
 
-import matplotlib.pyplot as plt
-fig, axes = plt.subplots()
-axes.set_aspect("equal")
-colors = firedrake.tripcolor(h, axes=axes)
-fig.colorbar(colors)
-plt.show()
-
-u, M, τ = z.subfunctions
-fig, axes = plt.subplots()
-axes.set_aspect("equal")
-colors = firedrake.tripcolor(u, axes=axes)
-fig.colorbar(colors)
-plt.show()
+# Save the results to disk
+with firedrake.CheckpointFile("mismip.h5", "w") as chk:
+    u, M, τ = z.subfunctions
+    chk.save_function(h, name="thickness")
+    chk.save_function(u, name="velocity")
+    chk.save_function(M, name="membrane_stress")
+    chk.save_function(τ, name="basal_stress")