From 80ca89a4f406b63f1faf92b095c6ff2a87535d80 Mon Sep 17 00:00:00 2001 From: David Ham Date: Fri, 4 Apr 2025 12:04:01 +0100 Subject: [PATCH 1/6] further updated notebooks --- 01-spd-helmholtz.ipynb | 28 +----- 02-poisson.ipynb | 6 +- 03-elasticity.ipynb | 20 ++-- 05-mixed-poisson.ipynb | 4 +- 06-pde-constrained-optimisation.ipynb | 16 ++-- 07-geometric-multigrid.ipynb | 2 +- 08-composable-solvers.ipynb | 9 +- 09-hybridisation.ipynb | 23 ++++- 10-sum-factorisation.ipynb | 16 +++- 11-extract-adjoint-solutions.ipynb | 132 ++++++++++++++++++++++++-- 12-HPC_demo.ipynb | 2 +- 11 files changed, 203 insertions(+), 55 deletions(-) diff --git a/01-spd-helmholtz.ipynb b/01-spd-helmholtz.ipynb index 6cc2779..50ebb21 100644 --- a/01-spd-helmholtz.ipynb +++ b/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)" ] diff --git a/02-poisson.ipynb b/02-poisson.ipynb index e6be072..607794d 100644 --- a/02-poisson.ipynb +++ b/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", diff --git a/03-elasticity.ipynb b/03-elasticity.ipynb index dd9d859..68f6091 100644 --- a/03-elasticity.ipynb +++ b/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", diff --git a/05-mixed-poisson.ipynb b/05-mixed-poisson.ipynb index 2739243..1115d7d 100644 --- a/05-mixed-poisson.ipynb +++ b/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()" ] diff --git a/06-pde-constrained-optimisation.ipynb b/06-pde-constrained-optimisation.ipynb index b970b75..f345922 100644 --- a/06-pde-constrained-optimisation.ipynb +++ b/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." ] }, { @@ -132,10 +133,6 @@ "metadata": {}, "outputs": [], "source": [ - "import os\n", - "if not os.path.isfile(\"stokes-control.msh\"):\n", - " # If the mesh is not available locally, we download it.\n", - " !curl -O https://raw.githubusercontent.com/firedrakeproject/notebooks/refs/heads/main/stokes-control.msh\n", "mesh = Mesh(\"stokes-control.msh\")" ] }, @@ -440,6 +437,13 @@ "\n", "How does it affect the solution before and after optimisation? " ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { diff --git a/07-geometric-multigrid.ipynb b/07-geometric-multigrid.ipynb index caf5cc6..70075a3 100644 --- a/07-geometric-multigrid.ipynb +++ b/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", diff --git a/08-composable-solvers.ipynb b/08-composable-solvers.ipynb index e50b6ae..3f07f29 100644 --- a/08-composable-solvers.ipynb +++ b/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, diff --git a/09-hybridisation.ipynb b/09-hybridisation.ipynb index 8e51d0f..a7f2052 100644 --- a/09-hybridisation.ipynb +++ b/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, diff --git a/10-sum-factorisation.ipynb b/10-sum-factorisation.ipynb index c157df9..c2c666d 100644 --- a/10-sum-factorisation.ipynb +++ b/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, diff --git a/11-extract-adjoint-solutions.ipynb b/11-extract-adjoint-solutions.ipynb index 0c72c58..a17812c 100644 --- a/11-extract-adjoint-solutions.ipynb +++ b/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,6 +331,20 @@ "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, @@ -337,7 +352,10 @@ "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", @@ -375,11 +393,113 @@ "axs[0, 0].set_title('Forward solution');\n", "axs[0, 1].set_title('Adjoint solution');" ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "type(solve_blocks[6].get_dependencies()[0].adj_value)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "[type(sol.get_dependencies()[0].adj_value) for sol in solve_blocks]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "len(forward_solutions)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "get_working_tape().visualise(\"foo.pdf\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "adjoint_solution" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "Jhat = ReducedFunctional(J, Control(nu))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "Jhat.derivative()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "taylor_test?" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "taylor_test(Jhat, nu, Function(nu).assign(1.0))" + ] + }, + { + "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, diff --git a/12-HPC_demo.ipynb b/12-HPC_demo.ipynb index fa5ed95..d53a650 100644 --- a/12-HPC_demo.ipynb +++ b/12-HPC_demo.ipynb @@ -3,7 +3,7 @@ { "cell_type": "code", "execution_count": null, - "id": "662e50fc", + "id": "9aa45217", "metadata": {}, "outputs": [], "source": [ From f6830f4dc7e7033d713b2040f089c9077890ab06 Mon Sep 17 00:00:00 2001 From: David Ham Date: Fri, 4 Apr 2025 13:19:33 +0100 Subject: [PATCH 2/6] fixes --- 06-pde-constrained-optimisation.ipynb | 4 ++ 11-extract-adjoint-solutions.ipynb | 88 --------------------------- 12-HPC_demo.ipynb | 2 +- 3 files changed, 5 insertions(+), 89 deletions(-) diff --git a/06-pde-constrained-optimisation.ipynb b/06-pde-constrained-optimisation.ipynb index f345922..fc42bb7 100644 --- a/06-pde-constrained-optimisation.ipynb +++ b/06-pde-constrained-optimisation.ipynb @@ -133,6 +133,10 @@ "metadata": {}, "outputs": [], "source": [ + "import os\n", + "if not os.path.isfile(\"stokes-control.msh\"):\n", + " # If the mesh is not available locally, we download it.\n", + " !curl -O https://raw.githubusercontent.com/firedrakeproject/notebooks/refs/heads/main/stokes-control.msh\n", "mesh = Mesh(\"stokes-control.msh\")" ] }, diff --git a/11-extract-adjoint-solutions.ipynb b/11-extract-adjoint-solutions.ipynb index a17812c..3e64384 100644 --- a/11-extract-adjoint-solutions.ipynb +++ b/11-extract-adjoint-solutions.ipynb @@ -393,94 +393,6 @@ "axs[0, 0].set_title('Forward solution');\n", "axs[0, 1].set_title('Adjoint solution');" ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "type(solve_blocks[6].get_dependencies()[0].adj_value)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "[type(sol.get_dependencies()[0].adj_value) for sol in solve_blocks]" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "len(forward_solutions)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "get_working_tape().visualise(\"foo.pdf\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "adjoint_solution" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "Jhat = ReducedFunctional(J, Control(nu))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "Jhat.derivative()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "taylor_test?" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "taylor_test(Jhat, nu, Function(nu).assign(1.0))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] } ], "metadata": { diff --git a/12-HPC_demo.ipynb b/12-HPC_demo.ipynb index d53a650..2289ea4 100644 --- a/12-HPC_demo.ipynb +++ b/12-HPC_demo.ipynb @@ -3,7 +3,7 @@ { "cell_type": "code", "execution_count": null, - "id": "9aa45217", + "id": "6aaadb71", "metadata": {}, "outputs": [], "source": [ From 11ce5a2fb3e3a1f0e030c8f9e485df13067aac9a Mon Sep 17 00:00:00 2001 From: David Ham Date: Sun, 6 Apr 2025 18:44:46 +0100 Subject: [PATCH 3/6] new adjoint notebook --- 13_adjoint_calculations.ipynb | 629 ++++++++++++++++++++++++++++++++++ 1 file changed, 629 insertions(+) create mode 100644 13_adjoint_calculations.ipynb diff --git a/13_adjoint_calculations.ipynb b/13_adjoint_calculations.ipynb new file mode 100644 index 0000000..f877d01 --- /dev/null +++ b/13_adjoint_calculations.ipynb @@ -0,0 +1,629 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "c40a1722", + "metadata": {}, + "outputs": [], + "source": [ + "try:\n", + " import firedrake\n", + "except ImportError:\n", + " !wget \"https://fem-on-colab.github.io/releases/firedrake-install-release-real.sh\" -O \"/tmp/firedrake-install.sh\" && bash \"/tmp/firedrake-install.sh\"\n", + " import firedrake" + ] + }, + { + "cell_type": "markdown", + "id": "e7b5808f-f0ae-4e40-bcca-eff67f251fc9", + "metadata": {}, + "source": [ + "# Solving adjoint problems\n", + "\n", + "Firedrake has built-in capabilities for differentiating the solution of a PDE with respect to any of its inputs. These are documented in [the manual](https://www.firedrakeproject.org/adjoint.html). Here we will explore these capabilities using a the same Burgers equation example that we already explored when we learned how to solve nonlinear problems. As before, we will define the equation in one spatial dimension and avoid considering boundary conditions by selecting a periodic domain:\n", + "\n", + "$$\n", + "\\frac{\\partial u}{\\partial t} + u \\frac{\\partial u}{\\partial x} - \\nu \\frac{\\partial^2 u}{\\partial x^2} = 0.\n", + "$$\n", + "\n", + "Applying finite element in time and implicit Euler in space, the problem becomes at each timestep find $u^{n+1}\\in V$ such that:\n", + "\n", + "$$\n", + "\\int_\\Omega \\frac{u^{n+1} - u^n}{\\Delta t} v + u^{n+1} \\frac{\\partial u^{n+1}}{\\partial x} v + \\nu \\frac{\\partial u^{n+1}}{\\partial x}\\frac{\\partial v}{\\partial x} \\, \\mathrm{d}x = 0 \\quad \\forall v \\in V.\n", + "$$\n", + "\n", + "As usual, we start by setting up plotting:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "11709a56-6feb-4cac-9a98-3327e8c65a56", + "metadata": {}, + "outputs": [], + "source": [ + "# Code in this cell makes plots appear an appropriate size and resolution in the browser window\n", + "%config InlineBackend.figure_format = 'svg'\n", + "\n", + "import matplotlib.pyplot as plt\n", + "\n", + "plt.rcParams['figure.figsize'] = (11, 6)" + ] + }, + { + "cell_type": "markdown", + "id": "e46e290c-1b65-42c5-8f73-686331d28fdc", + "metadata": {}, + "source": [ + "... and importing Firedrake:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "efa08a55-d8f3-422c-9e52-789fcacd2cbe", + "metadata": {}, + "outputs": [], + "source": [ + "from firedrake import *\n", + "from firedrake.__future__ import interpolate" + ] + }, + { + "cell_type": "markdown", + "id": "8047d692-358e-4f52-b678-505ec3e65533", + "metadata": {}, + "source": [ + "This is where the code is slightly different because we're solving an adjoint problem. Firedrake uses [pyadjoint](https://pyadjoint.org) to deliver its adjoint capability, but the relevant parts of the Pyadjoint interface are exposed via the `firedrake.adjoint` module." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "24e4fc04-bc5a-4f57-a80d-da3fb6be1677", + "metadata": {}, + "outputs": [], + "source": [ + "from firedrake.adjoint import *" + ] + }, + { + "cell_type": "markdown", + "id": "a9cd6fc6-1d30-4eb7-8232-9b0014dd7b2a", + "metadata": {}, + "source": [ + "Executing the adjoint problem depends on having recorded all of the operations on which the results of the PDE depend. We switch on the recording thus:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ae90fb73-f2dd-4026-880d-59b37f086c73", + "metadata": {}, + "outputs": [], + "source": [ + "continue_annotation()" + ] + }, + { + "cell_type": "markdown", + "id": "2725a8f3-c7b1-480e-92d2-71a605f40d84", + "metadata": {}, + "source": [ + "We next define the mesh, initial conditions and residual for the Burgers equation as we have done in previous notebooks." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d8e381b1-8d8e-4960-b04e-3b5f6edee6ba", + "metadata": {}, + "outputs": [], + "source": [ + "n = 100\n", + "with stop_annotating():\n", + " mesh = PeriodicIntervalMesh(n, length=2)\n", + "x = SpatialCoordinate(mesh)[0]\n", + "nu = Constant(1e-2)\n", + "V = FunctionSpace(mesh, \"Lagrange\", 2)\n", + "u_n1 = Function(V, name=\"u^{n+1}\")\n", + "u_n = Function(V, name=\"u^{n}\")\n", + "v = TestFunction(V)\n", + "u_init = Function(V).interpolate(sin(2*pi*x))\n", + "u_n.assign(u_init)\n", + "dt = 1.0 / n\n", + "F = (((u_n1 - u_n)/dt) * v +\n", + " u_n1 * u_n1.dx(0) * v + \n", + " nu*u_n1.dx(0)*v.dx(0))*dx" + ] + }, + { + "cell_type": "markdown", + "id": "1f9cb135-762c-428f-ad43-78c5d80fa348", + "metadata": {}, + "source": [ + "Now we timestep the equations as usual:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "725c7b79-29d5-482e-b606-761f7ec04f9f", + "metadata": {}, + "outputs": [], + "source": [ + "problem = NonlinearVariationalProblem(F, u_n1)\n", + "solver = NonlinearVariationalSolver(problem)\n", + "results = [Function(u_n)]\n", + "timesteps = 50\n", + "for t in ProgressBar(\"Time step\").iter(range(timesteps)):\n", + " solver.solve()\n", + " u_n.assign(u_n1)\n", + " results.append(Function(u_n))" + ] + }, + { + "cell_type": "markdown", + "id": "37638133-5468-421c-b663-7b2f7d7128cd", + "metadata": {}, + "source": [ + "The basic form of adjoint computiation computes the derivative of a scalar outcome of a PDE with respect to one or more inputs. Here we'd just like to see how the solution of the PDE at the end of the time interval depends on the initial conditions, so we choose the squared $L^2$ norm of the final solution as our functional:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ee93684f-cb5e-4867-b36b-22084c6ca73e", + "metadata": {}, + "outputs": [], + "source": [ + "J = assemble(u_n*u_n*dx)" + ] + }, + { + "cell_type": "markdown", + "id": "f894db4d-af9f-429c-a3d7-c317a2ba78dc", + "metadata": {}, + "source": [ + "We've now completed the forward solution. We want to stop recording so that any computations we do on the results aren't part of the tape:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "dbed6910-f7d6-4cc0-8e03-0dc1a4bfd068", + "metadata": {}, + "outputs": [], + "source": [ + "pause_annotation()" + ] + }, + { + "cell_type": "markdown", + "id": "bc7cdd0d-3df2-4265-9568-e22cdc039ea5", + "metadata": {}, + "source": [ + "Just as with the forward calculation, it's helpful to have a visualisation of the execution of the model. However, the adjoint calculations will be executed automatically, so we can't directly attach a progress bar. Instead, we pass a progress bar class to the tape and let it apply it automatically:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3ea75d66-9c55-4f22-9db2-7a13c8e41148", + "metadata": {}, + "outputs": [], + "source": [ + "get_working_tape().progress_bar = ProgressBar\n" + ] + }, + { + "cell_type": "markdown", + "id": "7f67d11c-5a84-4cb5-a365-d98a3e849027", + "metadata": {}, + "source": [ + "We select the initial conditions of the simulation as the control. The way we achieve this in code is by defining a [`ReducedFunctional`](https://www.firedrakeproject.org/adjoint.html#reduced-functionals) with the initial conditions as the [`Control`](https://pyadjoint.org/documentation/pyadjoint_api.html#pyadjoint.Control):" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7807c438-6175-4789-9624-52307943011e", + "metadata": {}, + "outputs": [], + "source": [ + "Jhat = ReducedFunctional(J, Control(u_init))" + ] + }, + { + "cell_type": "markdown", + "id": "229b7676-d483-47fa-9a45-0f7aa3b1ebd0", + "metadata": {}, + "source": [ + "The reduced functional is the core object in any adjoint calculation. For example, we can compute the derivative of a ReducedFunctional as follows. Notice the progress bar. It has 154 entries in it in contrast to the 50 timesteps. This is because the adjoint calculation records and replays individual Firedrake operations such as `solve` and `assign` rather than whole timesteps." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "43496003-814b-49d4-91e8-ce38a335e60f", + "metadata": {}, + "outputs": [], + "source": [ + "dJ = Jhat.derivative()" + ] + }, + { + "cell_type": "markdown", + "id": "00b6b619-14fa-47e7-96b4-f7cce2de1559", + "metadata": {}, + "source": [ + "Let's look at the derivative we've calculated. Can you intuitively work out why this is the right derivative?" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cdc905b6-48bc-43e6-89c5-18bac8ea7e1d", + "metadata": {}, + "outputs": [], + "source": [ + "# NBVAL_IGNORE_OUTPUT\n", + "from firedrake.pyplot import plot\n", + "\n", + "fig, axes = plt.subplots()\n", + "\n", + "paths=plot(dJ, axes=axes)" + ] + }, + { + "cell_type": "markdown", + "id": "3bedff1a-68b6-4bf1-b4d1-d7da00e6be07", + "metadata": {}, + "source": [ + "## Exercise 1\n", + "\n", + "The reduced functional defines a function:\n", + "\n", + "$$\n", + "\\hat{J}(m) = J(u(m), m)\n", + "$$\n", + "\n", + "Where $m$ is the control. As well as differentiating the functional, it's possible to evaluate it for a new intial condition. This is achieved by calling the reduced functional:\n", + "\n", + "```python3\n", + "Jhat(u_new)\n", + "```\n", + "\n", + "### a.\n", + "Try evaluating `Jhat` for a new initial condition. You'll still need to use a periodic funcition, but you could change the frequency of the wave, or try a more adventurous but still periodic function.\n", + "\n", + "### b. \n", + "You can access the solution of the last solve in the reduced functional with `get_solve_blocks()[-1].get_outputs()[0].saved_output`. See how the solution changes for different values of the initial condition.\n", + "\n", + "### c.\n", + "The gradient returned by `Jhat.derivative()` is linearised around the solution for the last state at which the reduced functional was evaluated. Plot the gradient for different initial conditions and see if you can intuit why it looks like it does.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "24821568-b9b3-4c26-a8c9-3dc948a42cc1", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "679e0876-e9fc-4af5-9154-bd931dea324d", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "id": "48a49a11-626e-457e-9c85-7d8824566b52", + "metadata": {}, + "source": [ + "## Taylor tests\n", + "\n", + "How do we know that the gradient we have calculated above is correct? The [Taylor test](https://www.firedrakeproject.org/adjoint.html#taylor-tests) provides us a very sensitive test of this. Almost any error in the computation of the gradient will cause its convergence rate to drop below 2. We evaluate the Taylor test at our original functional and with a perturbation of a constant function:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f20c41ff-de35-4cc0-b8cd-3598b4d7ebf1", + "metadata": {}, + "outputs": [], + "source": [ + "taylor_test(Jhat, u_init, assemble(interpolate(Constant(1.), V))) " + ] + }, + { + "cell_type": "markdown", + "id": "f27e3478-3ef4-41af-bb3a-ae9ae58f05eb", + "metadata": {}, + "source": [ + "## Examining the adjoint calculation\n", + "\n", + "When we first studied the Burgers equation we plotted an animation of the solution over time. As a reminder:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ffc4579b-f1d0-45bf-972c-7ad9ff2ad0f9", + "metadata": {}, + "outputs": [], + "source": [ + "# NBVAL_IGNORE_OUTPUT\n", + "from firedrake.pyplot import plot\n", + "from matplotlib.animation import FuncAnimation\n", + "\n", + "fig, axes = plt.subplots()\n", + "\n", + "def animate(u):\n", + " axes.clear()\n", + " plot(u, axes=axes)\n", + " axes.set_ylim((-1., 1.))\n", + "\n", + "interval = 4e3 * float(dt)\n", + "animation = FuncAnimation(fig, animate, frames=results, interval=interval)\n", + "\n", + "from IPython.display import HTML\n", + "HTML(animation.to_jshtml())" + ] + }, + { + "cell_type": "markdown", + "id": "1b7e37e9-5b11-4ae7-b5cd-e2cd9b4ba69a", + "metadata": {}, + "source": [ + "In order to understand how the adjoint solution works, it's helpful to do the same for the adjoint. We can extract the adjoint solutions to the PDE for each timestep as follows:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "23fe145a-df41-4362-b24f-54199baa8299", + "metadata": {}, + "outputs": [], + "source": [ + "adjoints = [block.get_dependencies()[0].adj_value for block in get_solve_blocks()[::-1]]" + ] + }, + { + "cell_type": "markdown", + "id": "748827e1-4b50-4411-bab1-c32179051faf", + "metadata": {}, + "source": [ + "We deliberately extracted the blocks in reverse order because the adjoint is solved in reverse order. The adjoint value of the control is the last solve value in the sequence so we add that:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a2cfd7de-6956-488f-a387-c83af6cb54a8", + "metadata": {}, + "outputs": [], + "source": [ + "adjoints += [Jhat.controls[0].adj_value]" + ] + }, + { + "cell_type": "markdown", + "id": "8e265b61-cf41-4014-8a56-ca99eac81328", + "metadata": {}, + "source": [ + "Now, if we take a look at those variables..." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2a092582-a161-45f7-8ff5-e334e67ec80f", + "metadata": {}, + "outputs": [], + "source": [ + "adjoints" + ] + }, + { + "cell_type": "markdown", + "id": "f14c44b4-810c-4ac9-9f57-c1358172fad8", + "metadata": {}, + "source": [ + "The adjoint variables are, as expected, cofunctions. To recover the adjoint solutions in a form we can visualise, we take the Riesz represntation to recover Functions:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c76e7526-d1b8-4040-ba0a-3a7b4673a5f1", + "metadata": {}, + "outputs": [], + "source": [ + "adjoints = [a.riesz_representation(riesz_map=\"L2\") for a in adjoints if a] " + ] + }, + { + "cell_type": "markdown", + "id": "5327cf8c-fe20-4022-a42f-d5bb0b6cac9e", + "metadata": {}, + "source": [ + "Now we can visualise the adjoint solution:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0f3723f8-8df1-4095-a6bc-f2bd1589380e", + "metadata": {}, + "outputs": [], + "source": [ + "fig, axes = plt.subplots()\n", + "\n", + "def animate(u):\n", + " axes.clear()\n", + " plot(u, axes=axes)\n", + " axes.set_ylim((-2.5, 2.5))\n", + "\n", + "interval = 4e3 * float(dt)\n", + "animation = FuncAnimation(fig, animate, frames=adjoints, interval=interval)\n", + "\n", + "from IPython.display import HTML\n", + "HTML(animation.to_jshtml())" + ] + }, + { + "cell_type": "markdown", + "id": "405540f6-367e-4f2d-be39-46b88d35cd35", + "metadata": {}, + "source": [ + "## Exercise 2\n", + "\n", + "Go back to the Poisson equation example, and solve Poisson's equation with Dirichlet boundary conditions. Find $u$ such that:\n", + "$$\n", + "\\begin{gather*}\n", + "-\\nabla^2 u = f \\qquad x\\in \\Omega\\\\\n", + "u = 0 \\qquad x \\in \\Gamma\n", + "\\end{gather*}\n", + "$$\n", + "Choose a suitable forcing function $f$ and evaluate the functional:\n", + "$$\n", + "J = \\int_\\Omega u^2 + f^2\\,\\mathrm{d}x\n", + "$$\n", + "Create a reduced functional and conduct a Taylor test in order to establish that your implementation is correct.\n", + "Before you start, you will need to remove the recording of the Burgers equation from the tape:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "17eabff0-9f62-4861-a701-a059626e251b", + "metadata": {}, + "outputs": [], + "source": [ + "get_working_tape().clear_tape()\n", + "continue_annotation()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6ff6d0fa-49ef-4ccc-989c-071da4d3c451", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "50169141-e31c-4b1d-a951-1bb2c89567ef", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "49ab1bc0-9706-4759-af26-8ee8c24f2ada", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "id": "ce79c50c-133e-4014-9d13-8e95bbf7fa9d", + "metadata": {}, + "source": [ + "## Visualising the tape\n", + "\n", + "An exceptionally useful aid to understanding the computation of the adjoint is to produce a visual representation of the tape. To produce a managable sized figure, we'll rerun the Burgers equation with just 5 timesteps:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9c1d602e-d5f0-4f34-a08f-8cb4fa959a02", + "metadata": {}, + "outputs": [], + "source": [ + "tape=get_working_tape()\n", + "continue_annotation()\n", + "tape.clear_tape()\n", + "u_n.assign(u_init)\n", + "timesteps = 5\n", + "for t in ProgressBar(\"Time step\").iter(range(timesteps)):\n", + " solver.solve()\n", + " u_n.assign(u_n1)\n", + "J = assemble(u_n*u_n*dx)\n", + "pause_annotation()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "31926022-2d1a-4726-a12a-0ea78f7a0538", + "metadata": {}, + "outputs": [], + "source": [ + "tape.visualise(\"tape.pdf\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "78e43809-5941-4ffd-a52d-c0ebe6dfb465", + "metadata": {}, + "outputs": [], + "source": [ + "!pip install pdf2image" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a722ed53-7931-478f-a718-41c6a61f0810", + "metadata": {}, + "outputs": [], + "source": [ + "from pdf2image import convert_from_path\n", + "\n", + "images = convert_from_path(\"tape.pdf\")\n", + "images[0] # first page" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7f581a32-8a82-4588-aa4a-fb1f689ab381", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "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.13.1" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} From a4ab531d8a2045b7f35170ae16f907d8a6f0247e Mon Sep 17 00:00:00 2001 From: David Ham Date: Thu, 10 Apr 2025 08:00:25 +0200 Subject: [PATCH 4/6] package installs --- 13_adjoint_calculations.ipynb | 37 +++++++++++++++++++++++------------ 1 file changed, 25 insertions(+), 12 deletions(-) diff --git a/13_adjoint_calculations.ipynb b/13_adjoint_calculations.ipynb index f877d01..ea4cccf 100644 --- a/13_adjoint_calculations.ipynb +++ b/13_adjoint_calculations.ipynb @@ -3,7 +3,7 @@ { "cell_type": "code", "execution_count": null, - "id": "c40a1722", + "id": "5165fa2d", "metadata": {}, "outputs": [], "source": [ @@ -526,14 +526,6 @@ "outputs": [], "source": [] }, - { - "cell_type": "code", - "execution_count": null, - "id": "49ab1bc0-9706-4759-af26-8ee8c24f2ada", - "metadata": {}, - "outputs": [], - "source": [] - }, { "cell_type": "markdown", "id": "ce79c50c-133e-4014-9d13-8e95bbf7fa9d", @@ -563,24 +555,45 @@ "pause_annotation()" ] }, + { + "cell_type": "markdown", + "id": "d0aafb5e-3f62-4b39-bc51-a9aabca49c01", + "metadata": {}, + "source": [ + "We now need to install a few packages that Colab don't install automatically:" + ] + }, { "cell_type": "code", "execution_count": null, - "id": "31926022-2d1a-4726-a12a-0ea78f7a0538", + "id": "7c1833ed-c09b-41a4-aff5-94c44989d1d8", "metadata": {}, "outputs": [], "source": [ - "tape.visualise(\"tape.pdf\")" + "!apt update\n", + "!apt install graphviz libgraphviz-dev poppler-utils" ] }, { "cell_type": "code", "execution_count": null, "id": "78e43809-5941-4ffd-a52d-c0ebe6dfb465", + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "!pip install graphviz pdf2image" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "31926022-2d1a-4726-a12a-0ea78f7a0538", "metadata": {}, "outputs": [], "source": [ - "!pip install pdf2image" + "tape.visualise(\"tape.pdf\")" ] }, { From f896d3ae17f56f37daaca412d0992cee8b5414eb Mon Sep 17 00:00:00 2001 From: David Ham Date: Thu, 10 Apr 2025 08:11:09 +0200 Subject: [PATCH 5/6] correct package name --- 13_adjoint_calculations.ipynb | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/13_adjoint_calculations.ipynb b/13_adjoint_calculations.ipynb index ea4cccf..d622023 100644 --- a/13_adjoint_calculations.ipynb +++ b/13_adjoint_calculations.ipynb @@ -3,7 +3,7 @@ { "cell_type": "code", "execution_count": null, - "id": "5165fa2d", + "id": "719efb41", "metadata": {}, "outputs": [], "source": [ @@ -583,7 +583,7 @@ }, "outputs": [], "source": [ - "!pip install graphviz pdf2image" + "!pip install pygraphviz pdf2image" ] }, { From e51c130bff9a8549eecfb53307ca414ee448bbef Mon Sep 17 00:00:00 2001 From: David Ham Date: Fri, 11 Apr 2025 08:03:36 +0200 Subject: [PATCH 6/6] extra package --- 13_adjoint_calculations.ipynb | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/13_adjoint_calculations.ipynb b/13_adjoint_calculations.ipynb index d622023..6d6011e 100644 --- a/13_adjoint_calculations.ipynb +++ b/13_adjoint_calculations.ipynb @@ -3,7 +3,7 @@ { "cell_type": "code", "execution_count": null, - "id": "719efb41", + "id": "c1993e37", "metadata": {}, "outputs": [], "source": [ @@ -571,7 +571,7 @@ "outputs": [], "source": [ "!apt update\n", - "!apt install graphviz libgraphviz-dev poppler-utils" + "!apt install graphviz libgraphviz-dev poppler-utils fonts-dejavu" ] }, { @@ -615,6 +615,16 @@ "id": "7f581a32-8a82-4588-aa4a-fb1f689ab381", "metadata": {}, "outputs": [], + "source": [ + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5ac64346-bd72-4a1c-ac0c-21ec13cbf4ae", + "metadata": {}, + "outputs": [], "source": [] } ],