This tutorial depends on step-22, step-64.
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.
Introduction
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.
Weak Form
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*}
Linear system and preconditioning
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.
GPU Acceleration
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:
- Setup
- Generate and refine the mesh (CPU)
- Distribute DoFs, compute constraints (CPU)
- Assemble right-hand side vector (CPU)
- Create mesh hierarchy for multigrid (CPU)
- Setup
Portable::MatrixFree operator: evaluate inverse Jacobians, mapping, cell geometry (CPU)
- Setup multigrid operators, smoothers, transfer (CPU)
- Solve
- GMRES algorithm (CPU logic, operations on GPU)
- Vector operations: inner products, ... (GPU)
- Matrix-vector product:
Portable::MatrixFree (GPU)
- GMG preconditioner: smoother, transfer, coarse solve (GPU)
- Postprocess
- Move solution to host memory
- Compute errors (currently CPU only)
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.
Test case
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 commented program
#include <deal.II/base/conditional_ostream.h>
#include <deal.II/base/timer.h>
#include <deal.II/distributed/tria.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/lac/la_parallel_block_vector.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/solver_gmres.h>
#include <deal.II/matrix_free/operators.h>
#include <deal.II/matrix_free/portable_fe_evaluation.h>
#include <deal.II/matrix_free/portable_matrix_free.h>
#include <deal.II/multigrid/mg_coarse.h>
#include <deal.II/multigrid/mg_matrix.h>
#include <deal.II/multigrid/mg_smoother.h>
#include <deal.II/multigrid/mg_transfer_global_coarsening.h>
#include <deal.II/multigrid/mg_transfer_matrix_free.h>
#include <deal.II/multigrid/multigrid.h>
#include <deal.II/multigrid/portable_mg_transfer_global_coarsening.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/vector_tools_integrate_difference.h>
* * * struct InterferenceTaperTransform *
The index of the velocity and pressure DoFHandler into the vector of DoFHandlers inside Portable::MatrixFree and the blocks of the solution vector.
constexpr unsigned int velocity_dof_handler_index = 0;
constexpr unsigned int pressure_dof_handler_index = 1;
Problem Definition
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.
template <
int dim,
typename Number>
class VelocityRightHandSide :
public Function<dim, Number>
const unsigned int component = 0)
const override;
template <
int dim,
typename Number>
VelocityRightHandSide<dim, Number>::value(
const Point<dim> &p,
const unsigned int component)
const
const double pi2 = pi * pi;
return pi * cy * (16.0 * pi2 * sx * sx * sy - 4.0 * pi2 * sy - sx);
return pi * cx * (-16.0 * pi2 * sx * sy * sy + 4.0 * pi2 * sx - sy);
const double sx2 = sx * sx;
const double sy2 = sy * sy;
const double sz2 = sz * sz;
const double cx2 = cx * cx;
const double cy2 = cy * cy;
const double cz2 = cz * cz;
const double pi3 = pi * pi2;
return -16.0 * pi3 * sx2 * sy2 * sz * cz +
16.0 * pi3 * sx2 * sy * sz2 * cy -
4.0 * pi3 * sx2 * sy * cy * cz2 +
4.0 * pi3 * sx2 * sz * cy2 * cz - pi * sx * cy * cz +
4.0 * pi3 * sy2 * sz * cx2 * cz -
4.0 * pi3 * sy * sz2 * cx2 * cy;
return 16.0 * pi3 * sx2 * sy2 * sz * cz -
4.0 * pi3 * sx2 * sz * cy2 * cz -
16.0 * pi3 * sx * sy2 * sz2 * cx +
4.0 * pi3 * sx * sy2 * cx * cz2 +
4.0 * pi3 * sx * sz2 * cx * cy2 -
4.0 * pi3 * sy2 * sz * cx2 * cz - pi * sy * cx * cz;
return -16.0 * pi3 * sx2 * sy * sz2 * cy +
4.0 * pi3 * sx2 * sy * cy * cz2 +
16.0 * pi3 * sx * sy2 * sz2 * cx -
4.0 * pi3 * sx * sy2 * cx * cz2 -
4.0 * pi3 * sx * sz2 * cx * cy2 +
4.0 * pi3 * sy * sz2 * cx2 * cy - pi * sz * cx * cy;
template <
int dim,
typename Number>
class VelocitySolution :
public Function<dim, Number>
const unsigned int component = 0)
const override;
template <
int dim,
typename Number>
VelocitySolution<dim, Number>::value(
const Point<dim> &p,
const unsigned int component)
const
const double s2x =
std::sin(2.0 * pi * x);
const double s2y =
std::sin(2.0 * pi * y);
return pi * sx * sx * s2y;
return -pi * s2x * sy * sy;
const double s2z =
std::sin(2.0 * pi * z);
const double sx2 = sx * sx;
const double sy2 = sy * sy;
const double sz2 = sz * sz;
const double dphidx = pi * s2x * sy2 * sz2;
const double dphidy = pi * s2y * sx2 * sz2;
const double dphidz = pi * s2z * sx2 * sy2;
template <
int dim,
typename Number>
class PressureSolution :
public Function<dim, Number>
const unsigned int = 0) const override
virtual RangeNumberType value(const Point< dim > &p, const unsigned int component=0) const
#define AssertIndexRange(index, range)
* * * RotationFunction< dim, Number >::RotationFunction Number(dim)
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
The velocity operator
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().
template <
int dim,
int fe_degree,
typename Number>
class VelocityOperatorQuad
const auto gradient_u = fe_eval->get_gradient(q_point);
fe_eval->submit_gradient(gradient_u, q_point);
class VelocityCellOperator
static const unsigned int n_q_points =
::Utilities::pow(n_q_points_1d, dim);
data, velocity_dof_handler_index);
fe_u.read_dof_values(src);
VelocityOperatorQuad<dim, degree_u, Number> quad_operation;
data->for_each_quad_point(
[&](
const int q_point) { quad_operation(&fe_u, q_point); });
fe_u.distribute_local_to_global(dst);
* * Point< dim > operator()(const Point< dim > &p) const *
#define DEAL_II_HOST_DEVICE
std::vector< index_type > data
Kokkos::View< Number *, MemorySpace::Default::kokkos_space > DeviceVector
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.
typename Number = double,
int n_q_points_1d = degree_u + 1>
PortableMFVelocityOperator() =
default;
PortableMFVelocityOperator(
void initialize_dof_vector(VectorType &vec)
const
data->initialize_dof_vector(vec, velocity_dof_handler_index);
return data->get_vector_partitioner(velocity_dof_handler_index)->size();
get_matrix_diagonal_inverse() const
return inverse_diagonal_entries;
Assert(row == col, ExcNotImplemented());
Assert(inverse_diagonal_entries.get() !=
nullptr &&
inverse_diagonal_entries->m() > 0,
return 1.0 / (*inverse_diagonal_entries)(row, row);
void vmult(VectorType &dst,
const VectorType &src)
const
dst =
static_cast<Number>(0.);
VelocityCellOperator<dim, degree_u, degree_p, Number, n_q_points_1d>
data->cell_loop(velocity_operator, src, dst);
data->copy_constrained_values(src, dst, velocity_dof_handler_index);
void Tvmult(VectorType & ,
const VectorType & )
const
Assert(
data.get() !=
nullptr, ExcNotInitialized());
this->inverse_diagonal_entries =
std::make_shared<DiagonalMatrix<VectorType>>();
this->inverse_diagonal_entries->get_vector();
data->initialize_dof_vector(inverse_diagonal, velocity_dof_handler_index);
VelocityOperatorQuad<dim, degree_u, Number> velocity_operator_quad;
compute_diagonal<dim, degree_u, degree_u + 1, dim, Number>(
velocity_dof_handler_index);
Number *raw_diagonal = inverse_diagonal.get_values();
inverse_diagonal.locally_owned_size(),
ExcMessage(
"Diagonal entries of a positive definite operator "
raw_diagonal[i] = 1. / raw_diagonal[i];
std::shared_ptr<Portable::MatrixFree<dim, Number>>
data;
std::shared_ptr<DiagonalMatrix<VectorType>> inverse_diagonal_entries;
#define Assert(cond, exc)
#define AssertThrow(cond, exc)
Tpetra::Vector< Number, LO, GO, NodeType< MemorySpace > > VectorType
The Schur complement operator
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.
fe_eval->submit_value(fe_eval->get_value(q_point), q_point);
static const unsigned int n_q_points =
::Utilities::pow(n_q_points_1d, dim);
MassCellOperator<dim, degree_u, degree_p, Number, n_q_points_1d>::operator()(
data, pressure_dof_handler_index);
fe_p.read_dof_values(src);
MassOperatorQuad<dim, degree_u, degree_p, Number, n_q_points_1d>
data->for_each_quad_point(
[&](
const int &q_point) { quad_operation(&fe_p, q_point); });
fe_p.distribute_local_to_global(dst);
int n_q_points_1d = degree_u + 1>
return
data->get_vector_partitioner(pressure_dof_handler_index)->
size();
Assert(row == col, ExcNotImplemented());
Assert(inverse_diagonal_entries.get() !=
nullptr &&
inverse_diagonal_entries->m() > 0,
return 1.0 / (*inverse_diagonal_entries)(row, row);
void vmult(VectorType &dst,
const VectorType &src)
const
dst =
static_cast<Number>(0.);
MassCellOperator<dim, degree_u, degree_p, Number, n_q_points_1d>
data->cell_loop(mass_operator, src, dst);
data->copy_constrained_values(src, dst, pressure_dof_handler_index);
std::shared_ptr<DiagonalMatrix<VectorType>>
get_matrix_diagonal_inverse() const
return inverse_diagonal_entries;
this->inverse_diagonal_entries =
std::make_shared<DiagonalMatrix<VectorType>>();
this->inverse_diagonal_entries->get_vector();
MassOperatorQuad<dim, degree_u, degree_p, Number, n_q_points_1d>
compute_diagonal<dim, degree_p, n_q_points_1d, 1, Number>(
pressure_dof_handler_index);
Number *raw_diagonal = inverse_diagonal.get_values();
inverse_diagonal.locally_owned_size(),
ExcMessage(
"Diagonal entries of a positive definite operator "
raw_diagonal[i] = 1. / raw_diagonal[i];
std::shared_ptr<Portable::MatrixFree<dim, Number>>
data;
std::shared_ptr<DiagonalMatrix<VectorType>> inverse_diagonal_entries;
* * * * std::vector< Number > ThermoPlasticMaterial< dim, ViscoplasticYieldLaw, Number >::get_state_parameters const
unsigned int global_dof_index
The Stokes operator
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.
static const unsigned int n_q_points =
::Utilities::pow(n_q_points_1d, dim);
data, velocity_dof_handler_index);
data, pressure_dof_handler_index);
fe_u.read_dof_values(src.block(0));
fe_p.read_dof_values(src.block(1));
data->for_each_quad_point([&](
const int &q_point) {
const Number pressure_value = fe_p.get_value(q_point);
for (
unsigned int d = 0; d < dim; ++d)
velocity_term[d][d] -= pressure_value;
fe_u.submit_gradient(velocity_term, q_point);
const Number pressure_term = -
trace(gradient_u);
fe_p.submit_value(pressure_term, q_point);
fe_u.distribute_local_to_global(dst.block(0));
fe_p.distribute_local_to_global(dst.block(1));
int n_q_points_1d = degree_u + 1>
class PortableMFStokesOperator
PortableMFStokesOperator(
dst = static_cast<
Number>(0.);
StokesCellOperator<dim, degree_u, degree_p, Number, n_q_points_1d>
data->cell_loop(stokes_operator, src, dst);
data->copy_constrained_values(src, dst);
std::shared_ptr<Portable::MatrixFree<dim, Number>>
data;
constexpr Number trace(const SymmetricTensor< 2, dim2, Number > &)
The BT operator
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.
static const unsigned int n_q_points =
::Utilities::pow(n_q_points_1d, dim);
data, velocity_dof_handler_index);
data, pressure_dof_handler_index);
fe_p.read_dof_values(src.block(1));
data->for_each_quad_point([&](
const int &q_point) {
const Number pressure_value = fe_p.get_value(q_point);
fe_u.submit_divergence(-pressure_value, q_point);
fe_u.distribute_local_to_global(dst.block(0));
int n_q_points_1d = degree_u + 1>
class PortableMFBTOperator
dst = static_cast<
Number>(0.);
BTCellOperator<dim, degree_u, degree_p, Number, n_q_points_1d>
data->cell_loop(cell_operator, src, dst);
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:
data->set_constrained_values(0.0,
velocity_dof_handler_index);
std::shared_ptr<Portable::MatrixFree<dim, Number>>
data;
The Preconditioner BlockSchurPreconditioner
The 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.
template <
class AInvOperator,
BlockSchurPreconditioner(
const AInvOperator &A_inverse_operator,
const SInvOperator &S_inverse_operator,
const BTOperator &BT_operator);
void vmult(VectorType &dst,
const VectorType &src)
const;
const AInvOperator &A_inverse_operator;
const SInvOperator &S_inverse_operator;
const BTOperator &BT_operator;
template <
class AInvOperator,
BlockSchurPreconditioner<AInvOperator, SInvOperator, BTOperator, VectorType>
::
BlockSchurPreconditioner(
const AInvOperator &A_inverse_operator,
const SInvOperator &S_inverse_operator,
const BTOperator &BT_operator)
: A_inverse_operator(A_inverse_operator)
, S_inverse_operator(S_inverse_operator)
, BT_operator(BT_operator)
template <class AInvOperator,
BlockSchurPreconditioner<AInvOperator, SInvOperator, BTOperator, VectorType>::
vmult(VectorType &dst, const VectorType &src) 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
S_inverse_operator.vmult(dst.block(1), src.block(1));
Apply the top right block: tmp.u = -B^T * dst.p + src.u
BT_operator.vmult(tmp, dst);
tmp.block(0).sadd(-1.0, 1.0, src.block(0));
Finally the velocity block:
A_inverse_operator.vmult(dst.block(0), tmp.block(0));
The main class StokesProblem
The remaining part of this tutorial is the StokesProblem class that puts everything together.
template <
int dim,
int degree_p,
typename Number =
double>
static constexpr unsigned int degree_u = degree_p + 1;
std::shared_ptr<Portable::MatrixFree<dim, Number>> mf_data;
BlockVectorType solution;
template <
int dim,
int degree_p,
typename Number>
StokesProblem<dim, degree_p, Number>::StokesProblem()
, fe_u(
FE_Q<dim>(degree_p + 1), dim)
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.
template <
int dim,
int degree_p,
typename Number>
void StokesProblem<dim, degree_p, Number>::setup_dofs()
dof_u.distribute_dofs(fe_u);
dof_p.distribute_dofs(fe_p);
const IndexSet &owned_set_u = dof_u.locally_owned_dofs();
constraints_u.reinit(owned_set_u, relevant_set_u);
const IndexSet &owned_set_p = dof_p.locally_owned_dofs();
constraints_p.reinit(owned_set_p, relevant_set_p);
std::vector<const DoFHandler<dim> *> dof_handlers = {&dof_u, &dof_p};
std::vector<const AffineConstraints<Number> *> constraints = {
&constraints_u, &constraints_p};
mf_data = std::make_shared<Portable::MatrixFree<dim, Number>>();
mf_data->reinit(mapping, dof_handlers, constraints, quad, additional_data);
void make_hanging_node_constraints(const DoFHandler< dim, spacedim > &dof_handler, AffineConstraints< number > &constraints)
@ update_values
Shape function values.
@ update_gradients
Shape function gradients.
UpdateFlags mapping_update_flags
create the right hand side on the host and move to device:
mf_data->initialize_dof_vector(rhs_host);
VelocityRightHandSide<dim, Number>(),
mf_data->initialize_dof_vector(rhs);
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).
template <
int dim,
int degree_p,
typename Number>
void StokesProblem<dim, degree_p, Number>::solve()
PortableMFStokesOperator<dim, degree_u, degree_p, Number> stokes_operator(
mf_data->initialize_dof_vector(solution);
::Timer t(tria.get_mpi_communicator());
stokes_operator.vmult(solution, rhs);
const double time = t.wall_time();
const double mdofs_p_second =
1e-6 *
static_cast<double>(solution.size()) / time;
pcout <<
"Stokes operator: " << time <<
" s, MDoFs/s: " << mdofs_p_second
PortableMFVelocityOperator<dim, degree_u, degree_p, Number>;
SmootherPreconditionerType>;
const auto coarse_grid_triangulations =
const unsigned int max_level = coarse_grid_triangulations.size() - 1;
Do not go down to level 0, because this will lead to slower runtime as the problem becomes very small:
const unsigned int min_level =
std::min(3U, max_level - 1);
std::vector<std::shared_ptr<Portable::MatrixFree<dim, Number>>>
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
Prepare the operators and data structures on all levels of the multigrid hierarchy
auto &dof_handler = mg_dof_handlers[
level];
auto &constraint = mg_constraints[
level];
dof_handler.reinit(*coarse_grid_triangulations[
level]);
dof_handler.distribute_dofs(fe_u);
constraint.reinit(dof_handler.locally_owned_dofs(),
void make_zero_boundary_constraints(const DoFHandler< dim, spacedim > &dof, const types::boundary_id boundary_id, AffineConstraints< number > &zero_boundary_constraints, const ComponentMask &component_mask={})
@ update_JxW_values
Transformed quadrature weights.
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.
mf_data_levels.emplace_back(mf_data);
mf_data_levels.emplace_back(
mf_data_levels.back()->reinit(
mapping, dof_handler, constraint, quad, additional_data);
mg_matrices[
level].reinit(mf_data_levels.back());
transfer operator
mg_transfers[
level + 1].reinit_geometric_transfer(
mg_dof_handlers[
level + 1],
mg_constraints[
level + 1],
MGTransferType mg_transfer(mg_transfers, [&](
const auto l,
auto &vec) {
mg_matrices[l].initialize_dof_vector(vec);
smoother
mg_matrices[
level].compute_diagonal();
smoother_data[
level].preconditioner =
std::make_shared<SmootherPreconditionerType>(
*mg_matrices[
level].get_matrix_diagonal_inverse());
smoother_data[
level].constraints.copy_from(mg_constraints[
level]);
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.
smoother_data[
level].smoothing_range = 1e-3;
smoother_data[
level].eig_cg_n_iterations = 40;
constexpr unsigned int invalid_unsigned_int
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).
smoother_data[
level].smoothing_range = 5;
smoother_data[
level].degree = 4;
smoother_data[
level].eig_cg_n_iterations = 20;
mg_smoother.
initialize(mg_matrices, smoother_data);
void initialize(const MGLevelObject< MatrixType2 > &matrices, const typename RelaxationType::AdditionalData &additional_data=typename RelaxationType::AdditionalData())
Estimate and print the eigenvalue spectrum of the velocity block on each level. This spectrum is later used by the Chebyshev iteration.
pcout <<
"GMG velocity block smoothers:" << std::endl;
mg_matrices[
level].initialize_dof_vector(vec);
mg_smoother.smoothers[
level].estimate_eigenvalues(vec);
pcout <<
" level: " <<
level <<
" n_dofs: " << vec.size()
<<
", eigenvalue spectrum: [ "
<< eigenvalue_info.min_eigenvalue_estimate <<
", "
<< eigenvalue_info.max_eigenvalue_estimate <<
" ]" << std::endl;
coarse-grid solver
void initialize(const MGSmootherBase< VectorType > &coarse_smooth)
put everything together
::Timer timer_smoother;
::Timer timer_transfer;
::Timer timer_coarse;
::Timer timer_residual;
auto make_timer_lambda = [&](
::Timer &timer) {
return [&](
const bool before,
const unsigned int ) {
mg.connect_pre_smoother_step(make_timer_lambda(timer_smoother));
mg.connect_post_smoother_step(make_timer_lambda(timer_smoother));
mg.connect_residual_step(make_timer_lambda(timer_residual));
mg.connect_restriction(make_timer_lambda(timer_transfer));
mg.connect_prolongation(make_timer_lambda(timer_transfer));
mg.connect_coarse_solve(make_timer_lambda(timer_coarse));
APreconditionerType preconditioner_A(dof_u,
mg, mg_transfer);
PortableMFMassOperator<dim, degree_u, degree_p, Number> mass_operator(
mass_operator.compute_diagonal();
PortableMFMassOperator<dim, degree_u, degree_p, Number>,
SPreconditionerType preconditioner_schur;
typename SPreconditionerType::AdditionalData additional_data;
additional_data.smoothing_range = 15.;
additional_data.degree = 3;
additional_data.eig_cg_n_iterations = 10;
additional_data.constraints.copy_from(constraints_p);
additional_data.preconditioner =
mass_operator.get_matrix_diagonal_inverse();
preconditioner_schur.initialize(mass_operator, additional_data);
PortableMFBTOperator<dim, degree_u, degree_p, Number>;
BTOperatorType BT_operator(mf_data);
BlockSchurPreconditioner<APreconditionerType,
preconditioner(preconditioner_A, preconditioner_schur, BT_operator);
::Timer t(tria.get_mpi_communicator());
solver.solve(stokes_operator, solution, rhs, preconditioner);
pcout <<
"Solver converged in " << solver_control.last_step()
<<
" iterations in " << t.wall_time() <<
" seconds" << std::endl;
pcout <<
"Velocity block GMG timings:"
<<
"\n smoother: " << timer_smoother.wall_time()
<<
" s\n transfer: " << timer_transfer.wall_time()
<<
" s\n coarse : " << timer_coarse.wall_time()
<<
" s\n residual: " << timer_residual.wall_time() <<
" s"
The postprocess() function moves the solution to host memory and integrates the difference to the manufactured solution to compute errors.
template <
int dim,
int degree_p,
typename Number>
void StokesProblem<dim, degree_p, Number>::postprocess()
mf_data->initialize_dof_vector(solution_host);
solution_host.
block(0).import_elements(solution.block(0),
solution_host.block(1).import_elements(solution.block(1),
constraints_u.distribute(solution_host.block(0));
constraints_p.distribute(solution_host.block(1));
solution_host.update_ghost_values();
dof_p,
QGauss<dim>(degree_p + 2), solution_host.block(1), 0);
solution_host.block(1).add(-mean_pressure);
VelocitySolution<dim, Number>(),
PressureSolution<dim, Number>(),
pcout <<
"velocity error: " << u_l2 <<
" pressure error: " << p_l2
BlockType & block(const unsigned int i)
The run() function prints some statistics and then performs a familiar refinement loop.
template <
int dim,
int degree_p,
typename Number>
void StokesProblem<dim, degree_p, Number>::run()
pcout << std::setprecision(10);
pcout << "\nKokkos execution space: "
<<
Kokkos::DefaultExecutionSpace::name();
<< "dim: " << dim << '\n'
<< "Element: Q" << degree_u << "-Q" << degree_p <<
std::endl;
unsigned int n_refinements = 10;
for (
unsigned int i = 0; i < n_refinements; ++i)
pcout <<
"\nrefinement: " << i
<<
", n_dofs: " << dof_u.n_dofs() + dof_p.n_dofs() <<
" = "
<< dof_u.n_dofs() <<
" + " << dof_p.n_dofs() << std::endl;
* * for(const auto &cell :triangulation.active_cell_iterators())
static unsigned int n_threads()
constexpr bool running_in_debug_mode()
void hyper_cube(Triangulation< dim, spacedim > &tria, const double left=0., const double right=1., const bool colorize=false)
unsigned int n_mpi_processes(const MPI_Comm mpi_communicator)
The main() function
The only interesting bits here are the template arguments that specify dimension and polynomial degree to be used.
int main(
int argc,
char **argv)
{
const unsigned int dim = 3;
const unsigned int degree_p = 1;
StokesProblem<dim, degree_p, double> problem;
* * int main(int argc, char **argv)
Results
Validation
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.
Higher Order Throughput
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.
Performance Portability
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.
Multiple GPUs
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
salloc --gpus h100:8 --ntasks 8 --cpus-per-task 10 --mem-per-cpu 2gb --time 4:00:00
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.
Possibilities for extensions
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!
The plain program
namespace Step104
{
constexpr unsigned int velocity_dof_handler_index = 0;
constexpr unsigned int pressure_dof_handler_index = 1;
template <int dim, typename Number>
class VelocityRightHandSide :
public Function<dim, Number>
{
public:
VelocityRightHandSide()
{}
const unsigned int component = 0) const override;
};
template <int dim, typename Number>
VelocityRightHandSide<dim, Number>::value(
const Point<dim> &p,
const unsigned int component) const
{
const double x = p[0];
const double y = p[1];
const double pi2 = pi * pi;
if constexpr (dim == 2)
{
if (component == 0)
return pi * cy * (16.0 * pi2 * sx * sx * sy - 4.0 * pi2 * sy - sx);
else
return pi * cx * (-16.0 * pi2 * sx * sy * sy + 4.0 * pi2 * sx - sy);
}
else
{
const double z = p[2];
const double sx2 = sx * sx;
const double sy2 = sy * sy;
const double sz2 = sz * sz;
const double cx2 = cx * cx;
const double cy2 = cy * cy;
const double cz2 = cz * cz;
const double pi3 = pi * pi2;
if (component == 0)
return -16.0 * pi3 * sx2 * sy2 * sz * cz +
16.0 * pi3 * sx2 * sy * sz2 * cy -
4.0 * pi3 * sx2 * sy * cy * cz2 +
4.0 * pi3 * sx2 * sz * cy2 * cz - pi * sx * cy * cz +
4.0 * pi3 * sy2 * sz * cx2 * cz -
4.0 * pi3 * sy * sz2 * cx2 * cy;
else if (component == 1)
return 16.0 * pi3 * sx2 * sy2 * sz * cz -
4.0 * pi3 * sx2 * sz * cy2 * cz -
16.0 * pi3 * sx * sy2 * sz2 * cx +
4.0 * pi3 * sx * sy2 * cx * cz2 +
4.0 * pi3 * sx * sz2 * cx * cy2 -
4.0 * pi3 * sy2 * sz * cx2 * cz - pi * sy * cx * cz;
else
return -16.0 * pi3 * sx2 * sy * sz2 * cy +
4.0 * pi3 * sx2 * sy * cy * cz2 +
16.0 * pi3 * sx * sy2 * sz2 * cx -
4.0 * pi3 * sx * sy2 * cx * cz2 -
4.0 * pi3 * sx * sz2 * cx * cy2 +
4.0 * pi3 * sy * sz2 * cx2 * cy - pi * sz * cx * cy;
}
}
template <int dim, typename Number>
class VelocitySolution :
public Function<dim, Number>
{
public:
VelocitySolution()
{}
const unsigned int component = 0) const override;
};
template <int dim, typename Number>
VelocitySolution<dim, Number>::value(
const Point<dim> &p,
const unsigned int component) const
{
const double x = p[0];
const double y = p[1];
const double s2x =
std::sin(2.0 * pi * x);
const double s2y =
std::sin(2.0 * pi * y);
if constexpr (dim == 2)
{
if (component == 0)
return pi * sx * sx * s2y;
else
return -pi * s2x * sy * sy;
}
else
{
const double z = p[2];
const double s2z =
std::sin(2.0 * pi * z);
const double sx2 = sx * sx;
const double sy2 = sy * sy;
const double sz2 = sz * sz;
const double dphidx = pi * s2x * sy2 * sz2;
const double dphidy = pi * s2y * sx2 * sz2;
const double dphidz = pi * s2z * sx2 * sy2;
if (component == 0)
return dphidy - dphidz;
else if (component == 1)
return dphidz - dphidx;
else
return dphidx - dphidy;
}
}
template <int dim, typename Number>
class PressureSolution :
public Function<dim, Number>
{
public:
PressureSolution()
{}
const unsigned int = 0) const override
{
if constexpr (dim == 2)
else
}
};
template <int dim, int fe_degree, typename Number>
class VelocityOperatorQuad
{
public:
*fe_eval,
const int q_point) const
{
}
};
template <int dim,
int degree_u,
int degree_p,
int n_q_points_1d>
class VelocityCellOperator
{
public:
static const unsigned int n_q_points =
{
data, velocity_dof_handler_index);
fe_u.read_dof_values(src);
VelocityOperatorQuad<dim, degree_u, Number> quad_operation;
data->for_each_quad_point(
[&](const int q_point) { quad_operation(&fe_u, q_point); });
fe_u.distribute_local_to_global(dst);
}
};
template <int dim,
int degree_u,
int degree_p,
int n_q_points_1d = degree_u + 1>
{
public:
PortableMFVelocityOperator() = default;
PortableMFVelocityOperator(
{}
{
}
void initialize_dof_vector(VectorType &vec) const
{
data->initialize_dof_vector(vec, velocity_dof_handler_index);
}
{
return data->get_vector_partitioner(velocity_dof_handler_index)->size();
}
get_matrix_diagonal_inverse() const
{
return inverse_diagonal_entries;
}
{
(void)col;
Assert(row == col, ExcNotImplemented());
Assert(inverse_diagonal_entries.get() !=
nullptr &&
inverse_diagonal_entries->m() > 0,
ExcNotInitialized());
return 1.0 / (*inverse_diagonal_entries)(row, row);
}
void vmult(VectorType &dst, const VectorType &src) const
{
dst =
static_cast<Number>(0.);
VelocityCellOperator<dim, degree_u, degree_p, Number, n_q_points_1d>
velocity_operator;
data->cell_loop(velocity_operator, src, dst);
data->copy_constrained_values(src, dst, velocity_dof_handler_index);
}
void Tvmult(VectorType & , const VectorType & ) const
{
}
{
Assert(
data.get() !=
nullptr, ExcNotInitialized());
this->inverse_diagonal_entries =
std::make_shared<DiagonalMatrix<VectorType>>();
this->inverse_diagonal_entries->get_vector();
data->initialize_dof_vector(inverse_diagonal, velocity_dof_handler_index);
VelocityOperatorQuad<dim, degree_u, Number> velocity_operator_quad;
MatrixFreeTools::
compute_diagonal<dim, degree_u, degree_u + 1, dim, Number>(
inverse_diagonal,
velocity_operator_quad,
velocity_dof_handler_index);
Number *raw_diagonal = inverse_diagonal.get_values();
Kokkos::parallel_for(
"invert A diagonal",
inverse_diagonal.locally_owned_size(),
KOKKOS_LAMBDA(int i) {
ExcMessage("Diagonal entries of a positive definite operator "
"should be positive"));
raw_diagonal[i] = 1. / raw_diagonal[i];
});
}
private:
std::shared_ptr<Portable::MatrixFree<dim, Number>>
data;
std::shared_ptr<DiagonalMatrix<VectorType>> inverse_diagonal_entries;
};
template <int dim,
int degree_u,
int degree_p,
int n_q_points_1d>
class MassOperatorQuad
{
public:
const int q_point) const
{
}
};
template <int dim,
int degree_u,
int degree_p,
int n_q_points_1d>
class MassCellOperator
{
public:
static const unsigned int n_q_points =
};
template <int dim,
int degree_u,
int degree_p,
int n_q_points_1d>
MassCellOperator<dim, degree_u, degree_p, Number, n_q_points_1d>::operator()(
{
data, pressure_dof_handler_index);
fe_p.read_dof_values(src);
MassOperatorQuad<dim, degree_u, degree_p, Number, n_q_points_1d>
quad_operation;
data->for_each_quad_point(
[&](const int &q_point) { quad_operation(&fe_p, q_point); });
fe_p.distribute_local_to_global(dst);
}
template <int dim,
int degree_u,
int degree_p,
int n_q_points_1d = degree_u + 1>
{
public:
PortableMFMassOperator(
{}
{
return data->get_vector_partitioner(pressure_dof_handler_index)->size();
}
{
(void)col;
Assert(row == col, ExcNotImplemented());
Assert(inverse_diagonal_entries.get() !=
nullptr &&
inverse_diagonal_entries->m() > 0,
ExcNotInitialized());
return 1.0 / (*inverse_diagonal_entries)(row, row);
}
void vmult(VectorType &dst, const VectorType &src) const
{
dst =
static_cast<Number>(0.);
MassCellOperator<dim, degree_u, degree_p, Number, n_q_points_1d>
mass_operator;
data->cell_loop(mass_operator, src, dst);
data->copy_constrained_values(src, dst, pressure_dof_handler_index);
}
std::shared_ptr<DiagonalMatrix<VectorType>>
get_matrix_diagonal_inverse() const
{
return inverse_diagonal_entries;
}
{
this->inverse_diagonal_entries =
std::make_shared<DiagonalMatrix<VectorType>>();
this->inverse_diagonal_entries->get_vector();
MassOperatorQuad<dim, degree_u, degree_p, Number, n_q_points_1d>
quad_operation;
MatrixFreeTools::
compute_diagonal<dim, degree_p, n_q_points_1d, 1, Number>(
inverse_diagonal,
quad_operation,
pressure_dof_handler_index);
Number *raw_diagonal = inverse_diagonal.get_values();
Kokkos::parallel_for(
"invert Mass diagonal",
inverse_diagonal.locally_owned_size(),
KOKKOS_LAMBDA(int i) {
ExcMessage("Diagonal entries of a positive definite operator "
"should be positive"));
raw_diagonal[i] = 1. / raw_diagonal[i];
});
}
private:
std::shared_ptr<Portable::MatrixFree<dim, Number>>
data;
std::shared_ptr<DiagonalMatrix<VectorType>> inverse_diagonal_entries;
};
template <int dim,
int degree_u,
int degree_p,
int n_q_points_1d>
class StokesCellOperator
{
public:
static const unsigned int n_q_points =
{
data, velocity_dof_handler_index);
data, pressure_dof_handler_index);
fe_u.read_dof_values(src.
block(0));
fe_p.read_dof_values(src.
block(1));
data->for_each_quad_point([&](
const int &q_point) {
const Number pressure_value = fe_p.get_value(q_point);
for (
unsigned int d = 0;
d < dim; ++
d)
velocity_term[d][d] -= pressure_value;
fe_u.submit_gradient(velocity_term, q_point);
fe_p.submit_value(pressure_term, q_point);
});
fe_u.distribute_local_to_global(dst.
block(0));
fe_p.distribute_local_to_global(dst.
block(1));
}
};
template <
int dim,
int degree_u,
int degree_p,
int n_q_points_1d = degree_u + 1>
class PortableMFStokesOperator
{
public:
PortableMFStokesOperator(
{}
void vmult(VectorType &dst, const VectorType &src) const
{
dst =
static_cast<Number>(0.);
StokesCellOperator<dim, degree_u, degree_p, Number, n_q_points_1d>
stokes_operator;
data->cell_loop(stokes_operator, src, dst);
data->copy_constrained_values(src, dst);
}
private:
std::shared_ptr<Portable::MatrixFree<dim, Number>>
data;
};
template <int dim,
int degree_u,
int degree_p,
int n_q_points_1d>
class BTCellOperator
{
public:
static const unsigned int n_q_points =
{
data, velocity_dof_handler_index);
data, pressure_dof_handler_index);
fe_p.read_dof_values(src.
block(1));
data->for_each_quad_point([&](
const int &q_point) {
const Number pressure_value = fe_p.get_value(q_point);
fe_u.submit_divergence(-pressure_value, q_point);
});
fe_u.distribute_local_to_global(dst.
block(0));
}
};
template <
int dim,
int degree_u,
int degree_p,
int n_q_points_1d = degree_u + 1>
class PortableMFBTOperator
{
public:
PortableMFBTOperator(
{}
void vmult(VectorType &dst, const VectorType &src) const
{
dst =
static_cast<Number>(0.);
BTCellOperator<dim, degree_u, degree_p, Number, n_q_points_1d>
cell_operator;
data->cell_loop(cell_operator, src, dst);
data->set_constrained_values(0.0,
dst.block(0),
velocity_dof_handler_index);
}
private:
std::shared_ptr<Portable::MatrixFree<dim, Number>>
data;
};
template <class AInvOperator,
class SInvOperator,
class BTOperator,
{
public:
BlockSchurPreconditioner(const AInvOperator &A_inverse_operator,
const SInvOperator &S_inverse_operator,
const BTOperator &BT_operator);
void vmult(VectorType &dst, const VectorType &src) const;
private:
const AInvOperator &A_inverse_operator;
const SInvOperator &S_inverse_operator;
const BTOperator &BT_operator;
};
template <class AInvOperator,
class SInvOperator,
class BTOperator,
BlockSchurPreconditioner<AInvOperator, SInvOperator, BTOperator, VectorType>::
BlockSchurPreconditioner(const AInvOperator &A_inverse_operator,
const SInvOperator &S_inverse_operator,
const BTOperator &BT_operator)
: A_inverse_operator(A_inverse_operator)
, S_inverse_operator(S_inverse_operator)
, BT_operator(BT_operator)
{}
template <class AInvOperator,
class SInvOperator,
class BTOperator,
void
BlockSchurPreconditioner<AInvOperator, SInvOperator, BTOperator, VectorType>::
vmult(VectorType &dst, const VectorType &src) const
{
if (tmp.size() == 0)
tmp.reinit(src);
{
S_inverse_operator.vmult(dst.block(1), src.block(1));
dst.block(1) *= -1.0;
}
{
BT_operator.vmult(tmp, dst);
tmp.block(0).sadd(-1.0, 1.0, src.block(0));
}
A_inverse_operator.vmult(dst.block(0), tmp.block(0));
}
template <int dim, int degree_p, typename Number = double>
class StokesProblem
{
public:
static constexpr unsigned int degree_u = degree_p + 1;
StokesProblem();
using BlockVectorType =
private:
void setup_dofs();
void solve();
void postprocess();
std::shared_ptr<Portable::MatrixFree<dim, Number>> mf_data;
BlockVectorType solution;
BlockVectorType rhs;
};
template <int dim, int degree_p, typename Number>
StokesProblem<dim, degree_p, Number>::StokesProblem()
: tria(MPI_COMM_WORLD)
, mapping(1)
, fe_u(
FE_Q<dim>(degree_p + 1), dim)
, fe_p(degree_p)
, dof_u(tria)
, dof_p(tria)
{}
template <int dim, int degree_p, typename Number>
void StokesProblem<dim, degree_p, Number>::setup_dofs()
{
dof_u.distribute_dofs(fe_u);
dof_p.distribute_dofs(fe_p);
const IndexSet &owned_set_u = dof_u.locally_owned_dofs();
constraints_u.reinit(owned_set_u, relevant_set_u);
constraints_u.close();
const IndexSet &owned_set_p = dof_p.locally_owned_dofs();
constraints_p.reinit(owned_set_p, relevant_set_p);
constraints_p.close();
std::vector<const DoFHandler<dim> *> dof_handlers = {&dof_u, &dof_p};
std::vector<const AffineConstraints<Number> *> constraints = {
&constraints_u, &constraints_p};
mf_data = std::make_shared<Portable::MatrixFree<dim, Number>>();
mf_data->reinit(mapping, dof_handlers, constraints, quad, additional_data);
{
rhs_host;
mf_data->initialize_dof_vector(rhs_host);
dof_u,
VelocityRightHandSide<dim, Number>(),
constraints_u);
mf_data->initialize_dof_vector(rhs);
}
}
template <int dim, int degree_p, typename Number>
void StokesProblem<dim, degree_p, Number>::solve()
{
PortableMFStokesOperator<dim, degree_u, degree_p, Number> stokes_operator(
mf_data);
mf_data->initialize_dof_vector(solution);
{
::Timer t(tria.get_mpi_communicator());
stokes_operator.vmult(solution, rhs);
const double time = t.wall_time();
const double mdofs_p_second =
1
e-6 *
static_cast<double>(solution.size()) / time;
pcout << "Stokes operator: " << time << " s, MDoFs/s: " << mdofs_p_second
<< std::endl;
solution = 0.0;
}
solver_control,
using LevelMatrixType =
PortableMFVelocityOperator<dim, degree_u, degree_p, Number>;
SmootherPreconditionerType>;
using MGTransferType =
const auto coarse_grid_triangulations =
tria);
const unsigned int max_level = coarse_grid_triangulations.size() - 1;
const unsigned int min_level =
std::min(3U, max_level - 1);
max_level);
min_level, max_level);
std::vector<std::shared_ptr<Portable::MatrixFree<dim, Number>>>
mf_data_levels;
{
auto &dof_handler = mg_dof_handlers[
level];
auto &constraint = mg_constraints[
level];
dof_handler.reinit(*coarse_grid_triangulations[
level]);
dof_handler.distribute_dofs(fe_u);
constraint.reinit(dof_handler.locally_owned_dofs(),
constraint.close();
additional_data;
mf_data_levels.emplace_back(mf_data);
else
{
mf_data_levels.emplace_back(
mf_data_levels.back()->reinit(
mapping, dof_handler, constraint, quad, additional_data);
}
mg_matrices[
level].reinit(mf_data_levels.back());
}
mg_transfers[
level + 1].reinit_geometric_transfer(
mg_dof_handlers[
level + 1],
mg_constraints[
level + 1],
MGTransferType mg_transfer(mg_transfers, [&](const auto l, auto &vec) {
mg_matrices[
l].initialize_dof_vector(vec);
});
min_level, max_level);
{
mg_matrices[
level].compute_diagonal();
smoother_data[
level].preconditioner =
std::make_shared<SmootherPreconditionerType>(
*mg_matrices[
level].get_matrix_diagonal_inverse());
smoother_data[
level].constraints.copy_from(mg_constraints[
level]);
{
smoother_data[
level].smoothing_range = 1
e-3;
smoother_data[
level].eig_cg_n_iterations = 40;
}
else
{
smoother_data[
level].smoothing_range = 5;
smoother_data[
level].degree = 4;
smoother_data[
level].eig_cg_n_iterations = 20;
}
}
mg_smoother.
initialize(mg_matrices, smoother_data);
pcout << "GMG velocity block smoothers:" << std::endl;
{
mg_matrices[
level].initialize_dof_vector(vec);
auto eigenvalue_info =
pcout <<
" level: " <<
level <<
" n_dofs: " << vec.size()
<< ", eigenvalue spectrum: [ "
<< eigenvalue_info.min_eigenvalue_estimate << ", "
<< eigenvalue_info.max_eigenvalue_estimate << " ]" << std::endl;
}
mg_coarse,
mg_transfer,
mg_smoother,
mg_smoother,
min_level,
max_level);
{
auto make_timer_lambda = [&](
::Timer &timer) {
return [&](const bool before, const unsigned int ) {
if (before)
timer.start();
else
timer.stop();
};
};
mg.connect_pre_smoother_step(make_timer_lambda(timer_smoother));
mg.connect_post_smoother_step(make_timer_lambda(timer_smoother));
mg.connect_residual_step(make_timer_lambda(timer_residual));
mg.connect_restriction(make_timer_lambda(timer_transfer));
mg.connect_prolongation(make_timer_lambda(timer_transfer));
mg.connect_coarse_solve(make_timer_lambda(timer_coarse));
}
APreconditionerType preconditioner_A(dof_u,
mg, mg_transfer);
PortableMFMassOperator<dim, degree_u, degree_p, Number> mass_operator(
mf_data);
mass_operator.compute_diagonal();
PortableMFMassOperator<dim, degree_u, degree_p, Number>,
SPreconditionerType preconditioner_schur;
{
typename SPreconditionerType::AdditionalData additional_data;
additional_data.smoothing_range = 15.;
additional_data.degree = 3;
additional_data.eig_cg_n_iterations = 10;
additional_data.constraints.copy_from(constraints_p);
additional_data.preconditioner =
mass_operator.get_matrix_diagonal_inverse();
preconditioner_schur.initialize(mass_operator, additional_data);
}
using BTOperatorType =
PortableMFBTOperator<dim, degree_u, degree_p, Number>;
BTOperatorType BT_operator(mf_data);
BlockSchurPreconditioner<APreconditionerType,
SPreconditionerType,
BTOperatorType,
BlockVectorType>
preconditioner(preconditioner_A, preconditioner_schur, BT_operator);
::Timer t(tria.get_mpi_communicator());
solver.solve(stokes_operator, solution, rhs, preconditioner);
t.stop();
pcout << "Solver converged in " << solver_control.last_step()
<< " iterations in " << t.wall_time() << " seconds" << std::endl;
pcout << "Velocity block GMG timings:"
<<
"\n smoother: " << timer_smoother.
wall_time()
<<
" s\n transfer: " << timer_transfer.
wall_time()
<<
" s\n coarse : " << timer_coarse.
wall_time()
<<
" s\n residual: " << timer_residual.
wall_time() <<
" s" << std::endl;
}
template <int dim, int degree_p, typename Number>
void StokesProblem<dim, degree_p, Number>::postprocess()
{
solution_host;
mf_data->initialize_dof_vector(solution_host);
solution_host.
block(0).import_elements(solution.block(0),
solution_host.
block(1).import_elements(solution.block(1),
constraints_u.distribute(solution_host.
block(0));
constraints_p.distribute(solution_host.
block(1));
solution_host.
block(1).add(-mean_pressure);
VelocitySolution<dim, Number>(),
cellwise_errors_ul2,
quadrature_formula,
PressureSolution<dim, Number>(),
cellwise_errors_pl2,
quadrature_formula,
cellwise_errors_ul2,
cellwise_errors_pl2,
pcout << "velocity error: " << u_l2 << " pressure error: " << p_l2
<< std::endl;
}
template <int dim, int degree_p, typename Number>
void StokesProblem<dim, degree_p, Number>::run()
{
pcout << std::setprecision(10);
<< " threads each) in ";
pcout << "DEBUG mode";
else
pcout << "RELEASE mode";
pcout << "\nKokkos execution space: "
<< Kokkos::DefaultExecutionSpace::name();
pcout << '\n'
<< "dim: " << dim << '\n'
<< "Element: Q" << degree_u << "-Q" << degree_p << std::endl;
unsigned int n_refinements = 10;
for (unsigned int i = 0; i < n_refinements; ++i)
{
if (i == 0)
{
tria.refine_global(2);
}
else
{
tria.refine_global(1);
}
setup_dofs();
pcout << "\nrefinement: " << i
<< ", n_dofs: " << dof_u.n_dofs() + dof_p.n_dofs() << " = "
<< dof_u.n_dofs() << " + " << dof_p.n_dofs() << std::endl;
solve();
postprocess();
}
}
}
int main(
int argc,
char **argv)
{
using namespace Step104;
const unsigned int dim = 3;
const unsigned int degree_p = 1;
StokesProblem<dim, degree_p, double> problem;
problem.run();
}
void update_ghost_values() const
MGLevelObject< RelaxationType > smoothers
DeviceVector< Number > & block(unsigned int index)
void submit_gradient(const gradient_type &gradient, const int q_point)
void submit_value(const value_type &val_in, const int q_point)
value_type get_value(const int q_point) const
gradient_type get_gradient(const int q_point) const
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
constexpr T pow(const T base, const int iexp)
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
1.9.8