DO NOT MERGE

01-spd-helmholtz.ipynb

@@ -20,6 +20,7 @@
    "outputs": [],
    "source": [
     "# Code in this cell makes plots appear an appropriate size and resolution in the browser window\n",
+    "# If the following line fails then the user needs to install ipympl.\n",
     "%config InlineBackend.figure_format = 'svg'\n",
     "\n",
     "import matplotlib.pyplot as plt\n",
@@ -303,7 +304,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Since we know that the Helmholtz equation defines a symmetric problem, we instruct PETSc to employ the conjugate gradient method. We do not consider preconditioning, for now."
+    "We now solve the variational problem. Since we don't specify any solver parameters, the default direct solver will be used."
    ]
   },
   {
@@ -364,28 +365,6 @@
     "fig.colorbar(collection);"
    ]
   },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "# Exercises\n",
-    "\n",
-    "## Exercise 1: \n",
-    "\n",
-    "### 1a: use a higher approximation degree\n",
-    "\n",
-    "Solve the same problem, only this time, use a piecewise quadratic approximation space.\n",
-    "\n",
-    "- Hint: check the help for `FunctionSpace` to see how to specify the degree.\n",
-    "\n",
-    "### 1b: use a quadrilateral mesh\n",
-    "\n",
-    "Solve the same problem, but using quadrilateral, rather than triangular, cells.\n",
-    "\n",
-    "- Hint 1: check the help for `UnitSquareMesh` to see how to make a quadrilateral mesh\n",
-    "- Hint 2: To specify a piecewise continuous space on quadrilaterals, use the family name `\"Q\"`."
-   ]
-  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -427,7 +406,8 @@
     "L = f * v * dx\n",
     "\n",
     "uh = Function(V)\n",
-    "solve(a == L, uh)\n",
+    "solve(a == L, uh, solver_parameters={'ksp_type': 'cg',\n",
+    "                                     'pc_type': 'none'})\n",
     "\n",
     "u_exact = cos(2*pi*x)*cos(2*pi*y)"
    ]

02-poisson.ipynb

@@ -35,8 +35,10 @@
     "\n",
     "for some known function $f$. To have a well-posed problem, we must impose Dirichlet conditions over at least part of the domain boundary:\n",
     "\n",
-    "$$u(x) = g(x) \\quad \\forall x \\in \\Gamma_D,\\\\\n",
-    "\\nabla u(x)\\cdot \\vec{n} = h(x) \\quad \\forall x \\in \\Gamma_N.$$\n",
+    "$$\\begin{gather*}\n",
+    "u(x) = g(x) \\quad \\forall x \\in \\Gamma_D,\\\\\n",
+    "\\nabla u(x)\\cdot \\vec{n} = h(x) \\quad \\forall x \\in \\Gamma_N.\n",
+    "\\end{gather*}$$\n",
     "\n",
     "As before, the Neumann condition is imposed weakly by setting the boundary integral over the relevant part of the boundary. The Dirichlet condition is imposed strongly by modifying the function space from which $u$ is drawn.\n",
     "\n",

03-elasticity.ipynb

@@ -19,17 +19,17 @@
    "source": [
     "# Linear elasticity\n",
     "\n",
-    "*This work is adapted from an earlier FEniCS tutorial.*\n",
+    "*This work is adapted from a previous FEniCS tutorial no longer online*\n",
     "\n",
     "Having studied a few scalar-valued problems, we now move on to a vector-valued problem.  The equations of linear elasticity.  Here, we'll treat the isotropic case.\n",
     "\n",
     "For small deformations, the governing equation is:\n",
     "\n",
     "$$ -\\nabla \\cdot \\sigma = f \\text{ in } \\Omega, $$\n",
-    "with\n",
-    "$$ \\DeclareMathOperator{\\Tr}{Tr}\n",
-    "\\text{the stress tensor}\\quad \\sigma := \\lambda \\Tr(\\epsilon)\\mathbb{I} + 2\\mu\\epsilon\\\\\n",
-    "\\text{and the symmetric strain rate tensor}\\quad   \\epsilon := \\frac{1}{2}\\left(\\nabla u + (\\nabla u)^T\\right), $$\n",
+    "with the stress tensor:\n",
+    "$$ \\sigma := \\lambda \\Tr(\\epsilon)\\mathbb{I} + 2\\mu\\epsilon$$\n",
+    "and the symmetric strain rate tensor:\n",
+    "$$\\epsilon := \\frac{1}{2}\\left(\\nabla u + (\\nabla u)^T\\right), $$\n",
     "where $u$ is the unknown vector displacement field, and $\\mu$ and $\\lambda$ are the Lamè parameters.\n",
     "\n",
     "As before, the variational formulation consists of multiplying by a test function in some suitable finite element space, $v \\in V$, and integrating.  Note that this time, the solution $u$, and hence the test space $V$ are *vector*-valued (so multiplication actually means taking the inner product).\n",
@@ -110,7 +110,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "bc = DirichletBC(V, as_vector([0., 0.]), 1)"
+    "bc = DirichletBC(V, Constant([0, 0]), 1)"
    ]
   },
   {
@@ -264,7 +264,7 @@
     "    lambda_ = Constant(0.25)\n",
     "    Id = Identity(mesh.geometric_dimension()) # 2x2 Identity tensor\n",
     "    \n",
-    "    bc = DirichletBC(V, as_vector([0., 0.]), 1)\n",
+    "    bc = DirichletBC(V, Constant([0, 0]), 1)\n",
     "    u = TrialFunction(V)\n",
     "    v = TestFunction(V)\n",
     "    a = inner(sigma(u), epsilon(v))*dx\n",
@@ -345,7 +345,7 @@
     "\n",
     "    def sigma(u):\n",
     "        return lambda_*div(u)*Id + 2*mu*epsilon(u)    \n",
-    "    bc = DirichletBC(V, as_vector([0., 0.]), 1)\n",
+    "    bc = DirichletBC(V, Constant([0, 0]), 1)\n",
     "    u = TrialFunction(V)\n",
     "    v = TestFunction(V)\n",
     "    a = inner(sigma(u), epsilon(v))*dx\n",
@@ -356,8 +356,8 @@
     "    b0 = Function(V)\n",
     "    b1 = Function(V)\n",
     "    b2 = Function(V)\n",
-    "    b0.interpolate(as_vector([1., 0.]))\n",
-    "    b1.interpolate(as_vector([0., 1.]))\n",
+    "    b0.interpolate(Constant([1, 0]))\n",
+    "    b1.interpolate(Constant([0, 1]))\n",
     "    b2.interpolate(as_vector([-y, x]))\n",
     "    nullmodes = VectorSpaceBasis([b0, b1, b2])\n",
     "    # Make sure they're orthonormal.\n",

05-mixed-poisson.ipynb

@@ -115,7 +115,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "rt = FiniteElement(\"Raviart-Thomas\", triangle, 2, variant=\"integral\")\n",
+    "rt = FiniteElement(\"Raviart-Thomas\", triangle, 2)\n",
     "Sigma = FunctionSpace(mesh, rt)\n",
     "V = FunctionSpace(mesh, \"DG\", 1)"
    ]
@@ -206,7 +206,7 @@
     "\n",
     "tripcolor(uh, axes=axes[1])\n",
     "axes[1].set_aspect(\"equal\")\n",
-    "axes[1].set_title(\"$u$\")\n",
+    "axes[1].set_title(r\"$u$\")\n",
     "\n",
     "plt.tight_layout()"
    ]

06-pde-constrained-optimisation.ipynb

@@ -52,7 +52,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "from firedrake import *"
+    "from firedrake import *\n",
+    "from firedrake.__future__ import interpolate"
    ]
   },
   {
@@ -116,7 +117,7 @@
     "\\end{split}\n",
     "$$\n",
     "\n",
-    "Here, $u$ and $p$ are unknown velocity and pressure, $f$ is a prescribed inflow, $g$ is the control variable that we will optimise for and $\\alpha$ is a regularisation parameter.  This corresponds physically to minimising the loss of energy as heat by controlling the in/outflow on $\\Gamma_\\text{circ}$.  The regularisation parameter penalises too many non-zero control values."
+    "Here, $u$ and $p$ are unknown velocity and pressure, $f$ is a prescribed inflow, $g$ is the control variable that we will optimise for and $\\alpha$ is a regularisation parameter.  This corresponds physically to minimising the loss of energy as heat by controlling the in/outflow on $\\Gamma_\\text{circ}$.  The regularisation parameter penalises large control values."
    ]
   },
   {
@@ -440,6 +441,13 @@
     "\n",
     "How does it affect the solution before and after optimisation? "
    ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
   }
  ],
  "metadata": {

07-geometric-multigrid.ipynb

@@ -427,7 +427,7 @@
     "                       0)\n",
     "\n",
     "    value = as_vector([gbar*(1 - (12*t)**2), 0])\n",
-    "    bcs = [DirichletBC(W.sub(0), value, (1, 2)),\n",
+    "    bcs = [DirichletBC(W.sub(0), interpolate(value, V), (1, 2)),\n",
     "           DirichletBC(W.sub(0), zero(2), (3, 4))]\n",
     "    \n",
     "    a = (nu*inner(grad(u), grad(v)) - p*div(v) + q*div(u))*dx\n",

08-composable-solvers.ipynb

@@ -142,7 +142,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "nullspace = MixedVectorSpaceBasis(W, [W.sub(0), VectorSpaceBasis(constant=True)])"
+    "nullspace = MixedVectorSpaceBasis(W, [W.sub(0), VectorSpaceBasis(constant=True, comm=mesh.comm)])"
    ]
   },
   {
@@ -669,6 +669,13 @@
    "outputs": [],
    "source": []
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
   {
    "cell_type": "code",
    "execution_count": null,

09-hybridisation.ipynb

@@ -961,11 +961,32 @@
     "wh.assign(0.0)\n",
     "uD_solver_hybrid.solve()"
    ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
   }
  ],
  "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
   "language_info": {
-   "name": "python"
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.13.1"
   }
  },
  "nbformat": 4,

10-sum-factorisation.ipynb

@@ -406,8 +406,22 @@
   }
  ],
  "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
   "language_info": {
-   "name": "python"
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.13.1"
   }
  },
  "nbformat": 4,

11-extract-adjoint-solutions.ipynb

@@ -103,6 +103,7 @@
    "outputs": [],
    "source": [
     "from firedrake import *\n",
+    "from firedrake.__future__ import interpolate\n",
     "from firedrake.adjoint import *\n",
     "continue_annotation()"
    ]
@@ -122,7 +123,7 @@
    "source": [
     "# Define a simple mesh\n",
     "n = 32\n",
-    "mesh = UnitSquareMesh(n, n)\n",
+    "mesh = UnitSquareMesh(n, 1)\n",
     "\n",
     "# Define P2 function space and corresponding test function\n",
     "V = VectorFunctionSpace(mesh, \"Lagrange\", 2)\n",
@@ -134,7 +135,7 @@
     "\n",
     "# Assign initial condition\n",
     "x, y = SpatialCoordinate(mesh)\n",
-    "u_ = interpolate(as_vector([sin(pi*x), 0]), V)\n",
+    "u_ = assemble(interpolate(as_vector([sin(pi*x), 0]), V))\n",
     "u.assign(u_)\n",
     "\n",
     "# Set diffusivity constant\n",
@@ -246,7 +247,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "stop_annotating();"
+    "pause_annotation();"
    ]
   },
   {
@@ -330,14 +331,31 @@
     "For each solve block, `block`, the adjoint solution is stored as the attribute `block.adj_sol`."
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
   {
    "cell_type": "code",
    "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
     "for block in solve_blocks:\n",
-    "    assert block.adj_sol is not None"
+    "    assert block.adj_sol is not None\n",
+    "\n",
+    "def adj_solution(solve_block):\n",
+    "    return solve_block.get_dependencies()[0].adj_value"
    ]
   },
   {
@@ -364,7 +382,7 @@
     "for i in range(num_exports+1):\n",
     "    t = i*timesteps_per_export*dt\n",
     "    tricontourf(forward_solutions[i], axes=axs[i, 0])\n",
-    "    adjoint_solution = dJdu if i == num_exports else solve_blocks[timesteps_per_export*i].adj_sol\n",
+    "    adjoint_solution = dJdu if i == num_exports else adj_solution(solve_blocks[timesteps_per_export*i])\n",
     "    # Get the Riesz representer\n",
     "    adjoint_solution = dJdu.riesz_representation(riesz_map=\"H1\")\n",
     "    tricontourf(adjoint_solution, axes=axs[i, 1])\n",
@@ -378,8 +396,22 @@
   }
  ],
  "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
   "language_info": {
-   "name": "python"
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.13.1"
   }
  },
  "nbformat": 4,

12-HPC_demo.ipynb

@@ -3,7 +3,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "662e50fc",
+   "id": "6aaadb71",
    "metadata": {},
    "outputs": [],
    "source": [