Make it go for nunataks

icepack · Oct 7, 2023 · 4169dc7 · 4169dc7
1 parent 69d5c48
commit 4169dc7
Show file tree

Hide file tree

Showing 2 changed files with 34 additions and 8 deletions.
diff --git 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
@@ -10,13 +10,16 @@
 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")
 args = parser.parse_args()
 
 # 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)