step-104.h

Go to the documentation of this file.

 = 0) const override
 *       {
 *         const double pi = numbers::PI;
 *         if constexpr (dim == 2)
 *           return std::cos(pi * p[0]) * std::cos(pi * p[1]);
 *         else
 *           return std::cos(pi * p[0]) * std::cos(pi * p[1]) * std::cos(pi * p[2]);
 *       }
 *     };
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_104-Thevelocityoperator"></a> 
 * <h3>The velocity operator</h3>
 *   
 
 * 
 * The matrix-free operator for the velocity block @f$A@f$
 * given by @f$(\nabla u,\nabla v)@f$ 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().
 * 
 
 * 
 * 
 * @code
 *     template <int dim, int fe_degree, typename Number>
 *     class VelocityOperatorQuad
 *     {
 *     public:
 *       DEAL_II_HOST_DEVICE void operator()(
 *         Portable::FEEvaluation<dim, fe_degree, fe_degree + 1, dim, Number>
 *                  *fe_eval,
 *         const int q_point) const
 *       {
 *         const auto gradient_u = fe_eval->get_gradient(q_point);
 *         fe_eval->submit_gradient(gradient_u, q_point);
 *       }
 *     };
 *   
 *   
 *   
 *     template <int dim,
 *               int degree_u,
 *               int degree_p,
 *               typename Number,
 *               int n_q_points_1d>
 *     class VelocityCellOperator
 *     {
 *     public:
 *       static const unsigned int n_q_points =
 *         ::Utilities::pow(n_q_points_1d, dim);
 *   
 *       DEAL_II_HOST_DEVICE void
 *       operator()(const typename Portable::MatrixFree<dim, Number>::Data *data,
 *                  const Portable::DeviceVector<Number>                   &src,
 *                  Portable::DeviceVector<Number> &dst) const
 *       {
 *         Portable::FEEvaluation<dim, degree_u, n_q_points_1d, dim, Number> fe_u(
 *           data, velocity_dof_handler_index);
 *   
 *         fe_u.read_dof_values(src);
 *         fe_u.evaluate(EvaluationFlags::gradients);
 *   
 *         VelocityOperatorQuad<dim, degree_u, Number> quad_operation;
 *   
 *         data->for_each_quad_point(
 *           [&](const int q_point) { quad_operation(&fe_u, q_point); });
 *   
 *         fe_u.integrate(EvaluationFlags::gradients);
 *         fe_u.distribute_local_to_global(dst);
 *       }
 *     };
 *   
 *   
 * @endcode
 * 
 * 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 @f$A@f$ is symmetric anyway.
 * 
 * @code
 *     template <int dim,
 *               int degree_u,
 *               int degree_p,
 *               typename Number = double,
 *               typename VectorType =
 *                 LinearAlgebra::distributed::Vector<Number, MemorySpace::Default>,
 *               int n_q_points_1d = degree_u + 1>
 *     class PortableMFVelocityOperator : public EnableObserverPointer
 *     {
 *     public:
 *       PortableMFVelocityOperator() = default;
 *   
 *       PortableMFVelocityOperator(
 *         std::shared_ptr<Portable::MatrixFree<dim, Number>> data_in)
 *         : data(data_in)
 *       {}
 *   
 *       void reinit(std::shared_ptr<Portable::MatrixFree<dim, Number>> data_in)
 *       {
 *         data = data_in;
 *       }
 *   
 *       void initialize_dof_vector(VectorType &vec) const
 *       {
 *         data->initialize_dof_vector(vec, velocity_dof_handler_index);
 *       }
 *   
 *       types::global_dof_index m() const
 *       {
 *         return data->get_vector_partitioner(velocity_dof_handler_index)->size();
 *       }
 *   
 *       std::shared_ptr<DiagonalMatrix<
 *         LinearAlgebra::distributed::Vector<Number, MemorySpace::Default>>>
 *       get_matrix_diagonal_inverse() const
 *       {
 *         return inverse_diagonal_entries;
 *       }
 *   
 *       double el(const types::global_dof_index row,
 *                 const types::global_dof_index col) const
 *       {
 *         (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 & /* dst */, const VectorType & /* src */) const
 *       {
 *         AssertThrow(false, ExcNotImplemented());
 *       }
 *   
 *       void compute_diagonal()
 *       {
 *         Assert(data.get() != nullptr, ExcNotInitialized());
 *   
 *         this->inverse_diagonal_entries =
 *           std::make_shared<DiagonalMatrix<VectorType>>();
 *         VectorType &inverse_diagonal =
 *           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>(
 *             *data.get(),
 *             inverse_diagonal,
 *             velocity_operator_quad,
 *             EvaluationFlags::gradients,
 *             EvaluationFlags::gradients,
 *             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) {
 *             Assert(raw_diagonal[i] > 0.,
 *                    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;
 *     };
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_104-TheSchurcomplementoperator"></a> 
 * <h3>The Schur complement operator</h3>
 * 
 
 * 
 * The preconditioner requires a Schur complement
 * approximation, which is here given by a mass matrix in the
 * pressure space @f$(p,q)@f$. 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.
 * 
 * @code
 *     template <int dim,
 *               int degree_u,
 *               int degree_p,
 *               typename Number,
 *               int n_q_points_1d>
 *     class MassOperatorQuad
 *     {
 *     public:
 *       DEAL_II_HOST_DEVICE void operator()(
 *         Portable::FEEvaluation<dim, degree_p, n_q_points_1d, 1, Number> *fe_eval,
 *         const int q_point) const
 *       {
 *         fe_eval->submit_value(fe_eval->get_value(q_point), q_point);
 *       }
 *     };
 *   
 *   
 *   
 *     template <int dim,
 *               int degree_u,
 *               int degree_p,
 *               typename Number,
 *               int n_q_points_1d>
 *     class MassCellOperator
 *     {
 *     public:
 *       static const unsigned int n_q_points =
 *         ::Utilities::pow(n_q_points_1d, dim);
 *   
 *       DEAL_II_HOST_DEVICE void
 *       operator()(const typename Portable::MatrixFree<dim, Number>::Data *data,
 *                  const Portable::DeviceVector<Number>                   &src,
 *                  Portable::DeviceVector<Number> &dst) const;
 *     };
 *   
 *   
 *   
 *     template <int dim,
 *               int degree_u,
 *               int degree_p,
 *               typename Number,
 *               int n_q_points_1d>
 *     DEAL_II_HOST_DEVICE void
 *     MassCellOperator<dim, degree_u, degree_p, Number, n_q_points_1d>::operator()(
 *       const typename Portable::MatrixFree<dim, Number>::Data *data,
 *       const Portable::DeviceVector<Number>                   &src,
 *       Portable::DeviceVector<Number>                         &dst) const
 *     {
 *       Portable::FEEvaluation<dim, degree_p, n_q_points_1d, 1, Number> fe_p(
 *         data, pressure_dof_handler_index);
 *   
 *       fe_p.read_dof_values(src);
 *       fe_p.evaluate(EvaluationFlags::values);
 *   
 *       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.integrate(EvaluationFlags::values);
 *       fe_p.distribute_local_to_global(dst);
 *     }
 *   
 *   
 *   
 *     template <int dim,
 *               int degree_u,
 *               int degree_p,
 *               typename Number = double,
 *               typename VectorType =
 *                 LinearAlgebra::distributed::Vector<Number, MemorySpace::Default>,
 *               int n_q_points_1d = degree_u + 1>
 *     class PortableMFMassOperator : public EnableObserverPointer
 *     {
 *     public:
 *       PortableMFMassOperator(
 *         const std::shared_ptr<Portable::MatrixFree<dim, Number>> &data_in)
 *         : data(data_in)
 *       {}
 *   
 *       types::global_dof_index m() const
 *       {
 *         return data->get_vector_partitioner(pressure_dof_handler_index)->size();
 *       }
 *   
 *       Number el(const types::global_dof_index row,
 *                 const types::global_dof_index col) const
 *       {
 *         (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;
 *       }
 *   
 *       void compute_diagonal()
 *       {
 *         this->inverse_diagonal_entries =
 *           std::make_shared<DiagonalMatrix<VectorType>>();
 *         VectorType &inverse_diagonal =
 *           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>(
 *             *data.get(),
 *             inverse_diagonal,
 *             quad_operation,
 *             EvaluationFlags::values,
 *             EvaluationFlags::values,
 *             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) {
 *             Assert(raw_diagonal[i] > 0.,
 *                    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;
 *     };
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_104-TheStokesoperator"></a> 
 * <h3>The Stokes operator</h3>
 * 
 
 * 
 * The following set of classes provides the whole Stokes operator
 * @f{eqnarray*}{
 * \begin{bmatrix} A & B^T \\ B & 0 \end{bmatrix}.
 * @f}
 * 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.
 * 
 * @code
 *     template <int dim,
 *               int degree_u,
 *               int degree_p,
 *               typename Number,
 *               int n_q_points_1d>
 *     class StokesCellOperator
 *     {
 *     public:
 *       static const unsigned int n_q_points =
 *         ::Utilities::pow(n_q_points_1d, dim);
 *   
 *       DEAL_II_HOST_DEVICE void
 *       operator()(const typename Portable::MatrixFree<dim, Number>::Data *data,
 *                  const Portable::DeviceBlockVector<Number>              &src,
 *                  Portable::DeviceBlockVector<Number> &dst) const
 *       {
 *         Portable::FEEvaluation<dim, degree_u, n_q_points_1d, dim, Number> fe_u(
 *           data, velocity_dof_handler_index);
 *         Portable::FEEvaluation<dim, degree_p, n_q_points_1d, 1, Number> fe_p(
 *           data, pressure_dof_handler_index);
 *   
 *         fe_u.read_dof_values(src.block(0));
 *         fe_p.read_dof_values(src.block(1));
 *         fe_u.evaluate(EvaluationFlags::gradients);
 *         fe_p.evaluate(EvaluationFlags::values);
 *   
 *         data->for_each_quad_point([&](const int &q_point) {
 *           const Tensor<2, dim, Number> gradient_u = fe_u.get_gradient(q_point);
 *           const Number                 pressure_value = fe_p.get_value(q_point);
 *   
 *           Tensor<2, dim, Number> velocity_term = gradient_u;
 *           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.integrate(EvaluationFlags::gradients);
 *         fe_p.integrate(EvaluationFlags::values);
 *         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,
 *       typename Number = double,
 *       typename VectorType =
 *         LinearAlgebra::distributed::BlockVector<Number, MemorySpace::Default>,
 *       int n_q_points_1d = degree_u + 1>
 *     class PortableMFStokesOperator
 *     {
 *     public:
 *       PortableMFStokesOperator(
 *         const std::shared_ptr<Portable::MatrixFree<dim, Number>> &data_in)
 *         : data(data_in)
 *       {}
 *   
 *       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;
 *     };
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_104-TheBToperator"></a> 
 * <h3>The BT operator</h3>
 * The last ingredient of the block preconditioner is the action
 * of the block @f$B^T@f$ given by @f$-(p,\nabla \cdot v)@f$. The operator reads
 * pressure values @f$p@f$ 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.
 * 
 * @code
 *     template <int dim,
 *               int degree_u,
 *               int degree_p,
 *               typename Number,
 *               int n_q_points_1d>
 *     class BTCellOperator
 *     {
 *     public:
 *       static const unsigned int n_q_points =
 *         ::Utilities::pow(n_q_points_1d, dim);
 *   
 *       DEAL_II_HOST_DEVICE void
 *       operator()(const typename Portable::MatrixFree<dim, Number>::Data *data,
 *                  const Portable::DeviceBlockVector<Number>              &src,
 *                  Portable::DeviceBlockVector<Number> &dst) const
 *       {
 *         Portable::FEEvaluation<dim, degree_u, n_q_points_1d, dim, Number> fe_u(
 *           data, velocity_dof_handler_index);
 *         Portable::FEEvaluation<dim, degree_p, n_q_points_1d, 1, Number> fe_p(
 *           data, pressure_dof_handler_index);
 *   
 *         fe_p.read_dof_values(src.block(1));
 *         fe_p.evaluate(EvaluationFlags::values);
 *   
 *         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.integrate(EvaluationFlags::gradients);
 *         fe_u.distribute_local_to_global(dst.block(0));
 *       }
 *     };
 *   
 *   
 *   
 *     template <
 *       int dim,
 *       int degree_u,
 *       int degree_p,
 *       typename Number = double,
 *       typename VectorType =
 *         LinearAlgebra::distributed::BlockVector<Number, MemorySpace::Default>,
 *       int n_q_points_1d = degree_u + 1>
 *     class PortableMFBTOperator
 *     {
 *     public:
 *       PortableMFBTOperator(
 *         const std::shared_ptr<Portable::MatrixFree<dim, Number>> &data_in)
 *         : data(data_in)
 *       {}
 *   
 *   
 *   
 *       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);
 *   
 * @endcode
 * 
 * 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:
 * 
 * @code
 *         data->set_constrained_values(0.0,
 *                                      dst.block(0),
 *                                      velocity_dof_handler_index);
 *       }
 *   
 *     private:
 *       std::shared_ptr<Portable::MatrixFree<dim, Number>> data;
 *     };
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_104-ThePreconditionercodeBlockSchurPreconditionercode"></a> 
 * <h3>The Preconditioner <code>BlockSchurPreconditioner</code></h3>
 *   
 
 * 
 * The following class implements the block preconditioner. The class
 * takes the types of the operators for @f$A^{-1}@f$, @f$S^{-1}@f$, and
 * @f$B^T@f$ as template arguments. This is the same preconditioner used
 * in @ref step_32 "step-32" and @ref step_56 "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.
 * 
 * @code
 *     template <class AInvOperator,
 *               class SInvOperator,
 *               class BTOperator,
 *               class VectorType>
 *     class BlockSchurPreconditioner : public EnableObserverPointer
 *     {
 *     public:
 *       /**
 *        * @brief Constructor
 *        * @param A_inverse_operator Approximation of the inverse of the velocity block.
 *        * @param S_inverse_operator Approximation of the inverse Schur complement.
 *        * @param BT_operator Operator for the B^T block of the Stokes system. Note
 *        * that this operator is exactly applied while @p A_inverse_operator and @p S_inverse_operator
 *        * are approximations for the purpose of the block preconditioner.
 *        */
 *       BlockSchurPreconditioner(const AInvOperator &A_inverse_operator,
 *                                const SInvOperator &S_inverse_operator,
 *                                const BTOperator   &BT_operator);
 *   
 *       /**
 *        * Matrix-vector product with this preconditioner object.
 *        */
 *       void vmult(VectorType &dst, const VectorType &src) const;
 *   
 *     private:
 *       mutable VectorType  tmp;
 *       const AInvOperator &A_inverse_operator;
 *       const SInvOperator &S_inverse_operator;
 *       const BTOperator   &BT_operator;
 *     };
 *   
 *   
 *   
 *     template <class AInvOperator,
 *               class SInvOperator,
 *               class BTOperator,
 *               class VectorType>
 *     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,
 *               class VectorType>
 *     void
 *     BlockSchurPreconditioner<AInvOperator, SInvOperator, BTOperator, VectorType>::
 *       vmult(VectorType &dst, const VectorType &src) const
 *     {
 * @endcode
 * 
 * Allocate the temporary vector on first use. Its content doesn't matter,
 * as it will be overwritten in the vmult() call.
 * 
 * @code
 *       if (tmp.size() == 0)
 *         tmp.reinit(src);
 *   
 * @endcode
 * 
 * First apply the Schur Complement inverse operator: dst.p = S^-1 * src.p
 * 
 * @code
 *       {
 *         S_inverse_operator.vmult(dst.block(1), src.block(1));
 *         dst.block(1) *= -1.0;
 *       }
 *   
 * @endcode
 * 
 * Apply the top right block: tmp.u = -B^T * dst.p + src.u
 * 
 * @code
 *       {
 *         BT_operator.vmult(tmp, dst);
 *         tmp.block(0).sadd(-1.0, 1.0, src.block(0));
 *       }
 *   
 * @endcode
 * 
 * Finally the velocity block:
 * 
 * @code
 *       A_inverse_operator.vmult(dst.block(0), tmp.block(0));
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_104-ThemainclasscodeStokesProblemcode"></a> 
 * <h3>The main class <code>StokesProblem</code></h3>
 *   
 
 * 
 * The remaining part of this tutorial is the StokesProblem class
 * that puts everything together.
 * 
 * @code
 *     template <int dim, int degree_p, typename Number = double>
 *     class StokesProblem
 *     {
 *     public:
 *       static constexpr unsigned int degree_u = degree_p + 1;
 *   
 *       StokesProblem();
 *   
 *       void run();
 *   
 *       using VectorType =
 *         LinearAlgebra::distributed::Vector<Number, MemorySpace::Default>;
 *       using BlockVectorType =
 *         LinearAlgebra::distributed::BlockVector<Number, MemorySpace::Default>;
 *   
 *     private:
 *       void setup_dofs();
 *   
 *       void solve();
 *   
 *       void postprocess();
 *   
 *       parallel::distributed::Triangulation<dim> tria;
 *   
 *       MappingQ<dim> mapping;
 *   
 *       FESystem<dim> fe_u;
 *       FE_Q<dim>     fe_p;
 *   
 *       DoFHandler<dim> dof_u;
 *       DoFHandler<dim> dof_p;
 *   
 *       AffineConstraints<Number> constraints_u;
 *       AffineConstraints<Number> constraints_p;
 *   
 *       std::shared_ptr<Portable::MatrixFree<dim, Number>> mf_data;
 *       BlockVectorType                                    solution;
 *       BlockVectorType                                    rhs;
 *       ConditionalOStream                                 pcout;
 *     };
 *   
 *   
 *     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)
 *       , pcout(std::cout, Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
 *     {}
 *   
 * @endcode
 * 
 * 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.
 * 
 * @code
 *     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();
 *       const IndexSet  relevant_set_u =
 *         DoFTools::extract_locally_relevant_dofs(dof_u);
 *       constraints_u.reinit(owned_set_u, relevant_set_u);
 *       DoFTools::make_hanging_node_constraints(dof_u, constraints_u);
 *       VectorTools::interpolate_boundary_values(
 *         dof_u, 0, Functions::ZeroFunction<dim, Number>(dim), constraints_u);
 *       constraints_u.close();
 *   
 *       const IndexSet &owned_set_p = dof_p.locally_owned_dofs();
 *       const IndexSet  relevant_set_p =
 *         DoFTools::extract_locally_relevant_dofs(dof_p);
 *   
 *       constraints_p.reinit(owned_set_p, relevant_set_p);
 *       DoFTools::make_hanging_node_constraints(dof_p, constraints_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>>();
 *   
 *       const QGauss<1> quad(degree_p + 2);
 *       typename Portable::MatrixFree<dim, Number>::AdditionalData additional_data;
 *       additional_data.mapping_update_flags = update_values | update_gradients;
 *       mf_data->reinit(mapping, dof_handlers, constraints, quad, additional_data);
 *   
 *       {
 * @endcode
 * 
 * create the right hand side on the host and move to device:
 * 
 
 * 
 * 
 * @code
 *         LinearAlgebra::distributed::BlockVector<Number, MemorySpace::Host>
 *           rhs_host;
 *         mf_data->initialize_dof_vector(rhs_host);
 *   
 *         VectorTools::create_right_hand_side(mapping,
 *                                             dof_u,
 *                                             QGauss<dim>(degree_u + 2),
 *                                             VelocityRightHandSide<dim, Number>(),
 *                                             rhs_host.block(0),
 *                                             constraints_u);
 *   
 *         mf_data->initialize_dof_vector(rhs);
 *         rhs.block(0).import_elements(rhs_host.block(0), VectorOperation::insert);
 *         rhs.block(1).import_elements(rhs_host.block(1), VectorOperation::insert);
 *       }
 *     }
 *   
 *   
 * @endcode
 * 
 * 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 @f$A@f$, to
 * approximate the action of @f$A^{-1}@f$.
 * We approximate the Schur Complement with a Chebyshev iteration
 * applied to the pressure mass matrix (without multigrid).
 * 
 * @code
 *     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 =
 *           1e-6 * static_cast<double>(solution.size()) / time;
 *         pcout << "Stokes operator: " << time << " s, MDoFs/s: " << mdofs_p_second
 *               << std::endl;
 *         solution = 0.0;
 *       }
 *   
 *       SolverControl solver_control(1000, 1e-8 * rhs.l2_norm());
 *   
 *       SolverGMRES<BlockVectorType> solver(
 *         solver_control,
 *         typename SolverGMRES<BlockVectorType>::AdditionalData(50, true));
 *   
 *       using LevelMatrixType =
 *         PortableMFVelocityOperator<dim, degree_u, degree_p, Number>;
 *       using SmootherPreconditionerType = DiagonalMatrix<VectorType>;
 *       using SmootherType               = PreconditionChebyshev<LevelMatrixType,
 *                                                  VectorType,
 *                                                  SmootherPreconditionerType>;
 *       using MGTransferType =
 *         MGTransferMatrixFree<dim, Number, MemorySpace::Default>;
 *   
 *       const auto coarse_grid_triangulations =
 *         MGTransferGlobalCoarseningTools::create_geometric_coarsening_sequence(
 *           tria);
 *   
 *       const unsigned int max_level = coarse_grid_triangulations.size() - 1;
 * @endcode
 * 
 * Do not go down to level 0, because this will lead to slower runtime as
 * the problem becomes very small:
 * 
 * @code
 *       const unsigned int min_level = std::min(3U, max_level - 1);
 *   
 *       MGLevelObject<DoFHandler<dim>> mg_dof_handlers(min_level, max_level);
 *       MGLevelObject<AffineConstraints<Number>> mg_constraints(min_level,
 *                                                               max_level);
 *       MGLevelObject<LevelMatrixType>           mg_matrices(min_level, max_level);
 *   
 *       MGLevelObject<Portable::MGTwoLevelTransfer<dim, VectorType>> mg_transfers(
 *         min_level, max_level);
 *   
 *       std::vector<std::shared_ptr<Portable::MatrixFree<dim, Number>>>
 *         mf_data_levels;
 *   
 * @endcode
 * 
 * Prepare the operators and data structures on all levels of the multigrid
 * hierarchy
 * 
 * @code
 *       for (unsigned int level = min_level; level <= max_level; ++level)
 *         {
 *           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(),
 *                             DoFTools::extract_locally_relevant_dofs(dof_handler));
 *   
 *           DoFTools::make_hanging_node_constraints(dof_handler, constraint);
 *           DoFTools::make_zero_boundary_constraints(dof_handler, constraint);
 *           constraint.close();
 *   
 *           typename Portable::MatrixFree<dim, Number>::AdditionalData
 *             additional_data;
 *           additional_data.mapping_update_flags =
 *             update_JxW_values | update_gradients;
 *   
 *           if (level == max_level)
 * @endcode
 * 
 * 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.
 * 
 * @code
 *             mf_data_levels.emplace_back(mf_data);
 *           else
 *             {
 *               const QGauss<1> quad(degree_p + 2);
 *               mf_data_levels.emplace_back(
 *                 std::make_shared<Portable::MatrixFree<dim, Number>>());
 *   
 *               mf_data_levels.back()->reinit(
 *                 mapping, dof_handler, constraint, quad, additional_data);
 *             }
 *   
 *           mg_matrices[level].reinit(mf_data_levels.back());
 *         }
 *   
 *       mg::Matrix<VectorType> mg_matrix(mg_matrices);
 *   
 * @endcode
 * 
 * transfer operator
 * 
 * @code
 *       for (unsigned int level = min_level; level < max_level; ++level)
 *         mg_transfers[level + 1].reinit_geometric_transfer(
 *           mg_dof_handlers[level + 1],
 *           mg_dof_handlers[level],
 *           mg_constraints[level + 1],
 *           mg_constraints[level]);
 *   
 *       MGTransferType mg_transfer(mg_transfers, [&](const auto l, auto &vec) {
 *         mg_matrices[l].initialize_dof_vector(vec);
 *       });
 *   
 * @endcode
 * 
 * smoother
 * 
 * @code
 *       MGLevelObject<typename SmootherType::AdditionalData> smoother_data(
 *         min_level, max_level);
 *   
 *       for (unsigned int level = min_level; level <= max_level; ++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]);
 *   
 *           if (level == min_level)
 *             {
 * @endcode
 * 
 * Use the Chebyshev iteration as an (approximate) solver on the
 * coarsest level. In this mode @p smoothing_range is a relative
 * target tolerance and must be strictly less than one; the number
 * of iterations is then chosen automatically by setting
 * @p degree to numbers::invalid_unsigned_int. We also use more
 * CG iterations for the eigenvalue estimate because when
 * @p min_level > 0, the coarse problem can still be reasonably
 * large and badly conditioned.
 * 
 * @code
 *               smoother_data[level].smoothing_range = 1e-3;
 *               smoother_data[level].degree = numbers::invalid_unsigned_int;
 *               smoother_data[level].eig_cg_n_iterations = 40;
 *             }
 *           else
 *             {
 * @endcode
 * 
 * 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).
 * 
 * @code
 *               smoother_data[level].smoothing_range     = 5;
 *               smoother_data[level].degree              = 4;
 *               smoother_data[level].eig_cg_n_iterations = 20;
 *             }
 *         }
 *   
 *       MGSmootherRelaxation<LevelMatrixType, SmootherType, VectorType> mg_smoother;
 *       mg_smoother.initialize(mg_matrices, smoother_data);
 *   
 * @endcode
 * 
 * Estimate and print the eigenvalue spectrum of the velocity block on each
 * level. This spectrum is later used by the Chebyshev iteration.
 * 
 * @code
 *       pcout << "GMG velocity block smoothers:" << std::endl;
 *       for (unsigned int level = min_level; level <= max_level; ++level)
 *         {
 *           VectorType vec;
 *           mg_matrices[level].initialize_dof_vector(vec);
 *           auto eigenvalue_info =
 *             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;
 *         }
 *   
 * @endcode
 * 
 * coarse-grid solver
 * 
 * @code
 *       MGCoarseGridApplySmoother<VectorType> mg_coarse;
 *       mg_coarse.initialize(mg_smoother);
 *   
 * @endcode
 * 
 * put everything together
 * 
 * @code
 *       Multigrid<VectorType> mg(mg_matrix,
 *                                mg_coarse,
 *                                mg_transfer,
 *                                mg_smoother,
 *                                mg_smoother,
 *                                min_level,
 *                                max_level);
 *   
 *   
 *       ::Timer timer_smoother;
 *       ::Timer timer_transfer;
 *       ::Timer timer_coarse;
 *       ::Timer timer_residual;
 *       {
 *         timer_smoother.reset();
 *         timer_transfer.reset();
 *         timer_coarse.reset();
 *         timer_residual.reset();
 *   
 *         auto make_timer_lambda = [&](::Timer &timer) {
 *           return [&](const bool before, const unsigned int /*level*/) {
 *             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));
 *       }
 *   
 *       using APreconditionerType = PreconditionMG<dim, VectorType, MGTransferType>;
 *       APreconditionerType preconditioner_A(dof_u, mg, mg_transfer);
 *   
 *   
 *       PortableMFMassOperator<dim, degree_u, degree_p, Number> mass_operator(
 *         mf_data);
 *       mass_operator.compute_diagonal();
 *   
 *       using SPreconditionerType = PreconditionChebyshev<
 *         PortableMFMassOperator<dim, degree_u, degree_p, Number>,
 *         VectorType>;
 *   
 *       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;
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * The postprocess() function moves the solution to host memory
 * and integrates the difference to the manufactured solution to
 * compute errors.
 * 
 * @code
 *     template <int dim, int degree_p, typename Number>
 *     void StokesProblem<dim, degree_p, Number>::postprocess()
 *     {
 *       LinearAlgebra::distributed::BlockVector<Number, MemorySpace::Host>
 *         solution_host;
 *       mf_data->initialize_dof_vector(solution_host);
 *   
 *       solution_host.block(0).import_elements(solution.block(0),
 *                                              VectorOperation::insert);
 *       solution_host.block(1).import_elements(solution.block(1),
 *                                              VectorOperation::insert);
 *   
 *       constraints_u.distribute(solution_host.block(0));
 *       constraints_p.distribute(solution_host.block(1));
 *       solution_host.update_ghost_values();
 *       const double mean_pressure = VectorTools::compute_mean_value(
 *         dof_p, QGauss<dim>(degree_p + 2), solution_host.block(1), 0);
 *       solution_host.block(1).add(-mean_pressure);
 *   
 *       const QGauss<dim> quadrature_formula(degree_u + 1);
 *   
 *       Vector<double> cellwise_errors_ul2(tria.n_active_cells());
 *       Vector<double> cellwise_errors_pl2(tria.n_active_cells());
 *   
 *       VectorTools::integrate_difference(dof_u,
 *                                         solution_host.block(0),
 *                                         VelocitySolution<dim, Number>(),
 *                                         cellwise_errors_ul2,
 *                                         quadrature_formula,
 *                                         VectorTools::L2_norm);
 *       VectorTools::integrate_difference(dof_p,
 *                                         solution_host.block(1),
 *                                         PressureSolution<dim, Number>(),
 *                                         cellwise_errors_pl2,
 *                                         quadrature_formula,
 *                                         VectorTools::L2_norm);
 *   
 *       const double u_l2 = VectorTools::compute_global_error(tria,
 *                                                             cellwise_errors_ul2,
 *                                                             VectorTools::L2_norm);
 *       const double p_l2 = VectorTools::compute_global_error(tria,
 *                                                             cellwise_errors_pl2,
 *                                                             VectorTools::L2_norm);
 *   
 *       pcout << "velocity error: " << u_l2 << " pressure error: " << p_l2
 *             << std::endl;
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * The run() function prints some statistics and then performs a familiar
 * refinement loop.
 * 
 * @code
 *     template <int dim, int degree_p, typename Number>
 *     void StokesProblem<dim, degree_p, Number>::run()
 *     {
 *       pcout << std::setprecision(10);
 *       pcout << "Running on " << Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD)
 *             << " MPI ranks (with " << MultithreadInfo::n_threads()
 *             << " threads each) in ";
 *       if constexpr (running_in_debug_mode())
 *         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)
 *             {
 *               GridGenerator::hyper_cube(tria);
 *               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();
 *         }
 *     }
 *   } // namespace Step104
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_104-Thecodemaincodefunction"></a> 
 * <h3>The <code>main()</code> function</h3>
 * 
 
 * 
 * The only interesting bits here are the template arguments that
 * specify dimension and polynomial degree to be used.
 * 
 * @code
 *   int main(int argc, char **argv)
 *   {
 *     using namespace Step104;
 *     Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv);
 *   
 *     const unsigned int                   dim      = 3;
 *     const unsigned int                   degree_p = 1;
 *     StokesProblem<dim, degree_p, double> problem;
 *     problem.run();
 *   }
 * @endcode
<a name="step_104-Results"></a><h1>Results</h1>
 
 
<a name="step_104-Validation"></a><h3>Validation</h3>
 
 
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.
 
<a name="step_104-HigherOrderThroughput"></a><h3>Higher Order Throughput</h3>
 
 
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.
 
 
 
<a name="step_104-PerformancePortability"></a><h3>Performance Portability</h3>
 
 
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 @ref step_37 "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.
 
<a name="step_104-MultipleGPUs"></a><h3>Multiple GPUs</h3>
 
 
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 @ref step_40 "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
@code
  mpirun -n 4 ./step-104
@endcode
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
<a href="https://docs.rcd.clemson.edu/palmetto/">Clemson Palmetto Cluster</a> a
job submission of
@code
  salloc --gpus h100:8 --ntasks 8 --cpus-per-task 10  --mem-per-cpu 2gb --time 4:00:00
@endcode
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.
 
<a name="step_104-Possibilitiesforextensions"></a><h3> Possibilities for extensions </h3>
 
 
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!
 *
 *
<a name="step_104-PlainProg"></a>
<h1> The plain program</h1>
@include "step-104.cc"
*/