From 4169dc71b76df501c1c6e07127c121d20f9985f2 Mon Sep 17 00:00:00 2001
From: Daniel Shapero <shapero.daniel@gmail.com>
Date: Sat, 7 Oct 2023 15:55:43 -0700
Subject: [PATCH] Make it go for nunataks

---
 demos/helheim/helheim.ipynb | 39 +++++++++++++++++++++++++++++--------
 demos/helheim/preprocess.py |  3 +++
 2 files changed, 34 insertions(+), 8 deletions(-)

diff --git a/demos/helheim/helheim.ipynb b/demos/helheim/helheim.ipynb
index acec93f..00e61b2 100644
--- a/demos/helheim/helheim.ipynb
+++ b/demos/helheim/helheim.ipynb
@@ -64,7 +64,7 @@
    "source": [
     "fig, ax = plt.subplots()\n",
     "ax.set_aspect(\"equal\")\n",
-    "firedrake.triplot(mesh, axes=ax)\n",
+    "firedrake.triplot(mesh, boundary_kw={\"colors\": (\"tab:red\", \"tab:green\", \"tab:blue\", \"tab:orange\")}, axes=ax)\n",
     "ax.legend();"
    ]
   },
@@ -191,8 +191,29 @@
     "τ_v = Constant(0.1)\n",
     "\n",
     "ϕ = Constant(0.9)\n",
-    "τ_f = firedrake.interpolate(ϕ * sqrt(inner(τ_d, τ_d)), Q)\n",
-    "u_f = firedrake.interpolate(sqrt(inner(u, u)), Q)"
+    "τ_min = Constant(0.01)\n",
+    "τ_f = firedrake.interpolate(max_value(τ_min, ϕ * sqrt(inner(τ_d, τ_d))), Q)\n",
+    "u_f = firedrake.interpolate(sqrt(inner(u, u)), Q)\n",
+    "\n",
+    "h_c = firedrake.Constant(50.0)\n",
+    "K_c = firedrake.Constant(3e3)\n",
+    "K = firedrake.interpolate(firedrake.conditional(h < h_c, K_c, u_f / τ_f ** m), Q)\n",
+    "\n",
+    "u_f = firedrake.Constant(1000.0)\n",
+    "τ_f = firedrake.interpolate((u_f / K) ** (1 / m), Q)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a161f4d4-0264-4952-bbf1-4a977d55145a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, ax = plt.subplots()\n",
+    "ax.set_aspect(\"equal\")\n",
+    "colors = firedrake.tripcolor(τ_f, axes=ax)\n",
+    "fig.colorbar(colors);"
    ]
   },
   {
@@ -202,7 +223,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "z.sub(2).interpolate(τ_d / 2);"
+    "z.sub(2).interpolate(ϕ * τ_d);"
    ]
   },
   {
@@ -212,7 +233,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "inflow_ids = (1, 2, 4, 5, 6)\n",
+    "inflow_ids = (1, 2, 4)\n",
     "outflow_ids = (3,)\n",
     "\n",
     "bc = firedrake.DirichletBC(Z.sub(0), u, inflow_ids)"
@@ -307,11 +328,13 @@
     "\n",
     "sparams = {\n",
     "    \"solver_parameters\": {\n",
-    "        \"snes_max_it\": 50,\n",
+    "        \"snes_max_it\": 200,\n",
     "        \"snes_monitor\": None,\n",
+    "        #\"snes_linesearch_monitor\": None,\n",
     "        \"snes_rtol\": 5e-7,\n",
-    "        \"snes_type\": \"newtonls\",\n",
-    "        \"snes_linesearch_type\": \"l2\",\n",
+    "        \"snes_type\": \"newtontr\",\n",
+    "        #\"snes_linesearch_type\": \"bt\",\n",
+    "        #\"snes_linesearch_max_it\": 8,\n",
     "        \"ksp_type\": \"gmres\",\n",
     "        \"pc_type\": \"lu\",\n",
     "        \"pc_factor_mat_solver_type\": \"mumps\",\n",
diff --git a/demos/helheim/preprocess.py b/demos/helheim/preprocess.py
index 8d6d7da..6607fc9 100644
--- a/demos/helheim/preprocess.py
+++ b/demos/helheim/preprocess.py
@@ -10,6 +10,7 @@
 import data
 
 parser = argparse.ArgumentParser()
+parser.add_argument("--outline")
 parser.add_argument("--num-levels", type=int, default=1)
 parser.add_argument("--regularization", type=float, default=80.0)
 parser.add_argument("--output")
@@ -17,6 +18,8 @@
 
 # Fetch the glacier outline, generate a mesh, and create a function space
 outline_filename = icepack.datasets.fetch_outline("helheim")
+if args.outline:
+    outline_filename = args.outline
 with open(outline_filename, "r") as outline_file:
     outline = geojson.load(outline_file)