12-HPC_demo.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "662e50fc",
   "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": "403628da",
   "metadata": {},
   "source": [
    "# HPC Demonstration\n",
    "\n",
    "In this notebook we build up a multigrid solver for an elliptic problem specifically designed for running Firedrake on a High Performance Computer (HPC). We will solve very large instances of the Poisson problem, demonstrating a range of different solver options and assessing their performance for a range of problem sizes. Additional supplimentary material is provided for running scripts on HPC.\n",
    "\n",
    "**Note:** The code in this notebook is designed to be run on a large HPC facility, as a result some cells may take a long time to run in an interactive notebook. We suggest _not_ re-running the notebook cells, but instead trying the exercises on a supercomputer.\n",
    "\n",
    "We start as always by importing Firedrake. We also define parprint to perform parallel printing as in this [demonstration](https://firedrakeproject.org/demos/parprint.py.html). The Python time module is imported to benchmark the different solvers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4c87c6a1",
   "metadata": {},
   "outputs": [],
   "source": [
    "from firedrake import *\n",
    "from firedrake.petsc import PETSc\n",
    "from time import time\n",
    "\n",
    "parprint = PETSc.Sys.Print"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "58bef8e4",
   "metadata": {},
   "source": [
    "## How big?\n",
    "\n",
    "The parameters `Nx`, `Nref` and `degree` defined below have been selected so that the simulation runs on a single core in a notebook. This is not the regime we want to think about in this tutorial, we want to think about very large problems. We will consider how each of these parameters affects the overall problem size."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "909b129e",
   "metadata": {},
   "outputs": [],
   "source": [
    "Nx = 8\n",
    "Nref = 2\n",
    "degree = 2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "370155eb",
   "metadata": {},
   "source": [
    "These three parameters determine the total number of degrees of freedom (DOFs) in our problem:\n",
    "- `Nx` defines our coarse grid in the mesh hierarchy, it is used to construct a coarse cube mesh.\n",
    "- `Nref` determines how many times the mesh is refined to create a mesh hierarchy.\n",
    "- `degree`, which we denote $k$, specifies the polynomial order of the basis functions used to approximate functions in our finite element space.\n",
    "\n",
    "The total number of DOFs is given by:\n",
    "$$\n",
    "n = (k \\times N_x \\times 2^{N_{ref}} + 1)^d\n",
    "$$\n",
    "where $d=3$ is the dimension of the domain in which we solve the problem.\n",
    "\n",
    "This small notebook example solves a problem with a large number of DOFs, but on HPC we want to solve problems _orders of magnitude larger still_, by the end of this notebook we will be considering problems larger than 30 000 000 DOFs.\n",
    "\n",
    "When solving problems using Firedrake in parallel, it's worth remembering that performance can be improved by adding more processes (MPI ranks) as long as the number of DOFs remains above [50 000 DOFs per core](https://firedrakeproject.org/parallelism.html#expected-performance-improvements)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "111ba871",
   "metadata": {},
   "source": [
    "## The equations\n",
    "We will consider the Poisson equation in a 3D domain $\\Omega = [0, 1]^3$:\n",
    "\n",
    "$$\n",
    "\\left\\{\n",
    "\\begin{aligned}\n",
    "\t-\\nabla^2 u &= f && \\text{on } \\Omega,\\\\\n",
    "\tu &= 0 && \\text{on } \\partial\\Omega,\n",
    "\\end{aligned}\n",
    "\\right.\n",
    "$$\n",
    "\n",
    "where $f$ is given by:\n",
    "\n",
    "$$\n",
    "f(x,y,z) = -\\frac{\\pi^2}{2}\n",
    "\\times\\left( 2\\cos(\\pi x) - \\cos\\left( \\frac{\\pi x}{2} \\right)\n",
    "- 2(a^2 + b^2)\\sin(\\pi x)\\tan \\left( \\frac{\\pi x}{4} \\right)  \\right)\n",
    "\\times\\sin(a\\pi y) \\sin(b\\pi z).\n",
    "$$\n",
    "\n",
    "We use this particular right hand side since it has corresponding analytic solution:\n",
    "\n",
    "$$\n",
    "u(x,y,z) =\n",
    "\\sin(\\pi x)\\tan\\left(\\frac{\\pi x}{4}\\right)\n",
    "\\sin(a\\pi y)\\sin(b\\pi z).\n",
    "$$\n",
    "Having an analytic solution allows us to compute the error in our computed solution as $e_h = \\|u_h - u\\|_{L^2}$. For this notebook we fix $a=1$ and $b=2$.\n",
    "\n",
    "The Poisson equation has the weak form: Find $u_h \\in V$ such that\n",
    "\n",
    "$$\n",
    "\\int_\\Omega \\nabla u_h\\cdot \\nabla v\\ dx = \\int_\\Omega f v\\ dx \\qquad \\forall v \\in V.\n",
    "$$\n",
    "\n",
    "For the discrete function space $V$ we initially consider piecewise quadratic Lagrange elements, that is \n",
    "```python\n",
    "V = FunctionSpace(mesh, \"CG\", 2)\n",
    "```\n",
    "\n",
    "It is straightforward to solve the equation using Firedrake by expressing this weak form in UFL.\n",
    "The Python code below generates a `problem` object of the desired size, a function `u_h` to store the solution and the analytic solution `truth` so we can compute the $L_2$ error norm, all of which we use throughout the rest of the notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "540607f9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create mesh and mesh hierarchy\n",
    "mesh = UnitCubeMesh(Nx, Nx, Nx)\n",
    "hierarchy = MeshHierarchy(mesh, Nref)\n",
    "mesh = hierarchy[-1]\n",
    "\n",
    "# Define the function space and print the DOFs\n",
    "V = FunctionSpace(mesh, \"Lagrange\", degree)\n",
    "dofs = V.dim()\n",
    "parprint('DOFs', dofs)\n",
    "\n",
    "u = TrialFunction(V)\n",
    "v = TestFunction(V)\n",
    "\n",
    "bcs = DirichletBC(V, zero(), (\"on_boundary\",))\n",
    "\n",
    "# Define the RHS and analytic solution\n",
    "x, y, z = SpatialCoordinate(mesh)\n",
    "\n",
    "a = Constant(1)\n",
    "b = Constant(2)\n",
    "exact = sin(pi*x)*tan(pi*x/4)*sin(a*pi*y)*sin(b*pi*z)\n",
    "truth = Function(V).interpolate(exact)\n",
    "f = -pi**2 / 2\n",
    "f *= 2*cos(pi*x) - cos(pi*x/2) - 2*(a**2 + b**2)*sin(pi*x)*tan(pi*x/4)\n",
    "f *= sin(a*pi*y)*sin(b*pi*z)\n",
    "\n",
    "# Define the problem using the bilinear form `a` and linear functional `L`\n",
    "a = dot(grad(u), grad(v))*dx\n",
    "L = f*v*dx\n",
    "u_h = Function(V)\n",
    "problem = LinearVariationalProblem(a, L, u_h, bcs=bcs)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8f66f00f",
   "metadata": {},
   "source": [
    "Creating a problem instance we can see there are just short of 275000 DOFs in this noteook example."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b1052160",
   "metadata": {},
   "source": [
    "## The Solver\n",
    "\n",
    "This table summarises the different solvers we will use:\n",
    "\n",
    "| Solver | Abbreviation |   Cost   | Information |\n",
    "|:-------|:------------:|:--------:|:------------|\n",
    "| LU     | LU           | O(n³)*   | Firedrake Default |\n",
    "| Conjugate Gradient + Algebraic Multigrid | CG + AMG | O(qn) | Sensible choice of KSP + PC |\n",
    "| Conjugate Gradient + Geometric Multigrid V-cycle | CG + GMG V-cycle | O(qn) | GMG in place of AMG |\n",
    "| Full Geometric Multigrid | CG + Full GMG | O(qn) |  |\n",
    "| Matrix free FMG with Telescoping | Matfree CG + telescoped full GMG | O(qn) | Reduced memory and communication |\n",
    "\n",
    "*See discussion at the end of the LU section\n",
    "\n",
    "The n in the above table is the problem size (number of DOFs) and q is the number of iterations taken by an iterative method. In this notebook we use multigrid preconditioners to try and minimise the number of iterations, q.\n",
    "\n",
    "We define a function to wrap the solve, so we can provide different solver options and to assess their performance, the run time is printed.\n",
    "This is a fairly crude way to profile our code, for a more in depth guide to profiling, take a look at the page on [optimising Firedrake performance](https://firedrakeproject.org/optimising.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7cfcd7c2",
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_solve(problem, parameters):\n",
    "    # Create a solver and time how long the solve takes\n",
    "    t = time()\n",
    "    solver = LinearVariationalSolver(problem, solver_parameters=parameters)\n",
    "    solver.solve()\n",
    "    parprint(\"Runtime :\", time() - t)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a2870011",
   "metadata": {},
   "source": [
    "## LU\n",
    "\n",
    "We can start by looking at the Firedrake's default solver options. If you don't specify any solver options a direct solver such as MUMPS will be used to perform an LU factorisation.\n",
    "\n",
    "Here we explicitly list the PETSc solver options so it's clear how the solver is set up. We also enable the `snes_view` so that PETSc prints the solver options it's using at runtime.\n",
    "\n",
    "**Warning:** This cell will take a long time (>2 minutes) to execute interactively!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "171d0ebf",
   "metadata": {},
   "outputs": [],
   "source": [
    "u_h.assign(0)\n",
    "lu_mumps = {\n",
    "    \"snes_view\": None,\n",
    "    \"ksp_type\": \"preonly\",\n",
    "    \"pc_type\": \"lu\",\n",
    "    \"pc_factor_mat_solver_type\": \"mumps\"\n",
    "}\n",
    "run_solve(problem, lu_mumps)\n",
    "parprint(\"Error   :\", errornorm(truth, u_h))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "43f08260",
   "metadata": {},
   "source": [
    "The above solve takes under a minute on a single Zen2 core of ARCHER2.\n",
    "\n",
    "Dense LU factorisations are very expensive, typically $O(n^3)$. Sparse LU with a state of the art solver like MUMPS or SuperLU_dist can do better, typically in $O(n^2)$ for a 3D finite element matrix or $O(n^{3/2})$ in 2D. For specific problems it may be possible to reduce that complexity even further.\n",
    "\n",
    "We can measure the computational cost of our solvers by increasing the problem size (the number of DOFs) and observing how this changes the solver run time. In the computational cost plots below you can see that the cost of LU factorisation is approximately $O(n^{5/3})$, but this cost grows far faster than the other solver methods.\n",
    "\n",
    "Direct solvers are very fast for small problems, which is why LU is the default solver in Firedrake. However, when $n$ gets large, direct solvers are no longer viable and should be avoided where possible.\n",
    "\n",
    "![](image/hpc_single.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "09b97914",
   "metadata": {},
   "source": [
    "## Iterative solvers\n",
    "\n",
    "An alternative to a direct solver is an iterative solver and PETSc gives us access to a large number of Krylov Subspace solvers (KSP). Since we have a symmetric problem, we can use the Conjugate Gradient (CG) method, which has computational cost $O(qn)$, where $q$ is the number of iterations for the method to converge. \n",
    "\n",
    "To reduce $q$ we can precondition the KSP, here we use PETSc's `gamg` Algebraic Multigrid (AMG) as a preconditioner.\n",
    "\n",
    "We assign 0 to the function `u_h` before we solve so that we aren't using the solution from the LU solve above as our initial guess for the CG solver."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2cd27a39",
   "metadata": {},
   "outputs": [],
   "source": [
    "u_h.assign(0)\n",
    "cg_amg = {\n",
    "    \"snes_view\": None,\n",
    "    \"ksp_type\": \"cg\",\n",
    "    \"pc_type\": \"gamg\",\n",
    "    \"pc_mg_log\": None\n",
    "}\n",
    "run_solve(problem, cg_amg)\n",
    "parprint(\"Error   :\", errornorm(truth, u_h))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d67053cd",
   "metadata": {},
   "source": [
    "Looking at the code where we defined the problem, in **The equations** section above, we have created a `MeshHierarchy` which allows for the use of Geometric Multigrid V-cycles to precondition the CG method within Firedrake. The solver options for this setup are shown below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "698bb7c1",
   "metadata": {},
   "outputs": [],
   "source": [
    "u_h.assign(0)\n",
    "cg_gmg_v = {\n",
    "    \"snes_view\": None,\n",
    "    \"ksp_type\": \"cg\",\n",
    "    \"pc_type\": \"mg\",\n",
    "    \"pc_mg_log\": None\n",
    "}\n",
    "run_solve(problem, cg_gmg_v)\n",
    "parprint(\"Error   :\", errornorm(truth, u_h))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c4955120",
   "metadata": {},
   "source": [
    "The CG solver with AMG or GMG V-cycles is significantly faster than the LU factorisation, but is still slower than using the full Geometric multigrid method, which we discuss in the next section.\n",
    "\n",
    "We can measure the weak scaling performance of the solvers by increasing the size of the problem in line with the number of processors. This is done approximately in the plot below, the number of DOFs per core is displayed underneath each data point. For a solver that weak scales perfectly, when we use twice as many cores to solve a problem twice as large, the total runtime should be the same and the lines in the plot should be approximately constant.\n",
    "\n",
    "In the weak scaling plot below, CG + GMG V-cycles weak scales for longer than the LU factorisation and CG + AMG does even better, but we also see that in _this_ setup we can do even better with full GMG methods.\n",
    "\n",
    "![](image/hpc_weak.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3f5b7c68",
   "metadata": {},
   "source": [
    "## Geometric Multigrid\n",
    "\n",
    "Using the multigrid hierarchy is possible to solve the Poisson problem using full multigrid sweeps (sometimes called F-cycles).\n",
    "\n",
    "By carefully choosing the number of smoothing steps (`mg_levels_ksp_max_it`) the number of CG iterations can be minimised."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9750d74b",
   "metadata": {},
   "outputs": [],
   "source": [
    "u_h.assign(0)\n",
    "fmg = {\n",
    "    \"snes_view\": None,\n",
    "    \"ksp_type\": \"cg\",\n",
    "    \"pc_type\": \"mg\",\n",
    "    \"pc_mg_log\": None,\n",
    "    \"pc_mg_type\": \"full\",\n",
    "    \"mg_levels_ksp_type\": \"chebyshev\",\n",
    "    \"mg_levels_ksp_max_it\": 2,\n",
    "    \"mg_levels_pc_type\": \"jacobi\",\n",
    "    \"mg_coarse_pc_type\": \"lu\",\n",
    "    \"mg_coarse_pc_factor_mat_solver_type\": \"mumps\"\n",
    "}\n",
    "run_solve(problem, fmg)\n",
    "parprint(\"Error   :\", errornorm(truth, u_h))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0e10eb44",
   "metadata": {},
   "source": [
    "Using full GMG gives a significant speed up over using multigrid V-cycles as a preconditioner.\n",
    "\n",
    "We can measure the strong scaling performance of the multigrid by choosing a large enough problem and seeing how long it takes to solve on different numbers of processes. In the plot below, the number of DOFs per core is displayed underneath each data point. For a solver that strong scales perfectly, when we use twice as many cores to solve the same size problem, the total runtime should be halved. This perfect scaling is plotted as a dashed line for comparison.\n",
    "\n",
    "The figure below shows what happens when we use this multigrid solver for a large problem. For this test we set `Nx = 10` and `Nref = 4` to make a problem with 33 076 161 DOFs and solve over multiple nodes.\n",
    "\n",
    "The full multigrid solver strong scales poorly beyond 2 nodes, the reason for this poor scaling is that the solver spends most of its time performing communication solving the problem on the coarse grid in a distributed manner. CG + AMG scales much better, but isn't as fast as using telescoping, which we discuss in the next section.  Designing a solver that is both fast and scalable for a given problem is often very challenging.\n",
    "\n",
    "![](image/hpc_strong.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ea794366",
   "metadata": {},
   "source": [
    "## Matrix free and telescoping\n",
    "\n",
    "In this section we show a final variation of the full multigrid solver above, which has advantages for larger problems and on HPC architectures.\n",
    "\n",
    "One key advantage of using geometric multigrid over algebraic multigrid is the ability to use matrix free methods. These methods never assemble the full finite element matrix, which for large problems gives a significant reduction in memory usage. More information on matrix free methods in Firedrake can be found in the [documentation](https://www.firedrakeproject.org/matrix-free.html). On the coarsest mesh of the multigrid hierarchy we can use the `firedrake.AssembledPC` to assemble the finite element matrix, which allows us to use a direct solver.\n",
    "\n",
    "The final set of solver options also deals with very large problems spread over multiple compute nodes. For a problem with a large multigrid hierarchy, the coarse grid problem is often so small that when it is solved over multiple nodes, the coarse solve spends all its time performing communication, which is slow.\n",
    "\n",
    "The solution is to let each node solve a local copy of the coarse grid problem, which avoids this communication. This functionality is enabled using the `telescope` preconditioner alongside the assembled preconditioner, as shown below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d11e2fad",
   "metadata": {},
   "outputs": [],
   "source": [
    "u_h.assign(0)\n",
    "telescope_factor = 1 # Set to number of nodes!\n",
    "fmg_matfree_telescope = {\n",
    "    \"snes_view\": None,\n",
    "    \"mat_type\": \"matfree\",\n",
    "    \"ksp_type\": \"cg\",\n",
    "    \"pc_type\": \"mg\",\n",
    "    \"pc_mg_log\": None,\n",
    "    \"pc_mg_type\": \"full\",\n",
    "    \"mg_levels_ksp_type\": \"chebyshev\",\n",
    "    \"mg_levels_ksp_max_it\": 2,\n",
    "    \"mg_levels_pc_type\": \"jacobi\",\n",
    "    \"mg_coarse_pc_type\": \"python\",\n",
    "    \"mg_coarse_pc_python_type\": \"firedrake.AssembledPC\",\n",
    "    \"mg_coarse_assembled\": {\n",
    "        \"mat_type\": \"aij\",\n",
    "        \"pc_type\": \"telescope\",\n",
    "        \"pc_telescope_reduction_factor\": telescope_factor,\n",
    "        \"pc_telescope_subcomm_type\": \"contiguous\",\n",
    "        \"telescope_pc_type\": \"lu\",\n",
    "        \"telescope_pc_factor_mat_solver_type\": \"mumps\"\n",
    "    }\n",
    "}\n",
    "run_solve(problem, fmg_matfree_telescope)\n",
    "parprint(\"Error   :\", errornorm(truth, u_h))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "14c92943",
   "metadata": {},
   "source": [
    "## Running on HPC\n",
    "\n",
    "To run these examples on HPC, the Firedrake code must be a Python script. Copy and paste these notebook cells into a text editor on the remote machine and save it as a Python script (extension `.py`).\n",
    "\n",
    "The code must run through a job scheduler using another script. An example job script suitable for running on ARCHER2 is provided below.\n",
    "\n",
    "To use this script change the account (`-A`) to your account, change the number of nodes (`--node=`) to the number of nodes you want to use and the time (`-t`) as appropriate, it is currently set to 10 _minutes_."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "edaf573d",
   "metadata": {},
   "source": [
    "```bash\n",
    "#!/bin/bash\n",
    "#SBATCH -p standard\n",
    "#SBATCH -A account\n",
    "#SBATCH -J firedrake\n",
    "#SBATCH --nodes=1\n",
    "#SBATCH --cpus-per-task=1\n",
    "#SBATCH --qos=standard\n",
    "#SBATCH -t 0:10:00\n",
    "\n",
    "export VENV_NAME=firedrake_08_2021\n",
    "export WORK=/work/e682/shared/firedrake_tarballs/firedrake_08_2021/\n",
    "export FIREDRAKE_TEMP=firedrake_tmp\n",
    "export LOCAL_BIN=$WORK\n",
    "\n",
    "myScript=\"HPC_demo.py\"\n",
    "\n",
    "module load epcc-job-env\n",
    "\n",
    "# Activate Firedrake venv (activate once on first node, extract once per node)\n",
    "source $LOCAL_BIN/firedrake_activate.sh\n",
    "srun --ntasks-per-node 1 $LOCAL_BIN/firedrake_activate.sh\n",
    "\n",
    "# Run Firedrake script\n",
    "srun --ntasks-per-node 128 $VIRTUAL_ENV/bin/python ${myScript}\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "abf8e1fc",
   "metadata": {},
   "source": [
    "If you named your jobscript `jobscript.slm`, then it can be submitted to the queue by running the following command on ARCHER2:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7d28917b",
   "metadata": {},
   "source": [
    "```bash\n",
    "sbatch jobscript.slm\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3bd2ea03",
   "metadata": {},
   "source": [
    "You can see your job's progress through the queue using:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1aad5278",
   "metadata": {},
   "source": [
    "``` bash\n",
    "squeue -u $USER\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7de92d5e",
   "metadata": {},
   "source": [
    "If you need to cancel a job for any reason, you can pass your job ID number as an argument to the scancel command:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1bd14f6b",
   "metadata": {},
   "source": [
    "``` bash\n",
    "scancel 123456\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "de06f2b8",
   "metadata": {},
   "source": [
    "Once your job has completed any output will be stored in files named `slurm-123456.out` and `slurm-123456.err`. The job ID `123456` is used as an example here, yours will be different each time you run a job."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "67c2827c",
   "metadata": {},
   "source": [
    "## Exercise\n",
    "\n",
    "Perform a convergence study for the Poisson problem above, using degree 2 Lagrange elements. To do this, solve the problem on a range of different mesh sizes. The cell diameter on the finest mesh in a multigrid hierarchy is given by $h = \\frac{\\sqrt{2}}{N}$, where $N = N_x \\times 2^{N_{ref}}$ is the number of cells along one edge of the cube on the finest grid.\n",
    "\n",
    "**Note:** If you're following along as part of a tutorial you will be assigned a single grid size and this exercise will be completed as a group."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9f284430",
   "metadata": {},
   "source": [
    "a)\n",
    "\n",
    "Choose an appropriate number of multigrid levels (`Nref`) and coarse grid size (`Nx`) for each mesh size `N`. For this exercise we will repeatedly double $N$ (to half the value of $h$), and measure the error for each solution. Use your answers to populate the table below:\n",
    "\n",
    "| N =  | 8 | 16 | 32 | 64 | 128 | 256 | 512 |\n",
    "|------|---|----|----|----|-----|-----|-----|\n",
    "| Nx   |   |    | 8  |    |     |     |     |\n",
    "| Nref |   |    | 2  |    |     |     |     |\n",
    "\n",
    "Throughout the exercise we have already entered appropriate values into the table. These values correspond to the case presented in the notebook."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "65fb7931",
   "metadata": {},
   "source": [
    "b)\n",
    " \n",
    "Calculate the number of DOFs for each problem size using the formula in the **How big?** section above. Use the total number of DOFs to work out how many processes would be appropriate for solving each problem size (try to pick a power of 2) and hence how many nodes you require for that simulation. Place all these values in the table:\n",
    "\n",
    "| N =       | 8 | 16 | 32     | 64 | 128 | 256 | 512 |\n",
    "|-----------|---|----|--------|----|-----|-----|-----|\n",
    "| DOFs      |   |    | 274625 |    |     |     |     |\n",
    "| Processes |   |    | 4      |    |     |     |     |\n",
    "| Nodes     |   |    | 1      |    |     |     |     |"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d88877d5",
   "metadata": {},
   "source": [
    "c)\n",
    "\n",
    "For each problem size (or your given problem size if you are in a group) we will execute a Python script on the HPC to solve the Poisson problem.\n",
    "\n",
    "Copy the cell containing the submission script in the **Running on HPC** above to your text editor on the HPC. Using your answer to (b), edit the lines `#SBATCH --nodes=1` to the number of nodes you require for your problem size and the parameter `--ntasks-per-node` in the line:\n",
    "```\n",
    "srun --ntasks-per-node 128 $VIRTUAL_ENV/bin/python ${myScript}\n",
    "```\n",
    "to the number of processes you require. Save the file as `jobscript.slm`.\n",
    "\n",
    "Next we must create a Firedrake script to run on HPC. If you are following as part of a tutorial a template will be provided, otherwise you can copy and paste code from cells in the notebook. Edit the values of `Nx` and `Nref` in the script to solve your selected problem size using your answer to (a). Ensure you save the files as `HPC_demo.py`\n",
    "\n",
    "Finally, submit the job to the queue using the command `sbatch jobscript.slm` on the HPC command line and, once the job has run, check the output files current directory and fill in the error in the table below:\n",
    "\n",
    "\n",
    "| N =   | 8 | 16 | 32       | 64 | 128 | 256 | 512 |\n",
    "|-------|---|----|----------|----|-----|-----|-----|\n",
    "| h     |   |    | 0.044    |    |     |     |     |\n",
    "| Error |   |    | 1.74E-05 |    |     |     |     |\n",
    "\n",
    "If you are performing the convergence study individually, continue editing the scripts to populate the rest of the table. Both the Python script and jobscript need to be changed to suit the problem size!"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e43035fd",
   "metadata": {},
   "source": [
    "d)\n",
    "\n",
    "Plot the error against h and measure the rate of convergence using `matplotlib`. If you are completing this as part of a tutorial submit your results from (c) to the instructor and they will combine the results and plot the graph."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "df5d73a6",
   "metadata": {},
   "source": [
    "**Hints:**\n",
    "- You don't need much compute power to solve small problems on coarse meshes, these will likely fit on one node.\n",
    "- Remember to make your job big enough for the number of processes that you run:\n",
    "    - Each MPI rank must own at least one cell in the mesh.\n",
    "    - Firedrake performs better when there are more than 50000 DOFs per rank."
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}