Start with basic notebook to test frontal ablation

Author: danshapero

notebooks/frontal-ablation.ipynb

@@ -0,0 +1,194 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a13bc8d5-3247-4792-931b-cca7574d1e34",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "from matplotlib.animation import FuncAnimation\n",
+    "from IPython.display import HTML\n",
+    "from tqdm.notebook import trange, tqdm\n",
+    "import firedrake\n",
+    "from firedrake import (\n",
+    "    Constant, inner, min_value, max_value, jump, avg, grad, dx, ds, dS\n",
+    ")\n",
+    "import irksome\n",
+    "from irksome import Dt"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "0b47295b-abd4-43cf-9e4c-618646526567",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "nx, ny = 32, 32\n",
+    "mesh = firedrake.UnitSquareMesh(nx, ny, diagonal=\"crossed\")\n",
+    "δ = 1.0 / nx"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "b9ac67fe-69cf-4cd4-8b3c-76b1bb7bcf19",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "degree = 0\n",
+    "element = firedrake.FiniteElement(\"DG\", \"triangle\", degree)\n",
+    "Q = firedrake.FunctionSpace(mesh, element)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "1afae3de-d7b9-49a1-9868-a4768f72d164",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "u_max = 1.0\n",
+    "u = firedrake.Constant((u_max, 0.0))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "36de950f-a4bf-4e0e-9b47-cae2a23fbd51",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "h_in = Constant(1.0)\n",
+    "x = firedrake.SpatialCoordinate(mesh)\n",
+    "L = Constant(0.25)\n",
+    "expr = firedrake.conditional(x[0] < L, 1.0, 0.0)\n",
+    "h_0 = firedrake.Function(Q).project(expr)\n",
+    "h = h_0.copy(deepcopy=True)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "de1cf5af-901c-4a6b-b537-9c070d337edd",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ϕ = firedrake.TestFunction(Q)\n",
+    "F_cells = (Dt(h) * ϕ - inner(h * u, grad(ϕ))) * dx\n",
+    "\n",
+    "ν = firedrake.FacetNormal(mesh)\n",
+    "f = max_value(0, h * inner(u, ν))\n",
+    "F_facets = jump(f) * jump(ϕ) * dS\n",
+    "\n",
+    "F_inflow = h_in * ϕ * min_value(0, inner(u, ν)) * ds\n",
+    "F_outflow = h * ϕ * max_value(0, inner(u, ν)) * ds\n",
+    "\n",
+    "F = F_cells + F_facets + F_inflow + F_outflow"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "9a458b42-4fd6-4e55-b0ce-83ebb0e8b99f",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "method = irksome.BackwardEuler()\n",
+    "t = Constant(0.0)\n",
+    "dt = Constant(0.5 * δ / u_max)\n",
+    "\n",
+    "params = {\n",
+    "    \"solver_parameters\": {\n",
+    "        \"snes_type\": \"newtonls\",\n",
+    "        \"ksp_type\": \"gmres\",\n",
+    "        \"pc_type\": \"lu\",\n",
+    "        \"pc_factor_mat_solver_type\": \"mumps\",\n",
+    "    },\n",
+    "}\n",
+    "\n",
+    "solver = irksome.TimeStepper(F, method, t, dt, h, **params)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "c34b4498-b788-4c07-a3f0-5c3a1aa58dba",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "hs = [h.copy(deepcopy=True)]\n",
+    "\n",
+    "final_time = 1.0\n",
+    "num_steps = int(final_time / float(dt))\n",
+    "for step in trange(num_steps):\n",
+    "    solver.advance()\n",
+    "    hs.append(h.copy(deepcopy=True))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a3a9e4fb-fa7d-4afc-ac4b-06d4e93c66d6",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%%capture\n",
+    "\n",
+    "fig, ax = plt.subplots()\n",
+    "ax.set_aspect(\"equal\")\n",
+    "colors = firedrake.tripcolor(\n",
+    "    hs[0], axes=ax, num_sample_points=1, shading=\"gouraud\"\n",
+    ");\n",
+    "\n",
+    "fn_plotter = firedrake.FunctionPlotter(mesh, num_sample_points=1)\n",
+    "def animate(h):\n",
+    "    colors.set_array(fn_plotter(h))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "1816a106-b3e7-4112-ac60-f073ea8d9120",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "animation = FuncAnimation(fig, animate, tqdm(hs), interval=1e3/10)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "446a0a61-0164-4367-b425-437b8f1bcabe",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "HTML(animation.to_html5_video())"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "firedrake",
+   "language": "python",
+   "name": "firedrake"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.11.9"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}