<commit_before>{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Creative Commons CC-BY 4.0, code under MIT license (c) 2016 <NAME>, <NAME>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using AmgX to Accelerate PETSc-Based CFD Codes\n",
    "=================================================="
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "**<NAME> & <NAME>**, George Washington University.\n",
    "\n",
    "This technical report supplements the presentation by <NAME> at the [2016 GPU Technology Conference](http://www.gputechconf.com), on Thu. April 7th: see session [S6355](http://mygtc.gputechconf.com/quicklink/W7yvg9)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "-------------------"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[PetIBM](https://github.com/barbagroup/PetIBM) is a parallel code to study flow around immersed bodies, solving the Navier-Stokes equations on a Cartesian mesh. The body can be embedded in the fluid mesh via the so-called _immersed boundary method_ (IBM), which brings many advantages in regards to mesh generation. The equations are discretized with a finite-difference scheme, via a pressure-projection method, resulting in a Poisson-like linear system. PetIBM uses the [PETSc](https://www.mcs.anl.gov/petsc/) library to solve the linear system of equations in parallel on CPU hardware, and now we want to accelerate it with the help of GPU computing. \n",
    "\n",
    "Our group previously developed another code, cuIBM, which implements the same methods as PetIBM, except that it runs on a sigle CPU and a single GPU, using the [CUSP](http://cusplibrary.github.io) library to solve the linear system. The limited memory on a single GPU restricts the size of CFD problems that we can solve using cuIBM to about 4 million mesh points. In three dimensions, we often need to solve larger problems, which led us to develop PetIBM. Our next step was to add the capability of exploiting multi-GPU systems. \n",
    "\n",
    "In this work, we set out to accelerate PetIBM using the NVIDIA [AmgX](https://developer.nvidia.com/amgx) library on multi-GPUs. We wrote a wrapper code to couple AmgX and PETSc, so that we require few modifications in PetIBM. Our AmgX wrapper is able to handle the situation where the PETSc-based code launches more MPI processes than the number of GPUs available, achieving heterogeneous speed-up. It is open source and we hope others may find it useful to accelerate their PETSc-based parallel programs with multiple GPUs."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Brief overview of PetIBM\n",
    "-----------------------------------------"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We all know that solving Poisson problems is computationally demanding. It is particularly challenging in unsteady computational fluid dynamics (CFD), where we need to solve the Poisson system over and over again during a time-marching iteration. The situation in PetIBM is tougher than typical unsteady CFD problems because the coefficient matrix is not a standard discrete Laplace operator: it is augmented by terms derived from the presence of the body. We are in fact solving *modified Poisson equation* at every time step, and this is responsible for about 90% of the run time."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The particular immersed-boundary formulation in PetIBM is due to Taira and Colonius [1]. The governing equations for the fluid velocity, $\\bf{u}$, and pressure, $p$, are the following:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\left\\{\n",
    "\\begin{aligned}\n",
    "& \\frac{\\partial \\bf{u}}{\\partial t} + \\bf{u}\\cdot\\nabla\\bf{u} = \n",
    "-\\nabla\\it{p}+\\frac{1}{Re}\\nabla^2\\bf{u}+\\int_s{\\bf{f}(\\bf{\\xi}(\\it{s}, \\it{t}))\\delta(\\bf{\\xi}-\\bf{x})}d\\it{s}, \\\\\n",
    "& \\nabla\\cdot\\bf{u}=0\\it{,} \\\\\n",
    "& \\bf{u}(\\bf{\\xi}(\\it{s}, \\it{t})) = \n",
    "\\int_\\bf{x}{\\bf{u}(\\bf{x})\\delta(\\bf{x}-\\bf{\\xi})}d\\bf{x}=\\bf{u}_\\it{B}(\\bf{\\xi}(\\it{s}, \\it{t}))\n",
    "\\end{aligned}\n",
    "\\right.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the first equation, Navier-Stokes is augmented by a source term that represents the effect of the body on the fluid, via the no-slip boundary condition. Here, $\\bf{f}$ is a singular force distribution along the solid boundary. The second equation enforces the conservation of mass in an incompressible fluid, and the third equation states that the fluid velocity is the same as the body velocity, $\\bf{u}_\\it{B}$ at the boundary (no slip). After discretization, we can write the above equations as a block system of linear equations, in which each submatrix represents an operator or a term in the PDE:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\left[\n",
    "\\begin{matrix}\n",
    "A & G & E^T \\\\\n",
    "G^T & 0 & 0 \\\\\n",
    "E & 0 & 0\n",
    "\\end{matrix}\n",
    "\\right]\n",
    "\\left(\n",
    "\\begin{matrix}\n",
    "q^{n+1} \\\\\n",
    "\\phi \\\\\n",
    "\\tilde{f}\n",
    "\\end{matrix}\n",
    "\\right)\n",
    "=\n",
    "\\left(\n",
    "\\begin{matrix}\n",
    "r^n \\\\\n",
    "0 \\\\\n",
    "u_B^{n+1}\n",
    "\\end{matrix}\n",
    "\\right)\n",
    "+\n",
    "\\left(\n",
    "\\begin{matrix}\n",
    "bc_1 \\\\\n",
    "-bc_2 \\\\\n",
    "0\n",
    "\\end{matrix}\n",
    "\\right)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, $q^{n+1}$ are the momentum fluxes at cell boundaries on the next time step, $n+1$, while $r^n$ are the convective terms on the previous time step; $\\phi$ is a vector of pressure values and $\\tilde{f}$ is a vector of force values at the boundary points. We can further define some larger submatrices and perform a block-LU decomposition:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "Q\\equiv\\left[G, E^T\\right],\\ \n",
    "\\lambda\\equiv\\left(\\begin{matrix}\\phi \\\\ \\tilde{f}\\end{matrix}\\right),\\ \n",
    "r_1\\equiv r^n + bc_1,\\ \n",
    "r_2\\equiv\\left(\\begin{matrix}-bc_2 \\\\ u_B^{n+1}\\end{matrix}\\right)\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\Rightarrow\n",
    "\\left[\n",
    "\\begin{matrix}\n",
    "A & 0 \\\\\n",
    "Q^T & -Q^TB^NQ\n",
    "\\end{matrix}\n",
    "\\right]\n",
    "\\left[\n",
    "\\begin{matrix}\n",
    "I & B^NQ \\\\\n",
    "0 & I\n",
    "\\end{matrix}\n",
    "\\right]\n",
    "\\left(\n",
    "\\begin{matrix}\n",
    "q^{n+1} \\\\\n",
    "\\lambda\n",
    "\\end{matrix}\n",
    "\\right)\n",
    "=\n",
    "\\left(\n",
    "\\begin{matrix}\n",
    "r_1 \\\\\n",
    "r_2\n",
    "\\end{matrix}\n",
    "\\right)\n",
    "+\n",
    "\\left(\n",
    "\\begin{matrix}\n",
    "-\\frac{\\Delta t^N}{2^N}\\left(LM^{-1}\\right)Q\\lambda \\\\\n",
    "0\n",
    "\\end{matrix}\n",
    "\\right)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above matrix system corresponds to solving three algebraic equations:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{aligned}\n",
    "& Aq^*=r_1 \\\\\n",
    "& Q^TB^NQ\\lambda=Q^Tq^*-r_2 \\\\\n",
    "& q^{n+1}=q^*-B^NQ\\lambda\n",
    "\\end{aligned}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, $q^{*}$ is the intermediate solution for the fluxes, used in the pressure projection step. The second equation is a *modified Poisson equation*, corresponding to a Poisson equation plus terms that are due to the presence of the body. The coefficient matrix in the modified Poisson eqauation, $Q^TB^NQ$, has a non-zero pattern like that shown in Figure 1, which corresponds to the situation where a circular cylinder is immersed in a fluid."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"plots/cylinder3dRe40QTBNQ.png\" style=\"width:300px;height:300px;\" />\n",
    "\n",
    "##### Figure 1. Non-zero pattern of the modified Poisson coefficient matrix."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The upper-left block in the coefficient matrix corresponds to the typical Poisson system, i.e., the discrete Laplace operator. The bottom, bottom-right corner and the right side of the matrix in Figure 1 are the modification terms due to the presence of the body."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In PetIBM, solving this modified Poisson equation takes over **90%** of the total run time. We will introduce GPU computing to accelerate only this part of the overall calculation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## AmgX\n",
    "------------"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[AmgX](https://developer.nvidia.com/amgx) is a library that provides linear solvers and preconditioners with capability of using multiple GPUs. It is developed and supported by NVIDIA, having entered beta production in January 2014. A tutorial, user's manual and license information are available on the [NVIDIA Developer website](https://developer.nvidia.com/amgx).\n",
    "\n",
    "Here we list some of the features of AmgX:\n",
    "\n",
    "* Krylov methods:\n",
    "    * CG, GMRES, BiCGStab, ... etc\n",
    "* Multigrid (MG) preconditioners:\n",
    "    * Classical algebraic multigrid (AMG)\n",
    "    * Unsmoothed aggregation AMG\n",
    "* Smoothers:\n",
    "    * Block-Jacobi, Gauss-Seidel, incomplete LU, dense LU, ... etc\n",
    "* Cycles:\n",
    "    * V, W, F, FG, CGF\n",
    "* Multiple GPUs on single node / multiple nodes\n",
    "    * MPI (OpenMPI) / MPI Direct\n",
    "    * Single MPI rank <=> single GPU\n",
    "    * Multiple MPI ranks <=> single GPU\n",
    "* C API\n",
    "* Unified Virtual Addressing\n",
    "\n",
    "The various combinations of Krylov methods and multigrid preconditioners are very popular in CFD applications. In fact, the AmgX library was developed specifically aimed for CFD and similar computational physics applications, in a partnership between NVIDIA and ANSYS, a major provider of commercial CFD software.\n",
    "\n",
    "[Benchmarks](https://devblogs.nvidia.com/parallelforall/amgx-v1-0-enabling-reservoir-simulation-with-classical-amg/) published by NVIDIA in the Developer Blog (May 15, 2014) compared against the \"gold standard\" classical AMG implementation: BoomerAMG in the HYPRE package. They showed speedups of 2x in the setup phase and 4–6x in the solve phase, when comparing a single K40 GPU versus 2x10-core Intel Ivy Bridge CPUs, on problems of size 5 million and 10 million. But results on multi-GPUs were still unclear."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Our AmgX wrapper\n",
    "-----------------------"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Why did we need to write an AmgX \"wrapper\" code?\n",
    "\n",
    "AmgX requires the application code to provide it with the coefficient matrix, and the right-hand-side vector of the linear system. The application can create these data structures in the host CPU (then copy them to the device), or in-place in the GPU.\n",
    "\n",
    "But our application code uses the PETSc library. Both AmgX and PETSc have their own data structures for parallel linear algebra. PetIBM is also written in C++ with some object-oriented design, while AmgX offers a pure-C API. Thus, to use AmgX directly in PetIBM, we would need to add or modify a lot of code, which we would rather avoid. \n",
    "\n",
    "Also, we may later on want to develop some other software using PETSc, and we may also want to use AmgX with it. So a wrapper to couple AmgX and PETSc could help us or others in the future."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The philosophy of our wrapper is simple: except when initializing and finalizing, users should only need to call `setA` and `solve` in most cases. In unsteady CFD applications, with time-marching iterations, users would simply call `solve` repeatedly.  But there should be no need to be concerned about data conversion: `A`, `x`, and `rhs` are all PETSc data structures, and remain so, on the application side. We also would like configuration to be done simply through an input script during runtime.\n",
    "\n",
    "The wrapper code equips the PETSc application software with access to AmgX, without the burden of unpacking data structures from one into the other. Without it, to use AmgX from PETSc one would have to access the raw data from the sparse matrix in PETSc, then redistribute as needed across MPI processes to consolidate the data in a way that can be sent to the GPU device. Next, the application would need to initialize AmgX and allocate memory on the GPU—using the  API only, users would need about 10 different AmgX functions for just this part. Users would need more code to set up the matrix and vectors and send from host to device, then set up the solver, solve and finalize.  With our wrapper, users only need to declare an instance of the solver with `AmgXWrapper solver` and then call `solver.initialize(...)`. The diagram in Figure 2 shows the simple work flow provided by our wrapper."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"plots/wrapper.png\" style=\"width:500px;height:300px;\" />\n",
    "\n",
    "##### Figure 2. Diagram of our AmgX wrapper work flow."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numerical tests\n",
    "-----------------------------------------------"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### First Test: 2D Cylinder Flow, Re=40"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So far, the wrapper is just a wrapper. It's just an interface between PETSc ans AmgX and should work well. Let's see how it works. First we run tests with a 2D cylinder flow (Re=40). Mesh size is 2.25M.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"plots/rawAmgXMultiRanks_Cylinder2d_2.png\" style=\"width:500px;height:500px;\" />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The far left bar is the result from original CPU-version PetIBM with 6 CPU cores. We can see that solving modified Poisson equation, which is the blue portion, takes over 90% of total time in this case. \n",
    "\n",
    "The middle bar is the result from using the AmgX-version PetIBM with 1 GPU and 1 CPU core. We get 4 times speed-up on solving modified Poisson equation and 3 times speed-up for total time. We think it looks great, but not great enough. From the parts that are solved on CPUs, which are the portions with colors other than blue, we find now it takes more time to solve them. The reason is obvious: in this case we only use 1 CPU core to solve these CPU-based parts, while 6 CPU cores were used to solve them in the far left bar. Let's see how if we use 6 CPU cores to solve the CPU-based parts and 1 GPU to solve the modified Poisson system.\n",
    "\n",
    "The far right bar is the result. We indeed get the same performances on the parts solved or calculated with CPUs, but the speed-up of solving modified Poisson equation is not good as we expected. Actually, it's totally unacceptable for me. So, what happened?\n",
    "\n",
    "It's not difficult to answer it. Just think about what will happen if you solve a linear system with domain-decomposition methods and partition the system into 6 subdomains but solve them with only one CPU core or less than 6 CPU cores. Situation here is the same.\n",
    "\n",
    "Recall that AmgX solvers only solve the linear system. That means we need to prepare the linear systems on CPUs and using PETSc. So when we launch PetIBM with 6 CPU cores, the modified Poisson system is partitioned into 6 subdomains through some partitioning methods in PETSc. And when we call AmgX to solve the modified Poisson system, each CPU core creates an AmgX solver to solve the subdomain it holds. In our case, there are 6 subdomain-solvers but only one GPU. This means the 6 subdomain-solvers are competing GPU resources. And remember we are using domain-decomposition methods, so information exchanges occurs between these 6 subdomain-solvers and on one single GPU. Solving modified Poisson equation in the far right bar therefore is not as fast as we expected."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### More examples revealing this issue:\n",
    "\n",
    "Let's see more examples revealing this issue:\n",
    "\n",
    "#### 1. Solving 3D Poisson equations (6M unknowns; CG solver; Classical AMG preconditioner):\n",
    "\n",
    "<img src=\"plots/rawAmgXMultiRanks_Classical.png\" style=\"width:450px;height:300px;\" />\n",<commit_msg>add details on first test<commit_after>{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Creative Commons CC-BY 4.0, code under MIT license (c) 2016 <NAME>, <NAME>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using AmgX to Accelerate PETSc-Based CFD Codes\n",
    "=================================================="
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "**<NAME> & <NAME>**, George Washington University.\n",
    "\n",
    "This technical report supplements the presentation by <NAME> at the [2016 GPU Technology Conference](http://www.gputechconf.com), on Thu. April 7th: see session [S6355](http://mygtc.gputechconf.com/quicklink/W7yvg9)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "-------------------"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[PetIBM](https://github.com/barbagroup/PetIBM) is a parallel code to study flow around immersed bodies, solving the Navier-Stokes equations on a Cartesian mesh. The body can be embedded in the fluid mesh via the so-called _immersed boundary method_ (IBM), which brings many advantages in regards to mesh generation. The equations are discretized with a finite-difference scheme, via a pressure-projection method, resulting in a Poisson-like linear system. PetIBM uses the [PETSc](https://www.mcs.anl.gov/petsc/) library to solve the linear system of equations in parallel on CPU hardware, and now we want to accelerate it with the help of GPU computing. \n",
    "\n",
    "Our group previously developed another code, cuIBM, which implements the same methods as PetIBM, except that it runs on a sigle CPU and a single GPU, using the [CUSP](http://cusplibrary.github.io) library to solve the linear system. The limited memory on a single GPU restricts the size of CFD problems that we can solve using cuIBM to about 4 million mesh points. In three dimensions, we often need to solve larger problems, which led us to develop PetIBM. Our next step was to add the capability of exploiting multi-GPU systems. \n",
    "\n",
    "In this work, we set out to accelerate PetIBM using the NVIDIA [AmgX](https://developer.nvidia.com/amgx) library on multi-GPUs. We wrote a wrapper code to couple AmgX and PETSc, so that we require few modifications in PetIBM. Our AmgX wrapper is able to handle the situation where the PETSc-based code launches more MPI processes than the number of GPUs available, achieving heterogeneous speed-up. It is open source and we hope others may find it useful to accelerate their PETSc-based parallel programs with multiple GPUs."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Brief overview of PetIBM\n",
    "-----------------------------------------"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We all know that solving Poisson problems is computationally demanding. It is particularly challenging in unsteady computational fluid dynamics (CFD), where we need to solve the Poisson system over and over again during a time-marching iteration. The situation in PetIBM is tougher than typical unsteady CFD problems because the coefficient matrix is not a standard discrete Laplace operator: it is augmented by terms derived from the presence of the body. We are in fact solving *modified Poisson equation* at every time step, and this is responsible for about 90% of the run time."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The particular immersed-boundary formulation in PetIBM is due to Taira and Colonius [1]. The governing equations for the fluid velocity, $\\bf{u}$, and pressure, $p$, are the following:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\left\\{\n",
    "\\begin{aligned}\n",
    "& \\frac{\\partial \\bf{u}}{\\partial t} + \\bf{u}\\cdot\\nabla\\bf{u} = \n",
    "-\\nabla\\it{p}+\\frac{1}{Re}\\nabla^2\\bf{u}+\\int_s{\\bf{f}(\\bf{\\xi}(\\it{s}, \\it{t}))\\delta(\\bf{\\xi}-\\bf{x})}d\\it{s}, \\\\\n",
    "& \\nabla\\cdot\\bf{u}=0\\it{,} \\\\\n",
    "& \\bf{u}(\\bf{\\xi}(\\it{s}, \\it{t})) = \n",
    "\\int_\\bf{x}{\\bf{u}(\\bf{x})\\delta(\\bf{x}-\\bf{\\xi})}d\\bf{x}=\\bf{u}_\\it{B}(\\bf{\\xi}(\\it{s}, \\it{t}))\n",
    "\\end{aligned}\n",
    "\\right.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the first equation, Navier-Stokes is augmented by a source term that represents the effect of the body on the fluid, via the no-slip boundary condition. Here, $\\bf{f}$ is a singular force distribution along the solid boundary. The second equation enforces the conservation of mass in an incompressible fluid, and the third equation states that the fluid velocity is the same as the body velocity, $\\bf{u}_\\it{B}$ at the boundary (no slip). After discretization, we can write the above equations as a block system of linear equations, in which each submatrix represents an operator or a term in the PDE:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\left[\n",
    "\\begin{matrix}\n",
    "A & G & E^T \\\\\n",
    "G^T & 0 & 0 \\\\\n",
    "E & 0 & 0\n",
    "\\end{matrix}\n",
    "\\right]\n",
    "\\left(\n",
    "\\begin{matrix}\n",
    "q^{n+1} \\\\\n",
    "\\phi \\\\\n",
    "\\tilde{f}\n",
    "\\end{matrix}\n",
    "\\right)\n",
    "=\n",
    "\\left(\n",
    "\\begin{matrix}\n",
    "r^n \\\\\n",
    "0 \\\\\n",
    "u_B^{n+1}\n",
    "\\end{matrix}\n",
    "\\right)\n",
    "+\n",
    "\\left(\n",
    "\\begin{matrix}\n",
    "bc_1 \\\\\n",
    "-bc_2 \\\\\n",
    "0\n",
    "\\end{matrix}\n",
    "\\right)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, $q^{n+1}$ are the momentum fluxes at cell boundaries on the next time step, $n+1$, while $r^n$ are the convective terms on the previous time step; $\\phi$ is a vector of pressure values and $\\tilde{f}$ is a vector of force values at the boundary points. We can further define some larger submatrices and perform a block-LU decomposition:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "Q\\equiv\\left[G, E^T\\right],\\ \n",
    "\\lambda\\equiv\\left(\\begin{matrix}\\phi \\\\ \\tilde{f}\\end{matrix}\\right),\\ \n",
    "r_1\\equiv r^n + bc_1,\\ \n",
    "r_2\\equiv\\left(\\begin{matrix}-bc_2 \\\\ u_B^{n+1}\\end{matrix}\\right)\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\Rightarrow\n",
    "\\left[\n",
    "\\begin{matrix}\n",
    "A & 0 \\\\\n",
    "Q^T & -Q^TB^NQ\n",
    "\\end{matrix}\n",
    "\\right]\n",
    "\\left[\n",
    "\\begin{matrix}\n",
    "I & B^NQ \\\\\n",
    "0 & I\n",
    "\\end{matrix}\n",
    "\\right]\n",
    "\\left(\n",
    "\\begin{matrix}\n",
    "q^{n+1} \\\\\n",
    "\\lambda\n",
    "\\end{matrix}\n",
    "\\right)\n",
    "=\n",
    "\\left(\n",
    "\\begin{matrix}\n",
    "r_1 \\\\\n",
    "r_2\n",
    "\\end{matrix}\n",
    "\\right)\n",
    "+\n",
    "\\left(\n",
    "\\begin{matrix}\n",
    "-\\frac{\\Delta t^N}{2^N}\\left(LM^{-1}\\right)Q\\lambda \\\\\n",
    "0\n",
    "\\end{matrix}\n",
    "\\right)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above matrix system corresponds to solving three algebraic equations:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{aligned}\n",
    "& Aq^*=r_1 \\\\\n",
    "& Q^TB^NQ\\lambda=Q^Tq^*-r_2 \\\\\n",
    "& q^{n+1}=q^*-B^NQ\\lambda\n",
    "\\end{aligned}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, $q^{*}$ is the intermediate solution for the fluxes, used in the pressure projection step. The second equation is a *modified Poisson equation*, corresponding to a Poisson equation plus terms that are due to the presence of the body. The coefficient matrix in the modified Poisson eqauation, $Q^TB^NQ$, has a non-zero pattern like that shown in Figure 1, which corresponds to the situation where a circular cylinder is immersed in a fluid."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"plots/cylinder3dRe40QTBNQ.png\" style=\"width:300px;height:300px;\" />\n",
    "\n",
    "##### Figure 1. Non-zero pattern of the modified Poisson coefficient matrix."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The upper-left block in the coefficient matrix corresponds to the typical Poisson system, i.e., the discrete Laplace operator. The bottom, bottom-right corner and the right side of the matrix in Figure 1 are the modification terms due to the presence of the body."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In PetIBM, solving this modified Poisson equation takes over **90%** of the total run time. We will introduce GPU computing to accelerate only this part of the overall calculation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## AmgX\n",
    "------------"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[AmgX](https://developer.nvidia.com/amgx) is a library that provides linear solvers and preconditioners with capability of using multiple GPUs. It is developed and supported by NVIDIA, having entered beta production in January 2014. A tutorial, user's manual and license information are available on the [NVIDIA Developer website](https://developer.nvidia.com/amgx).\n",
    "\n",
    "Here we list some of the features of AmgX:\n",
    "\n",
    "* Krylov methods:\n",
    "    * CG, GMRES, BiCGStab, ... etc\n",
    "* Multigrid (MG) preconditioners:\n",
    "    * Classical algebraic multigrid (AMG)\n",
    "    * Unsmoothed aggregation AMG\n",
    "* Smoothers:\n",
    "    * Block-Jacobi, Gauss-Seidel, incomplete LU, dense LU, ... etc\n",
    "* Cycles:\n",
    "    * V, W, F, FG, CGF\n",
    "* Multiple GPUs on single node / multiple nodes\n",
    "    * MPI (OpenMPI) / MPI Direct\n",
    "    * Single MPI rank <=> single GPU\n",
    "    * Multiple MPI ranks <=> single GPU\n",
    "* C API\n",
    "* Unified Virtual Addressing\n",
    "\n",
    "The various combinations of Krylov methods and multigrid preconditioners are very popular in CFD applications. In fact, the AmgX library was developed specifically aimed for CFD and similar computational physics applications, in a partnership between NVIDIA and ANSYS, a major provider of commercial CFD software.\n",
    "\n",
    "[Benchmarks](https://devblogs.nvidia.com/parallelforall/amgx-v1-0-enabling-reservoir-simulation-with-classical-amg/) published by NVIDIA in the Developer Blog (May 15, 2014) compared against the \"gold standard\" classical AMG implementation: BoomerAMG in the HYPRE package. They showed speedups of 2x in the setup phase and 4–6x in the solve phase, when comparing a single K40 GPU versus 2x10-core Intel Ivy Bridge CPUs, on problems of size 5 million and 10 million. But results on multi-GPUs were still unclear."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Our AmgX wrapper\n",
    "-----------------------"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Why did we need to write an AmgX \"wrapper\" code?\n",
    "\n",
    "AmgX requires the application code to provide it with the coefficient matrix, and the right-hand-side vector of the linear system. The application can create these data structures in the host CPU (then copy them to the device), or in-place in the GPU.\n",
    "\n",
    "But our application code uses the PETSc library. Both AmgX and PETSc have their own data structures for parallel linear algebra. PetIBM is also written in C++ with some object-oriented design, while AmgX offers a pure-C API. Thus, to use AmgX directly in PetIBM, we would need to add or modify a lot of code, which we would rather avoid. \n",
    "\n",
    "Also, we may later on want to develop some other software using PETSc, and we may also want to use AmgX with it. So a wrapper to couple AmgX and PETSc could help us or others in the future."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The philosophy of our wrapper is simple: except when initializing and finalizing, users should only need to call `setA` and `solve` in most cases. In unsteady CFD applications, with time-marching iterations, users would simply call `solve` repeatedly.  But there should be no need to be concerned about data conversion: `A`, `x`, and `rhs` are all PETSc data structures, and remain so, on the application side. We also would like configuration to be done simply through an input script during runtime.\n",
    "\n",
    "The wrapper code equips the PETSc application software with access to AmgX, without the burden of unpacking data structures from one into the other. Without it, to use AmgX from PETSc one would have to access the raw data from the sparse matrix in PETSc, then redistribute as needed across MPI processes to consolidate the data in a way that can be sent to the GPU device. Next, the application would need to initialize AmgX and allocate memory on the GPU—using the  API only, users would need about 10 different AmgX functions for just this part. Users would need more code to set up the matrix and vectors and send from host to device, then set up the solver, solve and finalize.  With our wrapper, users only need to declare an instance of the solver with `AmgXWrapper solver` and then call `solver.initialize(...)`. The diagram in Figure 2 shows the simple work flow provided by our wrapper."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"plots/wrapper.png\" style=\"width:500px;height:300px;\" />\n",
    "\n",
    "##### Figure 2. Diagram of our AmgX wrapper work flow."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numerical tests\n",
    "-----------------------------------------------"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### First Test: 2D Cylinder Flow, Re=40"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With our AmgX wrapper code, we have access to the linear solvers on GPU provided by AmgX, from our PETSc application, PetIBM. As a first test, we used PetIBM to solve flow around a 2D cylinder at Reynolds number Re=40, on 2.25M mesh points, for 300 time steps. The detailed simulation settings can be found in the directory `data/cylinder2d_Re40`—the interested reader can repeat all of our tests (after installing PetIBM, AmgX and dependencies) using the corresponding data files for each run.\n",
    "\n",
    "The hardware used here was a workstation with a 6-core CPU (Intel i7-5930K) and up to 2 GPUs (NVIDIA K40c). First, we ran our original PetIBM code on the 6 CPU cores, then repeated the simulation using 1 CPU core and 1 GPU, and finally ran on all 6 CPU cores plus the GPU. The results in Figure 3 show the run time of each stage in the computation, for each system configuration.\n",
    "\n",
    "On the pure-CPU runs, both linear systems (velocity and pressure) are being solved with PETSc. In that case, we use a CG solver with block-Jacobi preconditioner for the velocity system, and a CG solver with GAMG preconditioner for the pressure system. With GPUs, only the pressure system is solved with AmgX; we use a CG solver with aggregation preconditioner in that case.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"plots/rawAmgXMultiRanks_Cylinder2d_2.png\" style=\"width:500px;height:500px;\" />\n",
    "\n",
    "##### Figure 3. Run time of PetIBM for flow around a 2D cylinder."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The far left bar is the result from original CPU-version PetIBM with 6 CPU cores. We can see that solving modified Poisson equation, which is the blue portion, takes over 90% of total time in this case. \n",
    "\n",
    "The middle bar is the result from using the AmgX-version PetIBM with 1 GPU and 1 CPU core. We get 4 times speed-up on solving modified Poisson equation and 3 times speed-up for total time. We think it looks great, but not great enough. From the parts that are solved on CPUs, which are the portions with colors other than blue, we find now it takes more time to solve them. The reason is obvious: in this case we only use 1 CPU core to solve these CPU-based parts, while 6 CPU cores were used to solve them in the far left bar. Let's see how if we use 6 CPU cores to solve the CPU-based parts and 1 GPU to solve the modified Poisson system.\n",
    "\n",
    "The far right bar is the result. We indeed get the same performances on the parts solved or calculated with CPUs, but the speed-up of solving modified Poisson equation is not good as we expected. Actually, it's totally unacceptable for me. So, what happened?\n",
    "\n",
    "It's not difficult to answer it. Just think about what will happen if you solve a linear system with domain-decomposition methods and partition the system into 6 subdomains but solve them with only one CPU core or less than 6 CPU cores. Situation here is the same.\n",
    "\n",
    "Recall that AmgX solvers only solve the linear system. That means we need to prepare the linear systems on CPUs and using PETSc. So when we launch PetIBM with 6 CPU cores, the modified Poisson system is partitioned into 6 subdomains through some partitioning methods in PETSc. And when we call AmgX to solve the modified Poisson system, each CPU core creates an AmgX solver to solve the subdomain it holds. In our case, there are 6 subdomain-solvers but only one GPU. This means the 6 subdomain-solvers are competing GPU resources. And remember we are using domain-decomposition methods, so information exchanges occurs between these 6 subdomain-solvers and on one single GPU. Solving modified Poisson equation in the far right bar therefore is not as fast as we expected."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### More examples revealing this issue:\n",
    "\n",
    "Let's see more examples revealing this issue:\n",
    "\n",
    "#### 1. Solving 3D Poisson equations (6M unknowns; CG solver; Classical AMG preconditioner):\n",
    "\n",
    "<img src=\"plots/rawAmgXMultiRanks_Classical.png\" style=\"width:450px;height:300px;\" />\n",
{'fraction_non_alphanumeric': 0.15286167425182734, 'fraction_numerical': 0.008495379947593435, 'mean_word_length': 2.9541934780237757, 'pattern_counts': {'":': 190, '<': 25, '<?xml version=': 0, '>': 25, 'https://': 12, 'lorem ipsum': 0, 'www.': 4, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '7527935', 'n_tokens_mistral': 11693, 'n_tokens_neox': 11049, 'n_words': 4444}