1d_periodic_Helmholtz.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "fS_kA1J9Ggg9"
   },
   "outputs": [],
   "source": [
    "try:\n",
    "    from firedrake import *\n",
    "except ImportError:\n",
    "    !wget \"https://fem-on-colab.github.io/releases/firedrake-install-real.sh\" -O \"/tmp/firedrake-install.sh\" && bash \"/tmp/firedrake-install.sh\"\n",
    "    from firedrake import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "PMzGFqKNJaLb"
   },
   "source": [
    "# Introduction to Firedrake using the 1D Helmholtz equation\n",
    "Firedrake is a library that provides tools for solving PDEs using the finite element method. We will see some of these tools in this notebook.\n",
    "\n",
    "The equation we will solve is the 1D Helmholtz equation\n",
    "\\begin{equation}\n",
    "-u_{xx}+ku = f,\n",
    "\\end{equation}\n",
    "where $k$ is a constant, $f$ is a given function and $u$ is the unknown function that we are solving for. We will solve this on the domain $0\\le x\\le 1$ and we will assume periodic boundary conditions. The first thing we need to do is set up a discretisation, or mesh, of this domain. We do this using one of the inbuilt Firedrake meshes. The input to this function is simply the number of cells in the domain, here 20:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "5NBbtrpSGgg-"
   },
   "outputs": [],
   "source": [
    "mesh = PeriodicUnitIntervalMesh(20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "ft-OREDuKwEb"
   },
   "source": [
    "We will now set up our right hand side function $f$ which is given by the expression\n",
    "\\begin{equation}\n",
    "f = (1+4\\pi^2)\\cos(2\\pi x).\n",
    "\\end{equation}\n",
    "This expression has been chosen because we know the analytic solution for $u$ in this case (via the \"method of manufactured solutions\"). Before we can set up a finite element function for $f$, we have to decide what finite element function space we are going to use. Here we will use a function space consisting of continous (\"CG\") piecewise linear functions (the \"1\" indicates the order of the polynomial - in this case linear)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "_khS6ydOKuim"
   },
   "outputs": [],
   "source": [
    "V = FunctionSpace(mesh, \"CG\", 1)\n",
    "f = Function(V)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "7ERJ9TzoLjxm"
   },
   "source": [
    "Now that we have a function, we can interpolate our chosen expression for $f$ and plot it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 282
    },
    "id": "_ozZE0Z-Ggg_",
    "outputId": "889d5b8d-27e8-4bc0-e83b-54c18f78d66f"
   },
   "outputs": [],
   "source": [
    "x = SpatialCoordinate(mesh)\n",
    "f.interpolate((1+4*pi**2)*cos(2*pi*x[0]))\n",
    "plot(f)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "0-euI4GsLz1o"
   },
   "source": [
    "Now we need to set up the weak form of our equation, which requires us to multiply by a test function $\\phi$ and integrate over the domain:\n",
    "\\begin{equation}\n",
    "\\int_0^1 (-\\phi u_{xx} + k\\phi u) \\, dx = \\int_0^1\\phi f \\, dx.\n",
    "\\end{equation}\n",
    "\n",
    "Recall that we chose to use a piecewise linear function space and note that we have a second order derivative on the left hand side. We can reduce the differentiability requirements on our solution $u$ by integrating the second order term by parts, giving:\n",
    "\\begin{equation}\n",
    "\\int_0^1 (\\phi_x u_x + k\\phi u) \\, dx = \\int_0^1\\phi f \\, dx.\n",
    "\\end{equation}\n",
    "\n",
    "We will take $k=1$ and set up our weak form by defining test and trial functions on our function space and writing $a$ for the left hand side and $L$ for the right hand side, we have:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "gk-jPiiAGghA"
   },
   "outputs": [],
   "source": [
    "phi = TestFunction(V)\n",
    "trial = TrialFunction(V)\n",
    "a = inner(grad(phi), grad(trial))*dx + inner(phi, trial)*dx\n",
    "L = inner(phi, f)*dx"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "CFCR-ap3QLLg"
   },
   "source": [
    "We need to set up a function to store the solution and then we can solve the linear system and plot the solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 282
    },
    "id": "Q-zpzRQCQHFp",
    "outputId": "ed02436b-f29b-43e6-a74a-5268de051b09"
   },
   "outputs": [],
   "source": [
    "u = Function(V)\n",
    "solve(a==L, u)\n",
    "plot(u)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "-H5iA8-2QYIb"
   },
   "source": [
    "We can also set up a function and interpolate the known analytical solution to it. We can plot this and also compute an error norm to see how well our approximation has worked."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 282
    },
    "id": "A5TE0aSWGghA",
    "outputId": "835caf28-b19b-4da1-b5e6-8984fa8788af"
   },
   "outputs": [],
   "source": [
    "u_exact = Function(V).interpolate(cos(2*pi*x[0]))\n",
    "plot(u_exact)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "iVQ6w1INGghB",
    "outputId": "c2a98898-e1fd-4411-e6f9-41979dd943d2"
   },
   "outputs": [],
   "source": [
    "print(errornorm(u, u_exact))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Y_f3z7tMTVhV"
   },
   "source": [
    "## Exercises:\n",
    "* Increase the resolution - try doubling the number of cells.\n",
    "* Increase the order of the function space to use piecewise quadratics by changing the line defining $V$ and passing in 2 instead of 1.\n",
    "* Solve this equation in 2D by changing the mesh to be PeriodicUnitSquareMesh(20, 20). If the analytical solution is\n",
    "\\begin{equation}\n",
    "u=\\cos(2\\pi x)\\cos(2\\pi y),\n",
    "\\end{equation}\n",
    "what does $f$ have to be for the two dimensional Helmholtz equation\n",
    "\\begin{equation}\n",
    "-\\nabla^2 u + u = f\n",
    "\\end{equation}\n",
    "to hold?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "bhhZ6mjYGghB"
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "colab": {
   "provenance": []
  },
  "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.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}