Make it work for n = 3

We need to prime the pump a bit by calling the time stepper first with n
= 1 and continue up to 3 but otherwise it works and with a longer
timestep than with clamping.

Author: danshapero

notebooks/rainier-monolithic.ipynb

@@ -248,10 +248,7 @@
     "ρ_I = Constant(ice_density)\n",
     "g = Constant(gravity)\n",
     "\n",
-    "\n",
-    "#n = Constant(glen_flow_law)\n",
-    "n = Constant(1.0)\n",
-    "\n",
+    "n = Constant(glen_flow_law)\n",
     "\n",
     "P = ρ_I * g * h\n",
     "S_n = inner(grad(s), grad(s))**((n - 1) / 2)\n",
@@ -412,7 +409,7 @@
     "\n",
     "sparams = {\n",
     "    \"solver_parameters\": {\n",
-    "        \"snes_monitor\": \":rainier-output-monolothic.txt\",\n",
+    "        \"snes_monitor\": \":rainier-output-initial.log\",\n",
     "        \"snes_type\": \"newtonls\",\n",
     "        \"snes_max_it\": 200,\n",
     "        \"snes_linesearch_type\": \"nleqerr\",\n",
@@ -462,13 +459,10 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "\"\"\"\n",
     "num_continuation_steps = 5\n",
     "for exponent in np.linspace(1.0, 3.0, num_continuation_steps):\n",
     "    n.assign(exponent)\n",
-    "    momentum_solver.solve()\n",
-    "\"\"\"\n",
-    "momentum_solver.solve()"
+    "    momentum_solver.solve()"
    ]
   },
   {
@@ -586,7 +580,7 @@
    "source": [
     "tableau = irksome.BackwardEuler()\n",
     "t = Constant(0.0)\n",
-    "dt = Constant(1.0 / 12)\n",
+    "dt = Constant(1.0 / 6)\n",
     "\n",
     "lower = firedrake.Function(W)\n",
     "upper = firedrake.Function(W)\n",
@@ -597,7 +591,7 @@
     "\n",
     "bparams = {\n",
     "    \"solver_parameters\": {\n",
-    "        \"snes_vi_monitor\": None,\n",
+    "        \"snes_monitor\": \":rainier-output-vi.log\",\n",
     "        \"snes_type\": \"vinewtonrsls\",\n",
     "        \"snes_max_it\": 200,\n",
     "        #\"snes_linesearch_type\": \"nleqerr\",\n",
@@ -621,6 +615,39 @@
     "The good part: run the model for 500 years of simulation time."
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "ffbd91d3-b0e4-416d-822d-2dd8f5f78a89",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "n.assign(1.0)\n",
+    "solver.advance(bounds=bounds)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "f360f81b-5c4f-4615-8937-b8f316099a0d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "n.assign(2.0)\n",
+    "solver.advance(bounds=bounds)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "3b9f978b-d4ca-4c74-82c8-cf159e527c50",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "n.assign(glen_flow_law)\n",
+    "solver.advance(bounds=bounds)"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -633,11 +660,9 @@
     "\n",
     "final_time = 500.0\n",
     "num_steps = int(final_time / float(dt))\n",
-    "progress_bar = trange(num_steps)\n",
-    "for step in progress_bar:\n",
+    "for step in trange(num_steps):\n",
     "    solver.advance(bounds=bounds)\n",
     "    h = w.subfunctions[3]\n",
-    "    progress_bar.set_description(str(h.dat.data_ro.min()))\n",
     "    s.interpolate(b + h)\n",
     "    a.interpolate(smb(s))\n",
     "\n",
@@ -678,7 +703,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "u, M, τ = z.subfunctions\n",
+    "u, M, τ, h = w.subfunctions\n",
     "\n",
     "fig, axes = plt.subplots()\n",
     "axes.set_aspect(\"equal\")\n",
@@ -725,6 +750,18 @@
     "from IPython.display import HTML\n",
     "HTML(animation.to_html5_video())"
    ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "53927034-2a55-4391-b0af-ac84f555f26a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "volumes = [firedrake.assemble(h * dx) for h in hs]\n",
+    "fig, ax = plt.subplots()\n",
+    "ax.plot(volumes);"
+   ]
   }
  ],
  "metadata": {