 |
deal.II version 9.8.0
|
Loading...
Searching...
No Matches
|
deal.II version 9.8.0
|
This tutorial depends on step-22, step-64.
| Table of contents | |
|---|---|
This program was contributed by Quang Hoang and Timo Heister, Clemson University.
We acknowledge support from the NVIDIA Academic Grant Program award "Sustainable NVIDIA GPU acceleration for open source FEM simulations in deal.II". Many thanks to Daniel Arndt, Rene Gassmöller, Martin Kronbichler, Ivan Prusak, and Bruno Turcksin for help with GPU support in deal.II and this tutorial.
This tutorial implements a matrix-free algorithm that works on GPUs using Kokkos to solve the Stokes problem
\begin{align} -\Delta \mathbf{u} + \nabla p &= \mathbf{f} \quad \text{in } \Omega, \\ -\nabla \cdot \mathbf{u} &= 0 \quad \text{in } \Omega, \\ \mathbf{u} &= \mathbf{0} \quad \text{on } \partial\Omega, \end{align}
where we want to solve for velocity \(\mathbf{u}\) and pressure \(p\). The Stokes equations model the behavior of slow-moving, viscous flow, which is often relevant in applications of biofluid dynamics, microfluidics, and flow of geomaterials like ice or rock.
The algorithm that is implemented in this tutorial is matrix-free, meaning we use no assembled matrices. Matrix-free methods have been instrumental in speeding up finite element computations, as memory bandwidth is often a massive bottleneck when using matrix-based methods. Matrix-free methods exploit the fact that simple computations such as addition and multiplication are relatively quick compared to memory bandwidth, so instead of storing data in a matrix and using them to apply the action of the linear operator, matrix vector products are computed on the fly. This approach is discussed in detail in step-37 using the MatrixFree framework that runs on CPUs only.
Here, we use the Portable::MatrixFree class instead, which is a similar framework but implemented using Kokkos to run on GPUs and first introduced in step-64. Just as in that tutorial, which solves the Helmholtz equation instead of the Stokes equations, we use Portable::MatrixFree and Portable::FEEvaluation, which make use of the Kokkos library to implement a portable algorithm that works well on both GPUs and CPUs. Finite element simulations on GPUs offer several benefits compared to traditional CPU calculations, for example a higher compute performance and memory bandwidth for a specified power budget.
Unlike step-64, this tutorial solves for two solution variables, namely velocity \(\mathbf{u}\) and pressure \(p\), which means we will need two DoFHandlers and Portable::FEEvaluation objects, and the solution will be stored in a block vector with two blocks. Unlike earlier tutorials that assemble matrices, like step-22, we do not build an FESystem that represents the combination of velocity and pressure. Instead, the Portable::MatrixFree object is initialized with a vector of DoFHandlers. Later, we can reference the velocity and pressure by referring to the DoFHandler index 0 (velocity) and index 1 (pressure). These indices will be passed to Portable::FEEvaluation objects when constructed.
While this program runs without a GPU, using a workstation with one or more GPUs or a GPU cluster is strongly recommended. Unfortunately, installing deal.II with GPU support is non-trivial. Please see the deal.II GPU wiki page for installation instructions.
The Stokes system and linear solvers for it have been discussed in various tutorials including step-22 (discretization, Schur complement solver), step-32 and step-56 (using the same block preconditioner), and also step-55 (simpler block preconditioner, combined with AMG). Instead of repeating this here, we refer to the derivation in these tutorials and only highlight the differences.
While the usage of stable Taylor-Hood elements is the same as in step-22 and step-56, the discretization is slightly different because we consider the standard Laplacian instead of the deformation tensor acting on the velocity. This leads to the weak form
\begin{eqnarray*} (\nabla \mathbf{u}, \nabla \mathbf{v})_{\Omega} - (p, \nabla \cdot \mathbf{v})_{\Omega} &= (\mathbf{f}, \mathbf{v})_{\Omega},\\ -(q, \nabla \cdot \mathbf{u})_{\Omega} &=&0. \end{eqnarray*}
Nevertheless, we arrive at the same block linear system
\begin{eqnarray*} \left(\begin{array}{cc} A & B^T \\ B & 0 \end{array}\right) \end{eqnarray*}
and use the same block preconditioner as in step-32 and step-56. We are left with the choice for the approximation of the inverse of the velocity or \(A\) block, the inverse of the Schur complement \(S\), the implementation of the \(B^T\) operator, and the Stokes operator, all of which we perform matrix-free.
Here, we choose to approximate the inverse of the velocity block by a single v-cycle using the matrix-free, global coarsening, geometric multigrid method. See [181] for the algorithmic details for the CPU implementation. As a smoother we use a Chebyshev iteration around the point-Jacobi scheme, i.e., the inverse of the diagonal of \(A\). This is discussed in step-37.
The inverse of the Schur complement \(S\) is approximated by a Chebyshev iteration applied to the pressure mass matrix (without using multigrid). We do not need to use a multigrid method because the mass matrix is much easier to solve, and because using multigrid on a significantly smaller problem would cause additional computational overhead.
While this is a slight variation on the preconditioners discussed in other tutorials, the same principles apply. Because we do not perform inner solves as done in step-32, the preconditioner is a linear operator. Therefore, we do not require the use of flexible GMRES. This saves precious device memory on the GPUs.
Not all parts of this example run on the GPU. Partly, this is because GPU support in deal.II is still early and some parts could and most certainly will be ported to run on GPUs. Additionally, some steps in this program do not lend themselves to GPU acceleration because they are hard or impossible to parallelize and because they require IO, complex logic, or dynamic memory allocation.
As in step-64, we use the template parameter MemorySpace::Default and MemorySpace::Host to indicate where the solution data of a particular vector is stored. MemorySpace::Host always indicates the host memory (CPU or main memory), while MemorySpace::Default refers to the memory of the default compute device (either CPUs or GPUs depending on how the Kokkos library used by deal.II was configured). While LinearAlgebra::distributed::Vector has a default value for the memory space, which is MemorySpace::Host, we always specify the memory space explicitly in this example. Any function operating on vectors stored on the compute device have to run on the GPU as a kernel and has to be marked with the DEAL_II_HOST_DEVICE precompiler macro to let the compiler know that this code has to be compiled for the programming model of the relevant compute device. This implies that the code in these functions cannot make use of functionality that is not available on this device. Before the CPU can operate on a vector, the data has to be copied into system memory (by copying into a vector with MemorySpace::Host).
In short, we currently perform setup and postprocessing on the CPU, while the linear solve runs on the GPU. The example performs the following steps:
Notice that roughly speaking the setup and postprocess happen on the CPU, with the linear solve step happening on the GPU. With complex problems, this will be the majority of the computational time.
Our test problem features a unit hypercube with the manufactured solution
\begin{align} u &= \pi \sin^2(\pi x) \sin(2\pi y), \\ v &= -\pi \sin(2\pi x) \sin^2(\pi y), \\ p &= \cos(\pi x) \cos(\pi y) \end{align}
in two dimensions and
\begin{align} u &= \pi \sin^2(\pi x)(\sin^2(\pi z)\sin(2\pi y)-\sin^2(\pi y)\sin(2\pi z)), \\ v &= \pi \sin^2(\pi y) (\sin^2(\pi x) \sin(2\pi z) - \sin^2(\pi z)\sin(2\pi x)), \\ w &= \pi \sin^2(\pi z) (\sin^2(\pi y) \sin(2\pi x) - \sin^2(\pi x)\sin(2\pi y)), \\ p &= \cos(\pi x)\cos(\pi y)\cos(\pi z) \end{align}
in three dimensions.
The index of the velocity and pressure DoFHandler into the vector of DoFHandlers inside Portable::MatrixFree and the blocks of the solution vector.
We start with the definition of the right-hand side and exact solution of the manufactured solution in 2 and 3 dimensions. In 2D the exact solution is given by
\begin{align*} u &= \pi \sin^2(\pi x) \sin(2 \pi y) \\ v &= -\pi \sin(2 \pi x) \sin^2(\pi y) \\ p &= \cos(\pi x) \cos(\pi y) \end{align*}
and in 3D by:
\begin{align*} u &= \pi \sin^2(\pi x)(\sin^2(\pi z)\sin(2\pi y)-\sin^2(\pi y)\sin(2\pi z)) \\ v &= \pi \sin^2(\pi y) (\sin^2(\pi x) \sin(2\pi z) - \sin^2(\pi z)\sin(2\pi x)) \\ w &= \pi \sin^2(\pi z) (\sin^2(\pi y) \sin(2\pi x) - \sin^2(\pi x)\sin(2\pi y)) \\ p &= \cos(\pi x)\cos(\pi y)\cos(\pi z) \end{align*}
The following classes define these in code.
The matrix-free operator for the velocity block \(A\) given by \((\nabla u,\nabla v)\) is defined by the class PortableMFVelocityOperator. It uses the class VelocityCellOperator, which is evaluated in parallel on each cell. On each cell, we define the action at each quadrature point with the small helper class VelocityOperatorQuad with operator().
This class finally provides the matrix-free operator for the velocity block. Note that we also compute the inverse diagonal of the operator, which is used in the Chebyshev smoother when we approximate A^{-1} with a GMG v-cycle. We note that Tvmult() is not implemented because it is not required for the smoother and \(A\) is symmetric anyway.
The preconditioner requires a Schur complement approximation, which is here given by a mass matrix in the pressure space \((p,q)\). This is implemented in a very similar way to the velocity block above. A notable difference is that we select the pressure by passing a dof_handler_index of 1 instead of 0 to the FEEvaluation class.
The following set of classes provides the whole Stokes operator
\begin{eqnarray*} \begin{bmatrix} A & B^T \\ B & 0 \end{bmatrix}. \end{eqnarray*}
While structured in a similar way (class PortableMFStokesOperator uses StokesCellOperator and a lambda function for the action at each quadrature point), we now operate on Portable::DeviceBlockVector and use two Portable::FEEvaluation objects, one for the velocity and one for the pressure.
Note that we don't need support for computing the diagonal, as this is not needed in the block preconditioner.
The last ingredient of the block preconditioner is the action of the block \(B^T\) given by \(-(p,\nabla \cdot v)\). The operator reads pressure values \(p\) and produces a velocity (you can think of it as a rectangular matrix block). Therefore, the implementation is similar to the Stokes operator in that we work with two Portable::FEEvaluation objects and operate on Portable::DeviceBlockVector.
Instead of copying constrained values, zero them out. The BT operator does not receive an input velocity to copy values from and zeroing out is the correct operation for boundary and hanging nodes for an update of the velocity:
BlockSchurPreconditionerThe following class implements the block preconditioner. The class takes the types of the operators for \(A^{-1}\), \(S^{-1}\), and \(B^T\) as template arguments. This is the same preconditioner used in step-32 and step-56.
We keep a temporary vector tmp inside this class to avoid reallocating memory in every vmult call. It has to be mutable because it is used in the vmult method that is declared as const.
Allocate the temporary vector on first use. Its content doesn't matter, as it will be overwritten in the vmult() call.
First apply the Schur Complement inverse operator: dst.p = S^-1 * src.p
Apply the top right block: tmp.u = -B^T * dst.p + src.u
Finally the velocity block:
StokesProblemThe remaining part of this tutorial is the StokesProblem class that puts everything together.
The setup_dofs() function distributes the two separate DoFHandlers for the velocity and pressure before initializing the MatrixFree object with an std::vector of both of them. We can later refer to the velocity using DoFHandler index 0 and the pressure using DoFHandler index 1.
create the right hand side on the host and move to device:
In the solve() function we set up the preconditioner and run the GMRES solver. For this, we construct the multigrid hierarchy for the GMG v-cycle with a Chebyshev iteration around the point-Jacobi scheme, i.e., the inverse of the diagonal of \(A\), to approximate the action of \(A^{-1}\). We approximate the Schur Complement with a Chebyshev iteration applied to the pressure mass matrix (without multigrid).
Do not go down to level 0, because this will lead to slower runtime as the problem becomes very small:
Prepare the operators and data structures on all levels of the multigrid hierarchy
On the finest level we can reuse the MatrixFree object from the Stokes operator. This way we can solve significantly larger problems before we run out of device memory.
transfer operator
smoother
Use the Chebyshev iteration as an (approximate) solver on the coarsest level. In this mode smoothing_range is a relative target tolerance and must be strictly less than one; the number of iterations is then chosen automatically by setting degree to numbers::invalid_unsigned_int. We also use more CG iterations for the eigenvalue estimate because when min_level > 0, the coarse problem can still be reasonably large and badly conditioned.
These values are chosen by experimentation for the problem at hand. We chose the smoothing range first. A good value will allow the smoother to effectively separate large and small scale oscillations in the residual and as such improve the convergence of the Chebyshev iteration and the multigrid method. Finally, the degree is chosen to minimize total runtime (a larger value increases the cost but improves the outer number of GMRES iterations).
Estimate and print the eigenvalue spectrum of the velocity block on each level. This spectrum is later used by the Chebyshev iteration.
coarse-grid solver
put everything together
The postprocess() function moves the solution to host memory and integrates the difference to the manufactured solution to compute errors.
The run() function prints some statistics and then performs a familiar refinement loop.
main() functionThe only interesting bits here are the template arguments that specify dimension and polynomial degree to be used.
When running the program on a single GPU (here an H100) and computing the L2 error norms between the finite element solution and the analytical solution discussed above, we can produce the following table of results:
| Refinement | DoFs | Iterations | Time (s) | Velocity Error | Pressure Error |
|---|---|---|---|---|---|
| 0 | 2,312 | 10 | 2.4e-02 | 4.5e-02 | 1.5e-01 |
| 1 | 15,468 | 20 | 6.0e-02 | 6.0e-03 | 1.2e-02 |
| 2 | 112,724 | 22 | 1.1e-01 | 7.6e-04 | 1.2e-03 |
| 3 | 859,812 | 23 | 2.5e-01 | 9.6e-05 | 2.3e-04 |
| 4 | 6,714,692 | 23 | 1.2e+00 | 1.2e-05 | 5.6e-05 |
| 5 | 53,070,468 | 24 | 8.8e+00 | 1.5e-06 | 1.4e-05 |
We see optimal convergence rates for L2 errors of velocity and pressure (3 and 2, respectively, which corresponds to best-approximation of the Q2 and Q1 spaces) and decent linear scaling with problem size: A factor 8 in increase in unknowns leads to a factor of 8 in solve time, at least on the finest refinement level.
GPUs profit from fine-grained parallelism. In fact, a modern GPU like an H100 requires a significant number of small parallel work items to saturate the device. An H100 has 132 SMs, each needing about 10-20 thread blocks to be occupied, with 256 threads per block. This would require at least 300,000 parallel threads to make use of the whole GPU. This is only a simplified estimate and only holds in the case the parallel task is compute-bound, while here, in fact, most kernels are typically memory-bound. Regardless, the matrix-free operator and smoother application is parallelized by a loop over cells and inner loops over 1d or 3d quadrature point loops, depending on the operation. This means that higher order methods offer more parallelism and they also need to stream in less memory per DoF.
Based on this argument, we would expect better performance for higher order methods. This is clearly visible when we take a look at the following results (run on an RTX 6000 using the finest mesh that fits into memory):
| Element | DoFs | Operator MDoFs/s |
|---|---|---|
| Q2Q1 | 53,070,468 | 1,746 |
| Q3Q2 | 188,174,468 | 1,627 |
| Q4Q3 | 58,112,836 | 2,895 |
| Q5Q4 | 116,203,076 | 3,160 |
The number of DoFs processed per second for a single application of the Stokes operator clearly increases significantly with polynomial degree. The whole solver is a bit harder to objectively compare as the higher order discretization requires more iterations and adjustments to GMG parameters like the Chebychev degree of the smoother. To keep things simple, we omit these results here.
The advantage of using Kokkos instead of writing a code in one of the vendor languages like CUDA, is that we can run on AMD, NVIDIA, Intel or other future GPU devices without any porting required. Furthermore, there are multiple host-parallel backends available in Kokkos and thus this tutorial also runs with the CPU-only configurations. This is perfect for code development and debugging, but of course significantly slower. By default, the Kokkos Serial backend does not use multiple cores or explicit SIMD instructions. As such, the CPU-based matrix-free deal.II code as in step-37 will be significantly faster.
As an example, here is a quick comparison between different workstation and server GPUs and the time to solution for the 53 million DoF problem size:
| Name | FP64 [TFLOPS] | VRAM [GB] | Bandwidth [TB/s] | Operator [MDoFs/s] | Solve Time [s] |
|---|---|---|---|---|---|
| AMD W7800 | 1.41 | 32 | 0.6 | 482 | 32 |
| NVIDIA RTX 6000 | 1.97 | 96 | 1.8 | 1,746 | 11 |
| NVIDIA H100 | 25.6 | 80 | 3.4 | 1,457 | 9 |
| 32 core 5975WX CPU | 1.4 | - | 0.05 | 72 | 220 |
The last line highlights the advantage of building on Kokkos over implementing CUDA or HIP code directly: Even without further optimization for CPUs (for example no SIMD vectorization is currently used) we can achieve decent performance even when running without a GPU. Here we are using the serial Kokkos backend and 32 MPI ranks. This is certainly fast enough for debugging and code development.
This example program can also be run on more than one GPU, which allows us to solve larger problems before running out of VRAM and to, at least ideally, achieve a faster time to solution.
The way one can use more than one GPU involves running with more than one MPI rank, like in step-40. This can be done by allocating one MPI rank per GPU. As an example, if you have single node or a workstation with 4 GPUs, you would run
This is typically enough for Kokkos to automatically pick a different GPU for each MPI rank. Kokkos has an API and command line options to control job placement manually, but this should not be necessary. Note that you can also oversubscribe a GPU cluster by running more MPI ranks than GPUs on a system. This will result in more than one rank sharing a GPU, which will typically degrade performance due to scheduling overhead and workload imbalance.
Note that because of this performance degradation the typical rule of using one MPI rank per physical CPU core is not applicable for a GPU enabled code. Part of the code that runs on the CPU (matrix-free setup, error estimation, etc.) can make use of multiple threads. This means for optimal performance you will need to allocate several threads per MPI rank if you are running on a GPU cluster. For example, on the Clemson Palmetto Cluster a job submission of
would allocate a node with 8 H100 GPUs, 8 MPI ranks, and 10 threads per rank.
As an example, here is a table with solve times running one of these nodes with up to 8 H100 GPUs:
| Refinement | DoFs | 1 GPU (s) | 2 GPU (s) | 4 GPU (s) | 8 GPU (s) |
|---|---|---|---|---|---|
| 0 | 2,312 | .065 | .074 | .078 | .077 |
| 1 | 15,468 | .16 | .21 | .21 | .21 |
| 2 | 112,724 | .49 | 1.10 | .98 | .81 |
| 3 | 859,812 | 2.9 | 4.1 | 3.4 | 3.3 |
| 4 | 6,714,692 | 25 | 30 | 18 | 14 |
| 5 | 53,070,468 | 210 | 280 | 160 | 100 |
| 6 | 421,991,684 | - | - | 1500 | 970 |
Notice that we can solve larger problems with 4 or more GPUs compared to a single GPU. On the one hand, we also have some reduction in time to solution, at least for the larger problem sizes. On the other hand, notice that time to solution does not scale linearly. This is likely due to additional communication overhead, which we will investigate and improve in the future.
You might have noticed that we used Number instead of double throughout the tutorial program. This is on purpose. Try changing the template argument to float in main(). With this change, all computations are now done in single instead of double precision. Especially on workstation GPUs, single precision computations are significantly faster (up to 64x). Even when the computations are memory-bound and not compute-bound, which is the case for most of the kernels in this example, memory usage is reduced by a factor of 2. Check what speedup you can observe and see whether you can solve bigger problems with the same amount of GPU memory.
Solving with lower precision comes with a trade-off, though. The integration of the operator and right-hand side is less accurate, which can lead to a degradation in accuracy of the solution, especially on finer meshes. By comparing the errors, determine the refinement level where floating-point error is larger than the error introduced by the discretization and linear solver tolerance.
Furthermore, be aware that our error computation using VectorTools::integrate_difference is also less precise when performed in single precision. It would be better to convert the solution to double precision and compute the error in double precision. The function LinearAlgebra::distributed::Vector::copy_locally_owned_data_from can do this conversion for you. Does that make a difference?
Finally, a more advanced strategy involves solving the linear system in double precision while running the multigrid preconditioner in single precision. This approach is called "mixed precision preconditioning". The method is very attractive on GPUs as it combines the accuracy of the double precision approach with the performance of the single precision approach. Considering we spent most of the computational effort inside the velocity block multigrid preconditioner, which does not need to be very accurate, this seems promising. The implementation is not difficult apart from having to carefully choose the correct number type in all GMG-related objects. For clarity, the mixed precision approach is not included in this tutorial, but it makes for a great exercise!
1.16.1