{ "cells": [ { "cell_type": "markdown", "id": "58eda45b-089d-4258-88fa-6ac36c24ed78", "metadata": {}, "source": [ "# Visualizing derivatives\n", "\n", "When we differentiate a functional with respect to a finite element discretized function in some discrete function space $V$, we do not obtain a function in $V$, but instead obtain an element of the associated dual space $V^*$. This notebook describes how a Riesz map can be used to construct a function, associated with the derivative, which can be visualized.\n", "\n", "This example makes use of tlm_adjoint with the [Firedrake](https://firedrakeproject.org/) backend, and we assume real spaces and a real build of Firedrake throughout.\n", "\n", "## Forward problem\n", "\n", "We consider the solution $u \\in V$ of a discretization of the Poisson equation subject to homogeneous Dirichlet boundary conditions,\n", "\n", "$$\\forall \\zeta \\in V_0 \\qquad \\int_\\Omega \\nabla \\zeta \\cdot \\nabla u = -\\int_\\Omega \\zeta m,$$\n", "\n", "where $V$ is a real $P_1$ continuous finite element space defining functions on the domain $\\Omega = \\left( 0, 1 \\right)^2$, with $m \\in V$, and where $V_0$ consists of the functions in $V$ which have zero trace. We define a functional\n", "\n", "$$J \\left( u \\right) = \\int_\\Omega \\left( 1 - x \\right)^4 u^2,$$\n", "\n", "where $x$ and $y$ denote Cartesian coordinates in $\\mathbb{R}^2$.\n", "\n", "We first solve the forward problem for\n", "\n", "$$m = \\mathcal{I} \\left[ \\sin \\left( \\pi x \\right) \\sin \\left( 2 \\pi y \\right) \\right],$$\n", "\n", "where $\\mathcal{I}$ maps to an element of $V$ through interpolation at mesh vertices." ] }, { "cell_type": "code", "execution_count": null, "id": "f977ffe4-03f5-426b-a437-aa012b4a125a", "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", "\n", "from firedrake import *\n", "from tlm_adjoint.firedrake import *\n", "from firedrake.pyplot import tricontourf\n", "\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "\n", "mesh = UnitSquareMesh(10, 10)\n", "X = SpatialCoordinate(mesh)\n", "space = FunctionSpace(mesh, \"Lagrange\", 1)\n", "test, trial = TestFunction(space), TrialFunction(space)\n", "m = Function(space, name=\"m\").interpolate(sin(pi * X[0]) * sin(2 * pi * X[1]))\n", "\n", "\n", "def forward(m):\n", " u = Function(space, name=\"u\")\n", " solve(inner(grad(trial), grad(test)) * dx == -inner(m, test) * dx,\n", " u, DirichletBC(space, 0.0, \"on_boundary\"))\n", "\n", " J = Functional(name=\"J\")\n", " J.assign(((1.0 - X[0]) ** 4) * u * u * dx)\n", " return u, J\n", "\n", "\n", "u, J = forward(m)\n", "\n", "print(f\"{J.value=}\")\n", "\n", "\n", "def plot_output(u, title):\n", " r = (u.dat.data_ro.min(), u.dat.data_ro.max())\n", " eps = (r[1] - r[0]) * 1.0e-12\n", " p = tricontourf(u, np.linspace(r[0] - eps, r[1] + eps, 32))\n", " plt.gca().set_title(title)\n", " plt.colorbar(p)\n", " plt.gca().set_aspect(1.0)\n", "\n", "\n", "plot_output(u, title=\"u\")" ] }, { "cell_type": "markdown", "id": "be99c89b-29c5-492e-b3bd-6792d3c3f4b2", "metadata": {}, "source": [ "## First order adjoint\n", "\n", "We can differentiate a functional respect to $m$ using `compute_gradient`. Specifically we compute the derivative of $\\hat{J} = J \\circ \\hat{u}$, where $\\hat{u}$ is the function which maps from the control $m$ to the solution to the discretized Poisson problem.\n", "\n", "When we compute a derivative of a functional with respect to a finite element discretized function using the adjoint method the result is not a function, but is instead a member of a dual space. Here we have $m \\in V$, and the derivative of $\\hat{J}$ with respect to $m$ is a member of the associated dual space $V^*$. In order to visualize this derivative we first need to map it to $V$. We can do this using a Riesz map. This is not uniquely defined, but here we choose to define a Riesz map using the $L^2$ inner product.\n", "\n", "Being precise the function we visualize, $g^\\sharp \\in V$, is defined such that\n", "\n", "$$\\forall \\zeta \\in V \\qquad \\left. \\frac{d \\hat{J} \\left( m + \\alpha \\zeta \\right)}{d \\alpha} \\right|_{\\alpha = 0} = \\int_\\Omega \\zeta g^\\sharp,$$\n", "\n", "where $\\alpha$ is a scalar." ] }, { "cell_type": "code", "execution_count": null, "id": "54f6e366-2cda-41f4-aa9c-d6c37d5bc26a", "metadata": {}, "outputs": [], "source": [ "reset_manager()\n", "\n", "start_manager()\n", "u, J = forward(m)\n", "stop_manager()\n", "\n", "print(f\"{J.value=}\")\n", "plot_output(u, title=\"u\")\n", "\n", "\n", "dJdm = compute_gradient(J, m)\n", "plot_output(dJdm.riesz_representation(\"L2\"), title=r\"$g^\\sharp$\")" ] }, { "cell_type": "markdown", "id": "16946cc4-64e1-438a-9f42-6a170decc5dc", "metadata": {}, "source": [ "## Hessian action\n", "\n", "Next we seek to differentiate $\\hat{J}$ twice with respect to $m$. However the second derivative defines a *bilinear* operator. This can be represented as a matrix – a Hessian matrix – but the number of elements in this matrix is equal to the *square* of the number of degrees of freedom for $m$.\n", "\n", "Instead of computing the full second derivative we can compute its action on a given direction $\\zeta \\in V$. The degrees of freedom associated with the result define the action of the Hessian matrix on a vector – specifically the action of the Hessian matrix on the vector consisting of the degrees of freedom for $\\zeta$.\n", "\n", "We do this in two stages. First we compute a directional derivative,\n", "\n", "$$\\left. \\frac{d \\hat{J} \\left( m + \\alpha \\zeta \\right)}{d \\alpha} \\right|_{\\alpha=0},$$\n", "\n", "computed using the tangent-linear method. We consider the case where $\\zeta = 1$." ] }, { "cell_type": "code", "execution_count": null, "id": "8e641946-4aeb-49e0-9415-6b86a822613e", "metadata": {}, "outputs": [], "source": [ "reset_manager()\n", "\n", "zeta = Function(space).interpolate(Constant(1.0))\n", "configure_tlm((m, zeta))\n", "\n", "start_manager()\n", "u, J = forward(m)\n", "stop_manager()\n", "\n", "print(f\"{J.value=}\")\n", "plot_output(u, title=\"u\")\n", "\n", "dJdm_zeta = var_tlm(J, (m, zeta))\n", "print(f\"{dJdm_zeta.value=}\")" ] }, { "cell_type": "markdown", "id": "74f20f32-c109-4506-8fbf-0db3e4cdffa1", "metadata": {}, "source": [ "Next we compute the derivative of *this* derivative using the adjoint method.\n", "\n", "Again we need to remember that the result is a member of the dual space $V^*$, and is not a function, so we again use a Riesz map to visualize it. Here we use the same Riesz map as before, defined using the $L^2$ inner product.\n", "\n", "Being precise the function we visualize, $h^\\sharp \\in V$, is defined such that\n", "\n", "$$\\forall \\chi \\in V \\qquad \\left. \\frac{\\partial^2 \\hat{J} \\left( m + \\alpha \\zeta + \\beta \\chi \\right)}{\\partial \\beta \\partial \\alpha} \\right|_{\\alpha = 0, \\beta = 0} = \\int_\\Omega \\chi h^\\sharp,$$\n", "\n", "where $\\zeta \\in V$ defines the direction on which the action is computed, and $\\alpha$ and $\\beta$ are scalars." ] }, { "cell_type": "code", "execution_count": null, "id": "6107d209-7731-42c3-8d84-f61b92d30232", "metadata": {}, "outputs": [], "source": [ "d2Jdm2_zeta = compute_gradient(dJdm_zeta, m)\n", "plot_output(d2Jdm2_zeta.riesz_representation(\"L2\"), title=r\"$h^\\sharp$\")" ] } ], "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.10.12" } }, "nbformat": 4, "nbformat_minor": 5 }