deal.II version GIT relicensing-1941-ga9013bea74 2024-10-08 06:50:00+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
The step-64 tutorial program

This tutorial depends on step-7, step-37.

Table of contents
  1. Introduction
  2. The commented program
  1. Results
  2. The plain program
This program was contributed by Bruno Turcksin and Daniel Arndt, Oak Ridge National Laboratory.

Introduction

This example shows how to implement a matrix-free method on the GPU using Kokkos for the Helmholtz equation with variable coefficients on a hypercube. The linear system will be solved using the conjugate gradient method and is parallelized with MPI.

In the last few years, heterogeneous computing in general and GPUs in particular have gained a lot of popularity. This is because GPUs offer better computing capabilities and memory bandwidth than CPUs for a given power budget. Among the architectures available in early 2019, GPUs are about 2x-3x as power efficient than server CPUs with wide SIMD for PDE-related tasks. GPUs are also the most popular architecture for machine learning. On the other hand, GPUs are not easy to program. This program explores the deal.II capabilities to see how efficiently such a program can be implemented.

We achieve a performance portable implementation by relying on the Kokkos library, which supports different devices including CPUs (serial, using OpenMP), NVidia GPUs (when configured with CUDA), AMD GPUs (with ROCm), and more. While we have tried for the interface of the matrix-free classes for the CPU and the GPU to be as close as possible, there are a few differences. However, the amount is fairly small and the use of Kokkos is limited to a few keywords.

The test case

In this example, we consider the Helmholtz problem

\begin{eqnarray*} - \nabla \cdot \nabla u + a(\mathbf x) u &=&1,\\ u &=& 0 \quad \text{on } \partial \Omega \end{eqnarray*}

where \(a(\mathbf x)\) is a variable coefficient.

We choose as domain \(\Omega=[0,1]^3\) and \(a(\mathbf x)=\frac{10}{0.05 + 2\|\mathbf x\|^2}\). Since the coefficient is symmetric around the origin but the domain is not, we will end up with a non-symmetric solution.

If you've made it this far into the tutorial, you will know how the weak formulation of this problem looks like and how, in principle, one assembles linear systems for it. Of course, in this program we will in fact not actually form the matrix, but rather only represent its action when one multiplies with it.

Moving data to and from the device

GPUs (we will use the term "device" from now on to refer to the GPU) have their own memory that is separate from the memory accessible to the CPU (we will use the term "host" from now on). A normal calculation on the device can be divided in three separate steps:

  1. the data is moved from the host to the device,
  2. the computation is done on the device,
  3. the result is moved back from the device to the host

The data movements can either be done explicitly by the user code or done automatically using UVM (Unified Virtual Memory). In deal.II, only the first method is supported. While it means an extra burden for the user, this allows for better control of data movement.

The data movement in deal.II is done using LinearAlgebra::ReadWriteVector. These vectors can be seen as buffers on the host that are used to either store data received from the device or to send data to the device. There are two types of vectors that can be used on the device:

If no memory space is specified, the default is MemorySpace::Host. MemorySpace::Default corresponds to data living on the device.

To copy a vector to/from the device, you can use import_elements(), which may also involve MPI communication:

IndexSet locally_owned_dofs, locally_relevant_dofs;
...fill the two IndexSet objects...
// Create the ReadWriteVector using an IndexSet instead of the size
LinearAlgebra::ReadWriteVector<double> owned_rw_vector(locally_owned_dofs);
...do something with the rw_vector...
// Move the data to the device:
LinearAlgebra::distributed::Vector<double, MemorySpace::Default>
distributed_vector_dev(locally_owned_dofs, MPI_COMM_WORLD);
distributed_vector_dev.import_elements(owned_rw_vector, VectorOperations::insert);
...do something with the dev_vector...
// Create a ReadWriteVector with a different IndexSet:
LinearAlgebra::ReadWriteVector<double>
relevant_rw_vector(locally_relevant_dofs);
// Move the data to the host, possibly using MPI communication:
relevant_rw_vector.import_elements(distributed_vector_dev, VectorOperations::insert);

The relevant_rw_vector is an object that stores a subset of all elements of the vector. Typically, these are the locally relevant DoFs, which implies that they overlap between different MPI processes. Consequently, the elements stored in that vector on one machine may not coincide with the ones stored by the GPU on that machine, requiring MPI communication to import them.

In all of these cases, while importing a vector, values can either be inserted (using VectorOperation::insert) or added to prior content of the vector (using VectorOperation::add).

Matrix-vector product implementation

The code necessary to evaluate the matrix-free operator on the device is very similar to the one on the host. However, there are a few differences, the main ones being that the local_apply() function in step-37 and the loop over quadrature points both need to be encapsulated in their own functors.

The commented program

First include the necessary files from the deal.II library known from the previous tutorials.

  #include <deal.II/base/conditional_ostream.h>
  #include <deal.II/base/quadrature_lib.h>
 
  #include <deal.II/dofs/dof_tools.h>
 
  #include <deal.II/fe/fe_q.h>
 
  #include <deal.II/grid/grid_generator.h>
  #include <deal.II/grid/tria.h>
 
  #include <deal.II/lac/affine_constraints.h>
  #include <deal.II/lac/la_parallel_vector.h>
  #include <deal.II/lac/precondition.h>
  #include <deal.II/lac/solver_cg.h>
 
  #include <deal.II/numerics/data_out.h>
  #include <deal.II/numerics/vector_tools.h>
 

The following ones include the data structures for the implementation of matrix-free methods on GPU:

  #include <deal.II/matrix_free/portable_fe_evaluation.h>
  #include <deal.II/matrix_free/portable_matrix_free.h>
  #include <deal.II/matrix_free/operators.h>
 
  #include <fstream>
 
 

As usual, we enclose everything into a namespace of its own:

  namespace Step64
  {
  using namespace dealii;
 
 

Class VaryingCoefficientFunctor

Next, we define a class that implements the varying coefficients we want to use in the Helmholtz operator. Later, we want to pass an object of this type to a Portable::MatrixFree object that expects the class to have an operator() that fills the values provided in the constructor for a given cell. This operator needs to run on the device, so it needs to be marked as DEAL_II_HOST_DEVICE for the compiler.

  template <int dim, int fe_degree>
  class VaryingCoefficientFunctor
  {
  public:
  VaryingCoefficientFunctor(double *coefficient)
  : coef(coefficient)
  {}
 
  operator()(const typename Portable::MatrixFree<dim, double>::Data *gpu_data,
  const unsigned int cell,
  const unsigned int q) const;
 
#define DEAL_II_HOST_DEVICE
Definition numbers.h:30

Since Portable::MatrixFree::Data doesn't know about the size of its arrays, we need to store the number of quadrature points and the number of degrees of freedom in this class to do necessary index conversions.

  static const unsigned int n_dofs_1d = fe_degree + 1;
  static const unsigned int n_local_dofs = Utilities::pow(n_dofs_1d, dim);
  static const unsigned int n_q_points = Utilities::pow(n_dofs_1d, dim);
 
  private:
  double *coef;
  };
 
 
 
constexpr T pow(const T base, const int iexp)
Definition utilities.h:966

The following function implements this coefficient. Recall from the introduction that we have defined it as \(a(\mathbf x)=\frac{10}{0.05 + 2\|\mathbf x\|^2}\)

  template <int dim, int fe_degree>
  VaryingCoefficientFunctor<dim, fe_degree>::operator()(
  const typename Portable::MatrixFree<dim, double>::Data *gpu_data,
  const unsigned int cell,
  const unsigned int q) const
  {
  const unsigned int pos = gpu_data->local_q_point_id(cell, n_q_points, q);
  const Point<dim> q_point = gpu_data->get_quadrature_point(cell, q);
 
  double p_square = 0.;
  for (unsigned int i = 0; i < dim; ++i)
  {
  const double coord = q_point[i];
  p_square += coord * coord;
  }
  coef[pos] = 10. / (0.05 + 2. * p_square);
  }
 
 
Definition point.h:111

Class HelmholtzOperatorQuad

The class HelmholtzOperatorQuad implements the evaluation of the Helmholtz operator at each quadrature point. It uses a similar mechanism as the MatrixFree framework introduced in step-37. In contrast to there, the actual quadrature point index is treated implicitly by converting the current thread index. As before, the functions of this class need to run on the device, so need to be marked as DEAL_II_HOST_DEVICE for the compiler.

  template <int dim, int fe_degree>
  class HelmholtzOperatorQuad
  {
  public:
  DEAL_II_HOST_DEVICE HelmholtzOperatorQuad(
  const typename Portable::MatrixFree<dim, double>::Data *gpu_data,
  double *coef,
  int cell)
  : gpu_data(gpu_data)
  , coef(coef)
  , cell(cell)
  {}
 
  DEAL_II_HOST_DEVICE void operator()(
  Portable::FEEvaluation<dim, fe_degree, fe_degree + 1, 1, double> *fe_eval,
  const int q_point) const;
 
  DEAL_II_HOST_DEVICE void set_matrix_free_data(
  {
  gpu_data = &data;
  }
 
  DEAL_II_HOST_DEVICE void set_cell(int new_cell)
  {
  cell = new_cell;
  }
 
  static const unsigned int n_q_points =
  ::Utilities::pow(fe_degree + 1, dim);
 
  static const unsigned int n_local_dofs = n_q_points;
 
  private:
  const typename Portable::MatrixFree<dim, double>::Data *gpu_data;
  double *coef;
  int cell;
  };
 
 

The Helmholtz problem we want to solve here reads in weak form as follows:

\begin{eqnarray*} (\nabla v, \nabla u)+ (v, a(\mathbf x) u) &=&(v,1) \quad \forall v. \end{eqnarray*}

If you have seen step-37, then it will be obvious that the two terms on the left-hand side correspond to the two function calls here:

  template <int dim, int fe_degree>
  DEAL_II_HOST_DEVICE void HelmholtzOperatorQuad<dim, fe_degree>::operator()(
  const int q_point) const
  {
  const unsigned int pos =
  gpu_data->local_q_point_id(cell, n_q_points, q_point);
 
  fe_eval->submit_value(coef[pos] * fe_eval->get_value(q_point), q_point);
  fe_eval->submit_gradient(fe_eval->get_gradient(q_point), q_point);
  }
 
 

Class LocalHelmholtzOperator

Finally, we need to define a class that implements the whole operator evaluation that corresponds to a matrix-vector product in matrix-based approaches.

  template <int dim, int fe_degree>
  class LocalHelmholtzOperator
  {
  public:

Again, the Portable::MatrixFree object doesn't know about the number of degrees of freedom and the number of quadrature points so we need to store these for index calculations in the call operator.

  static constexpr unsigned int n_dofs_1d = fe_degree + 1;
  static constexpr unsigned int n_local_dofs =
  Utilities::pow(fe_degree + 1, dim);
  static constexpr unsigned int n_q_points =
  Utilities::pow(fe_degree + 1, dim);
 
  LocalHelmholtzOperator(double *coefficient)
  : coef(coefficient)
  {}
 
  operator()(const unsigned int cell,
  const typename Portable::MatrixFree<dim, double>::Data *gpu_data,
  Portable::SharedData<dim, double> *shared_data,
  const double *src,
  double *dst) const;
 
  private:
  double *coef;
  };
 
 

This is the call operator that performs the Helmholtz operator evaluation on a given cell similar to the MatrixFree framework on the CPU. In particular, we need access to both values and gradients of the source vector and we write value and gradient information to the destination vector.

  template <int dim, int fe_degree>
  DEAL_II_HOST_DEVICE void LocalHelmholtzOperator<dim, fe_degree>::operator()(
  const unsigned int cell,
  const typename Portable::MatrixFree<dim, double>::Data *gpu_data,
  const double *src,
  double *dst) const
  {
  gpu_data, shared_data);
  fe_eval.read_dof_values(src);
  fe_eval.apply_for_each_quad_point(
  HelmholtzOperatorQuad<dim, fe_degree>(gpu_data, coef, cell));
  fe_eval.distribute_local_to_global(dst);
  }
 
 

Class HelmholtzOperator

The HelmholtzOperator class acts as wrapper for LocalHelmholtzOperator defining an interface that can be used with linear solvers like SolverCG. In particular, like every class that implements the interface of a linear operator, it needs to have a vmult() function that performs the action of the linear operator on a source vector.

  template <int dim, int fe_degree>
  class HelmholtzOperator : public Subscriptor
  {
  public:
  HelmholtzOperator(const DoFHandler<dim> &dof_handler,
  const AffineConstraints<double> &constraints);
 
  void
  &src) const;
 
  void initialize_dof_vector(
  const;
 
  void compute_diagonal();
 
  std::shared_ptr<DiagonalMatrix<
  get_matrix_diagonal_inverse() const;
 
 
 
  double el(const types::global_dof_index row,
  const types::global_dof_index col) const;
 
  private:
  std::shared_ptr<DiagonalMatrix<
  inverse_diagonal_entries;
  };
 
 
 

The following is the implementation of the constructor of this class. In the first part, we initialize the mf_data member variable that is going to provide us with the necessary information when evaluating the operator.

In the second half, we need to store the value of the coefficient for each quadrature point in every active, locally owned cell. We can ask the parallel triangulation for the number of active, locally owned cells but only have a DoFHandler object at hand. Since DoFHandler::get_triangulation() returns a Triangulation object, not a parallel::TriangulationBase object, we have to downcast the return value. This is safe to do here because we know that the triangulation is a parallel::distributed::Triangulation object in fact.

  template <int dim, int fe_degree>
  HelmholtzOperator<dim, fe_degree>::HelmholtzOperator(
  const DoFHandler<dim> &dof_handler,
  const AffineConstraints<double> &constraints)
  {
  const MappingQ<dim> mapping(fe_degree);
  const QGauss<1> quad(fe_degree + 1);
  mf_data.reinit(mapping, dof_handler, constraints, quad, additional_data);
 
 
  const unsigned int n_owned_cells =
  dynamic_cast<const parallel::TriangulationBase<dim> *>(
  &dof_handler.get_triangulation())
  ->n_locally_owned_active_cells();
  coef.reinit(Utilities::pow(fe_degree + 1, dim) * n_owned_cells);
 
  const VaryingCoefficientFunctor<dim, fe_degree> functor(coef.get_values());
  mf_data.evaluate_coefficients(functor);
  }
 
 
@ update_values
Shape function values.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.

The key step then is to use all of the previous classes to loop over all cells to perform the matrix-vector product. We implement this in the next function.

When applying the Helmholtz operator, we have to be careful to handle boundary conditions correctly. Since the local operator doesn't know about constraints, we have to copy the correct values from the source to the destination vector afterwards.

  template <int dim, int fe_degree>
  void HelmholtzOperator<dim, fe_degree>::vmult(
  const
  {
  dst = 0.;
  LocalHelmholtzOperator<dim, fe_degree> helmholtz_operator(
  coef.get_values());
  mf_data.cell_loop(helmholtz_operator, src, dst);
  mf_data.copy_constrained_values(src, dst);
  }
 
 
 
  template <int dim, int fe_degree>
  void HelmholtzOperator<dim, fe_degree>::initialize_dof_vector(
  {
  mf_data.initialize_dof_vector(vec);
  }
 
 
 
  template <int dim, int fe_degree>
  void HelmholtzOperator<dim, fe_degree>::compute_diagonal()
  {
  this->inverse_diagonal_entries.reset(
  &inverse_diagonal = inverse_diagonal_entries->get_vector();
  initialize_dof_vector(inverse_diagonal);
 
  HelmholtzOperatorQuad<dim, fe_degree> helmholtz_operator_quad(
  nullptr, coef.get_values(), -1);
 
  MatrixFreeTools::compute_diagonal<dim, fe_degree, fe_degree + 1, 0, double>(
  mf_data,
  inverse_diagonal,
  helmholtz_operator_quad,
 
  double *raw_diagonal = inverse_diagonal.get_values();
 
  Kokkos::parallel_for(
  inverse_diagonal.locally_owned_size(), KOKKOS_LAMBDA(int i) {
  Assert(raw_diagonal[i] > 0.,
  ExcMessage("No diagonal entry in a positive definite operator "
  "should be zero"));
  raw_diagonal[i] = 1. / raw_diagonal[i];
  });
  }
 
 
 
  template <int dim, int fe_degree>
  std::shared_ptr<DiagonalMatrix<
  HelmholtzOperator<dim, fe_degree>::get_matrix_diagonal_inverse() const
  {
  return inverse_diagonal_entries;
  }
 
 
 
  template <int dim, int fe_degree>
  types::global_dof_index HelmholtzOperator<dim, fe_degree>::m() const
  {
  return mf_data.get_vector_partitioner()->size();
  }
 
 
 
  template <int dim, int fe_degree>
  types::global_dof_index HelmholtzOperator<dim, fe_degree>::n() const
  {
  return mf_data.get_vector_partitioner()->size();
  }
 
 
 
  template <int dim, int fe_degree>
  double
  HelmholtzOperator<dim, fe_degree>::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);
  }
 
 
 
#define Assert(cond, exc)

Class HelmholtzProblem

This is the main class of this program. It defines the usual framework we use for tutorial programs. The only point worth commenting on is the solve() function and the choice of vector types.

  template <int dim, int fe_degree>
  class HelmholtzProblem
  {
  public:
  HelmholtzProblem();
 
  void run();
 
  private:
  void setup_system();
 
  void assemble_rhs();
 
  void solve();
 
  void output_results(const unsigned int cycle) const;
 
  MPI_Comm mpi_communicator;
 
 
  const FE_Q<dim> fe;
  DoFHandler<dim> dof_handler;
 
  IndexSet locally_owned_dofs;
  IndexSet locally_relevant_dofs;
 
  std::unique_ptr<HelmholtzOperator<dim, fe_degree>> system_matrix_dev;
 
Definition fe_q.h:554
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation

Since all the operations in the solve() function are executed on the graphics card, it is necessary for the vectors used to store their values on the GPU as well. LinearAlgebra::distributed::Vector can be told which memory space to use. It might be worth noticing that the communication between different MPI processes can be improved if the MPI implementation is GPU-aware and the configure flag DEAL_II_MPI_WITH_DEVICE_SUPPORT is enabled. (The value of this flag needs to be set at the time you call cmake when installing deal.II.)

In addition, we also keep a solution vector with CPU storage such that we can view and display the solution as usual.

The implementation of all the remaining functions of this class apart from Helmholtzproblem::solve() doesn't contain anything new and we won't further comment much on the overall approach.

  template <int dim, int fe_degree>
  HelmholtzProblem<dim, fe_degree>::HelmholtzProblem()
  : mpi_communicator(MPI_COMM_WORLD)
  , triangulation(mpi_communicator)
  , fe(fe_degree)
  , dof_handler(triangulation)
  , pcout(std::cout, Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
  {}
 
 
 
  template <int dim, int fe_degree>
  void HelmholtzProblem<dim, fe_degree>::setup_system()
  {
  dof_handler.distribute_dofs(fe);
 
  locally_owned_dofs = dof_handler.locally_owned_dofs();
  locally_relevant_dofs =
  system_rhs_dev.reinit(locally_owned_dofs, mpi_communicator);
 
  constraints.clear();
  constraints.reinit(locally_owned_dofs, locally_relevant_dofs);
  DoFTools::make_hanging_node_constraints(dof_handler, constraints);
  0,
  constraints);
  constraints.close();
 
  system_matrix_dev.reset(
  new HelmholtzOperator<dim, fe_degree>(dof_handler, constraints));
 
  ghost_solution_host.reinit(locally_owned_dofs,
  locally_relevant_dofs,
  mpi_communicator);
  system_matrix_dev->initialize_dof_vector(solution_dev);
  system_rhs_dev.reinit(solution_dev);
  }
 
 
 
void make_hanging_node_constraints(const DoFHandler< dim, spacedim > &dof_handler, AffineConstraints< number > &constraints)
IndexSet extract_locally_relevant_dofs(const DoFHandler< dim, spacedim > &dof_handler)
void interpolate_boundary_values(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const std::map< types::boundary_id, const Function< spacedim, number > * > &function_map, std::map< types::global_dof_index, number > &boundary_values, const ComponentMask &component_mask={})
STL namespace.

Unlike programs such as step-4 or step-6, we will not have to assemble the whole linear system but only the right hand side vector. This looks in essence like we did in step-4, for example, but we have to pay attention to using the right constraints object when copying local contributions into the global vector. In particular, we need to make sure the entries that correspond to boundary nodes are properly zeroed out. This is necessary for CG to converge. (Another solution would be to modify the vmult() function above in such a way that we pretend the source vector has zero entries by just not taking them into account in matrix-vector products. But the approach used here is simpler.)

At the end of the function, we can't directly copy the values from the host to the device but need to use an intermediate object of type LinearAlgebra::ReadWriteVector to construct the correct communication pattern necessary.

  template <int dim, int fe_degree>
  void HelmholtzProblem<dim, fe_degree>::assemble_rhs()
  {
  system_rhs_host(locally_owned_dofs,
  locally_relevant_dofs,
  mpi_communicator);
  const QGauss<dim> quadrature_formula(fe_degree + 1);
 
  FEValues<dim> fe_values(fe,
  quadrature_formula,
 
  const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
  const unsigned int n_q_points = quadrature_formula.size();
 
  Vector<double> cell_rhs(dofs_per_cell);
 
  std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
 
  for (const auto &cell : dof_handler.active_cell_iterators())
  if (cell->is_locally_owned())
  {
  cell_rhs = 0;
 
  fe_values.reinit(cell);
 
  for (unsigned int q_index = 0; q_index < n_q_points; ++q_index)
  {
  for (unsigned int i = 0; i < dofs_per_cell; ++i)
  cell_rhs(i) += (fe_values.shape_value(i, q_index) * 1.0 *
  fe_values.JxW(q_index));
  }
 
  cell->get_dof_indices(local_dof_indices);
  constraints.distribute_local_to_global(cell_rhs,
  local_dof_indices,
  system_rhs_host);
  }
  system_rhs_host.compress(VectorOperation::add);
 
  LinearAlgebra::ReadWriteVector<double> rw_vector(locally_owned_dofs);
  rw_vector.import_elements(system_rhs_host, VectorOperation::insert);
  system_rhs_dev.import_elements(rw_vector, VectorOperation::insert);
  }
 
 
 

This solve() function finally contains the calls to the new classes previously discussed. Here we don't use any preconditioner, i.e., precondition by the identity matrix, to focus just on the peculiarities of the Portable::MatrixFree framework. Of course, in a real application the choice of a suitable preconditioner is crucial but we have at least the same restrictions as in step-37 since matrix entries are computed on the fly and not stored.

After solving the linear system in the first part of the function, we copy the solution from the device to the host to be able to view its values and display it in output_results(). This transfer works the same as at the end of the previous function.

  template <int dim, int fe_degree>
  void HelmholtzProblem<dim, fe_degree>::solve()
  {
  system_matrix_dev->compute_diagonal();
 
  using PreconditionerType = PreconditionChebyshev<
  HelmholtzOperator<dim, fe_degree>,
  typename PreconditionerType::AdditionalData additional_data;
  additional_data.smoothing_range = 15.;
  additional_data.degree = 5;
  additional_data.eig_cg_n_iterations = 10;
  additional_data.preconditioner =
  system_matrix_dev->get_matrix_diagonal_inverse();
 
  PreconditionerType preconditioner;
  preconditioner.initialize(*system_matrix_dev, additional_data);
 
  SolverControl solver_control(system_rhs_dev.size(),
  1e-12 * system_rhs_dev.l2_norm());
  cg(solver_control);
  cg.solve(*system_matrix_dev, solution_dev, system_rhs_dev, preconditioner);
 
  pcout << " Solved in " << solver_control.last_step() << " iterations."
  << std::endl;
 
  LinearAlgebra::ReadWriteVector<double> rw_vector(locally_owned_dofs);
  rw_vector.import_elements(solution_dev, VectorOperation::insert);
  ghost_solution_host.import_elements(rw_vector, VectorOperation::insert);
 
  constraints.distribute(ghost_solution_host);
 
  ghost_solution_host.update_ghost_values();
  }
 

The output results function is as usual since we have already copied the values back from the GPU to the CPU.

While we're already doing something with the function, we might as well compute the \(L_2\) norm of the solution. We do this by calling VectorTools::integrate_difference(). That function is meant to compute the error by evaluating the difference between the numerical solution (given by a vector of values for the degrees of freedom) and an object representing the exact solution. But we can easily compute the \(L_2\) norm of the solution by passing in a zero function instead. That is, instead of evaluating the error \(\|u_h-u\|_{L_2(\Omega)}\), we are just evaluating \(\|u_h-0\|_{L_2(\Omega)}=\|u_h\|_{L_2(\Omega)}\) instead.

  template <int dim, int fe_degree>
  void HelmholtzProblem<dim, fe_degree>::output_results(
  const unsigned int cycle) const
  {
  DataOut<dim> data_out;
 
  data_out.attach_dof_handler(dof_handler);
  data_out.add_data_vector(ghost_solution_host, "solution");
  data_out.build_patches();
 
  data_out.set_flags(flags);
  data_out.write_vtu_with_pvtu_record(
  "./", "solution", cycle, mpi_communicator, 2);
 
  Vector<float> cellwise_norm(triangulation.n_active_cells());
  ghost_solution_host,
  cellwise_norm,
  QGauss<dim>(fe.degree + 2),
  const double global_norm =
  cellwise_norm,
  pcout << " solution norm: " << global_norm << std::endl;
  }
 
 
void attach_dof_handler(const DoFHandler< dim, spacedim > &)
double compute_global_error(const Triangulation< dim, spacedim > &tria, const InVector &cellwise_error, const NormType &norm, const double exponent=2.)
void integrate_difference(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const ReadVector< Number > &fe_function, const Function< spacedim, Number > &exact_solution, OutVector &difference, const Quadrature< dim > &q, const NormType &norm, const Function< spacedim, double > *weight=nullptr, const double exponent=2.)
DataOutBase::CompressionLevel compression_level

There is nothing surprising in the run() function either. We simply compute the solution on a series of (globally) refined meshes.

  template <int dim, int fe_degree>
  void HelmholtzProblem<dim, fe_degree>::run()
  {
  for (unsigned int cycle = 0; cycle < 7 - dim; ++cycle)
  {
  pcout << "Cycle " << cycle << std::endl;
 
  if (cycle == 0)
  triangulation.refine_global(1);
 
  setup_system();
 
  pcout << " Number of active cells: "
  << triangulation.n_global_active_cells() << std::endl
  << " Number of degrees of freedom: " << dof_handler.n_dofs()
  << std::endl;
 
  assemble_rhs();
  solve();
  output_results(cycle);
  pcout << std::endl;
  }
  }
  } // namespace Step64
 
 
void hyper_cube(Triangulation< dim, spacedim > &tria, const double left=0., const double right=1., const bool colorize=false)

The main() function

Finally for the main() function. By default, all the MPI ranks will try to access the device with number 0, which we assume to be the GPU device associated with the CPU on which a particular MPI rank runs. This works, but if we are running with MPI support it may be that multiple MPI processes are running on the same machine (for example, one per CPU core) and then they would all want to access the same GPU on that machine. If there is only one GPU in the machine, there is nothing we can do about it: All MPI ranks on that machine need to share it. But if there are more than one GPU, then it is better to address different graphic cards for different processes. The choice below is based on the MPI process id by assigning GPUs round robin to GPU ranks. (To work correctly, this scheme assumes that the MPI ranks on one machine are consecutive. If that were not the case, then the rank-GPU association may just not be optimal.) To make this work, MPI needs to be initialized before using this function.

  int main(int argc, char *argv[])
  {
  try
  {
  using namespace Step64;
 
  Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, 1);
 
  HelmholtzProblem<3, 3> helmholtz_problem;
  helmholtz_problem.run();
  }
  catch (std::exception &exc)
  {
  std::cerr << std::endl
  << std::endl
  << "----------------------------------------------------"
  << std::endl;
  std::cerr << "Exception on processing: " << std::endl
  << exc.what() << std::endl
  << "Aborting!" << std::endl
  << "----------------------------------------------------"
  << std::endl;
  return 1;
  }
  catch (...)
  {
  std::cerr << std::endl
  << std::endl
  << "----------------------------------------------------"
  << std::endl;
  std::cerr << "Unknown exception!" << std::endl
  << "Aborting!" << std::endl
  << "----------------------------------------------------"
  << std::endl;
  return 1;
  }
 
  return 0;
  }

Results

Since the main purpose of this tutorial is to demonstrate how to use the Portable::MatrixFree interface, not to compute anything useful in itself, we just show the expected output here:

Cycle 0
Number of active cells: 8
Number of degrees of freedom: 343
Solved in 10 iterations.
solution norm: 0.0205439
Cycle 1
Number of active cells: 64
Number of degrees of freedom: 2197
Solved in 14 iterations.
solution norm: 0.0205269
Cycle 2
Number of active cells: 512
Number of degrees of freedom: 15625
Solved in 29 iterations.
solution norm: 0.0205261
Cycle 3
Number of active cells: 4096
Number of degrees of freedom: 117649
Solved in 58 iterations.
solution norm: 0.0205261

One can make two observations here: First, the norm of the numerical solution converges, presumably to the norm of the exact (but unknown) solution. And second, the number of iterations keep increasing with each refinement of the mesh since we are only using a preconditioner based on the diagonal of the operator.

Possibilities for extensions

One could extend the tutorial to use multigrid with Chebyshev smoothers similar to step-37.

The plain program

/* ------------------------------------------------------------------------
*
* SPDX-License-Identifier: LGPL-2.1-or-later
* Copyright (C) 2019 - 2024 by the deal.II authors
*
* This file is part of the deal.II library.
*
* Part of the source code is dual licensed under Apache-2.0 WITH
* LLVM-exception OR LGPL-2.1-or-later. Detailed license information
* governing the source code and code contributions can be found in
* LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
*
* ------------------------------------------------------------------------
*
* Authors: Bruno Turcksin, Daniel Arndt, Oak Ridge National Laboratory, 2019
*/
#include <fstream>
namespace Step64
{
using namespace dealii;
template <int dim, int fe_degree>
class VaryingCoefficientFunctor
{
public:
VaryingCoefficientFunctor(double *coefficient)
: coef(coefficient)
{}
operator()(const typename Portable::MatrixFree<dim, double>::Data *gpu_data,
const unsigned int cell,
const unsigned int q) const;
static const unsigned int n_dofs_1d = fe_degree + 1;
static const unsigned int n_local_dofs = Utilities::pow(n_dofs_1d, dim);
static const unsigned int n_q_points = Utilities::pow(n_dofs_1d, dim);
private:
double *coef;
};
template <int dim, int fe_degree>
VaryingCoefficientFunctor<dim, fe_degree>::operator()(
const typename Portable::MatrixFree<dim, double>::Data *gpu_data,
const unsigned int cell,
const unsigned int q) const
{
const unsigned int pos = gpu_data->local_q_point_id(cell, n_q_points, q);
const Point<dim> q_point = gpu_data->get_quadrature_point(cell, q);
double p_square = 0.;
for (unsigned int i = 0; i < dim; ++i)
{
const double coord = q_point[i];
p_square += coord * coord;
}
coef[pos] = 10. / (0.05 + 2. * p_square);
}
template <int dim, int fe_degree>
class HelmholtzOperatorQuad
{
public:
DEAL_II_HOST_DEVICE HelmholtzOperatorQuad(
const typename Portable::MatrixFree<dim, double>::Data *gpu_data,
double *coef,
int cell)
: gpu_data(gpu_data)
, coef(coef)
, cell(cell)
{}
DEAL_II_HOST_DEVICE void operator()(
const int q_point) const;
DEAL_II_HOST_DEVICE void set_matrix_free_data(
{
gpu_data = &data;
}
DEAL_II_HOST_DEVICE void set_cell(int new_cell)
{
cell = new_cell;
}
static const unsigned int n_q_points =
::Utilities::pow(fe_degree + 1, dim);
static const unsigned int n_local_dofs = n_q_points;
private:
const typename Portable::MatrixFree<dim, double>::Data *gpu_data;
double *coef;
int cell;
};
template <int dim, int fe_degree>
DEAL_II_HOST_DEVICE void HelmholtzOperatorQuad<dim, fe_degree>::operator()(
const int q_point) const
{
const unsigned int pos =
gpu_data->local_q_point_id(cell, n_q_points, q_point);
fe_eval->submit_value(coef[pos] * fe_eval->get_value(q_point), q_point);
fe_eval->submit_gradient(fe_eval->get_gradient(q_point), q_point);
}
template <int dim, int fe_degree>
class LocalHelmholtzOperator
{
public:
static constexpr unsigned int n_dofs_1d = fe_degree + 1;
static constexpr unsigned int n_local_dofs =
Utilities::pow(fe_degree + 1, dim);
static constexpr unsigned int n_q_points =
Utilities::pow(fe_degree + 1, dim);
LocalHelmholtzOperator(double *coefficient)
: coef(coefficient)
{}
operator()(const unsigned int cell,
const typename Portable::MatrixFree<dim, double>::Data *gpu_data,
const double *src,
double *dst) const;
private:
double *coef;
};
template <int dim, int fe_degree>
DEAL_II_HOST_DEVICE void LocalHelmholtzOperator<dim, fe_degree>::operator()(
const unsigned int cell,
const typename Portable::MatrixFree<dim, double>::Data *gpu_data,
const double *src,
double *dst) const
{
gpu_data, shared_data);
fe_eval.read_dof_values(src);
HelmholtzOperatorQuad<dim, fe_degree>(gpu_data, coef, cell));
}
template <int dim, int fe_degree>
class HelmholtzOperator : public Subscriptor
{
public:
HelmholtzOperator(const DoFHandler<dim> &dof_handler,
const AffineConstraints<double> &constraints);
void
&src) const;
void initialize_dof_vector(
const;
std::shared_ptr<DiagonalMatrix<
get_matrix_diagonal_inverse() const;
double el(const types::global_dof_index row,
const types::global_dof_index col) const;
private:
std::shared_ptr<DiagonalMatrix<
inverse_diagonal_entries;
};
template <int dim, int fe_degree>
HelmholtzOperator<dim, fe_degree>::HelmholtzOperator(
const DoFHandler<dim> &dof_handler,
const AffineConstraints<double> &constraints)
{
const MappingQ<dim> mapping(fe_degree);
const QGauss<1> quad(fe_degree + 1);
mf_data.reinit(mapping, dof_handler, constraints, quad, additional_data);
const unsigned int n_owned_cells =
dynamic_cast<const parallel::TriangulationBase<dim> *>(
&dof_handler.get_triangulation())
->n_locally_owned_active_cells();
coef.reinit(Utilities::pow(fe_degree + 1, dim) * n_owned_cells);
const VaryingCoefficientFunctor<dim, fe_degree> functor(coef.get_values());
mf_data.evaluate_coefficients(functor);
}
template <int dim, int fe_degree>
void HelmholtzOperator<dim, fe_degree>::vmult(
const
{
dst = 0.;
LocalHelmholtzOperator<dim, fe_degree> helmholtz_operator(
coef.get_values());
mf_data.cell_loop(helmholtz_operator, src, dst);
mf_data.copy_constrained_values(src, dst);
}
template <int dim, int fe_degree>
void HelmholtzOperator<dim, fe_degree>::initialize_dof_vector(
{
mf_data.initialize_dof_vector(vec);
}
template <int dim, int fe_degree>
void HelmholtzOperator<dim, fe_degree>::compute_diagonal()
{
this->inverse_diagonal_entries.reset(
&inverse_diagonal = inverse_diagonal_entries->get_vector();
initialize_dof_vector(inverse_diagonal);
HelmholtzOperatorQuad<dim, fe_degree> helmholtz_operator_quad(
nullptr, coef.get_values(), -1);
MatrixFreeTools::compute_diagonal<dim, fe_degree, fe_degree + 1, 0, double>(
mf_data,
inverse_diagonal,
helmholtz_operator_quad,
double *raw_diagonal = inverse_diagonal.get_values();
Kokkos::parallel_for(
inverse_diagonal.locally_owned_size(), KOKKOS_LAMBDA(int i) {
Assert(raw_diagonal[i] > 0.,
ExcMessage("No diagonal entry in a positive definite operator "
"should be zero"));
raw_diagonal[i] = 1. / raw_diagonal[i];
});
}
template <int dim, int fe_degree>
std::shared_ptr<DiagonalMatrix<
HelmholtzOperator<dim, fe_degree>::get_matrix_diagonal_inverse() const
{
return inverse_diagonal_entries;
}
template <int dim, int fe_degree>
types::global_dof_index HelmholtzOperator<dim, fe_degree>::m() const
{
return mf_data.get_vector_partitioner()->size();
}
template <int dim, int fe_degree>
types::global_dof_index HelmholtzOperator<dim, fe_degree>::n() const
{
return mf_data.get_vector_partitioner()->size();
}
template <int dim, int fe_degree>
double
HelmholtzOperator<dim, fe_degree>::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);
}
template <int dim, int fe_degree>
class HelmholtzProblem
{
public:
HelmholtzProblem();
void run();
private:
void setup_system();
void assemble_rhs();
void solve();
void output_results(const unsigned int cycle) const;
MPI_Comm mpi_communicator;
const FE_Q<dim> fe;
DoFHandler<dim> dof_handler;
IndexSet locally_owned_dofs;
IndexSet locally_relevant_dofs;
std::unique_ptr<HelmholtzOperator<dim, fe_degree>> system_matrix_dev;
ghost_solution_host;
solution_dev;
system_rhs_dev;
};
template <int dim, int fe_degree>
HelmholtzProblem<dim, fe_degree>::HelmholtzProblem()
: mpi_communicator(MPI_COMM_WORLD)
, triangulation(mpi_communicator)
, fe(fe_degree)
, dof_handler(triangulation)
, pcout(std::cout, Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
{}
template <int dim, int fe_degree>
void HelmholtzProblem<dim, fe_degree>::setup_system()
{
dof_handler.distribute_dofs(fe);
locally_owned_dofs = dof_handler.locally_owned_dofs();
locally_relevant_dofs =
system_rhs_dev.reinit(locally_owned_dofs, mpi_communicator);
constraints.clear();
constraints.reinit(locally_owned_dofs, locally_relevant_dofs);
DoFTools::make_hanging_node_constraints(dof_handler, constraints);
0,
constraints);
constraints.close();
system_matrix_dev.reset(
new HelmholtzOperator<dim, fe_degree>(dof_handler, constraints));
ghost_solution_host.reinit(locally_owned_dofs,
locally_relevant_dofs,
mpi_communicator);
system_matrix_dev->initialize_dof_vector(solution_dev);
system_rhs_dev.reinit(solution_dev);
}
template <int dim, int fe_degree>
void HelmholtzProblem<dim, fe_degree>::assemble_rhs()
{
system_rhs_host(locally_owned_dofs,
locally_relevant_dofs,
mpi_communicator);
const QGauss<dim> quadrature_formula(fe_degree + 1);
FEValues<dim> fe_values(fe,
quadrature_formula,
const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
const unsigned int n_q_points = quadrature_formula.size();
Vector<double> cell_rhs(dofs_per_cell);
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
for (const auto &cell : dof_handler.active_cell_iterators())
if (cell->is_locally_owned())
{
cell_rhs = 0;
fe_values.reinit(cell);
for (unsigned int q_index = 0; q_index < n_q_points; ++q_index)
{
for (unsigned int i = 0; i < dofs_per_cell; ++i)
cell_rhs(i) += (fe_values.shape_value(i, q_index) * 1.0 *
fe_values.JxW(q_index));
}
cell->get_dof_indices(local_dof_indices);
constraints.distribute_local_to_global(cell_rhs,
local_dof_indices,
system_rhs_host);
}
system_rhs_host.compress(VectorOperation::add);
LinearAlgebra::ReadWriteVector<double> rw_vector(locally_owned_dofs);
rw_vector.import_elements(system_rhs_host, VectorOperation::insert);
system_rhs_dev.import_elements(rw_vector, VectorOperation::insert);
}
template <int dim, int fe_degree>
void HelmholtzProblem<dim, fe_degree>::solve()
{
system_matrix_dev->compute_diagonal();
using PreconditionerType = PreconditionChebyshev<
HelmholtzOperator<dim, fe_degree>,
typename PreconditionerType::AdditionalData additional_data;
additional_data.smoothing_range = 15.;
additional_data.degree = 5;
additional_data.eig_cg_n_iterations = 10;
additional_data.preconditioner =
system_matrix_dev->get_matrix_diagonal_inverse();
PreconditionerType preconditioner;
preconditioner.initialize(*system_matrix_dev, additional_data);
SolverControl solver_control(system_rhs_dev.size(),
1e-12 * system_rhs_dev.l2_norm());
cg(solver_control);
cg.solve(*system_matrix_dev, solution_dev, system_rhs_dev, preconditioner);
pcout << " Solved in " << solver_control.last_step() << " iterations."
<< std::endl;
LinearAlgebra::ReadWriteVector<double> rw_vector(locally_owned_dofs);
rw_vector.import_elements(solution_dev, VectorOperation::insert);
ghost_solution_host.import_elements(rw_vector, VectorOperation::insert);
constraints.distribute(ghost_solution_host);
ghost_solution_host.update_ghost_values();
}
template <int dim, int fe_degree>
void HelmholtzProblem<dim, fe_degree>::output_results(
const unsigned int cycle) const
{
DataOut<dim> data_out;
data_out.attach_dof_handler(dof_handler);
data_out.add_data_vector(ghost_solution_host, "solution");
data_out.build_patches();
data_out.set_flags(flags);
"./", "solution", cycle, mpi_communicator, 2);
Vector<float> cellwise_norm(triangulation.n_active_cells());
ghost_solution_host,
cellwise_norm,
QGauss<dim>(fe.degree + 2),
const double global_norm =
cellwise_norm,
pcout << " solution norm: " << global_norm << std::endl;
}
template <int dim, int fe_degree>
void HelmholtzProblem<dim, fe_degree>::run()
{
for (unsigned int cycle = 0; cycle < 7 - dim; ++cycle)
{
pcout << "Cycle " << cycle << std::endl;
if (cycle == 0)
triangulation.refine_global(1);
setup_system();
pcout << " Number of active cells: "
<< triangulation.n_global_active_cells() << std::endl
<< " Number of degrees of freedom: " << dof_handler.n_dofs()
<< std::endl;
assemble_rhs();
solve();
output_results(cycle);
pcout << std::endl;
}
}
} // namespace Step64
int main(int argc, char *argv[])
{
try
{
using namespace Step64;
Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, 1);
HelmholtzProblem<3, 3> helmholtz_problem;
helmholtz_problem.run();
}
catch (std::exception &exc)
{
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Exception on processing: " << std::endl
<< exc.what() << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
catch (...)
{
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Unknown exception!" << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
return 0;
}
void distribute_local_to_global(const InVector &local_vector, const std::vector< size_type > &local_dof_indices, OutVector &global_vector) const
void distribute(VectorType &vec) const
void set_flags(const FlagType &flags)
std::string write_vtu_with_pvtu_record(const std::string &directory, const std::string &filename_without_extension, const unsigned int counter, const MPI_Comm mpi_communicator, const unsigned int n_digits_for_counter=numbers::invalid_unsigned_int, const unsigned int n_groups=0) const
void add_data_vector(const VectorType &data, const std::vector< std::string > &names, const DataVectorType type=type_automatic, const std::vector< DataComponentInterpretation::DataComponentInterpretation > &data_component_interpretation={})
virtual void build_patches(const unsigned int n_subdivisions=0)
Definition data_out.cc:1062
void distribute_dofs(const FiniteElement< dim, spacedim > &fe)
const Triangulation< dim, spacedim > & get_triangulation() const
const IndexSet & locally_owned_dofs() const
types::global_dof_index n_dofs() const
size_type locally_owned_size() const
void integrate(const EvaluationFlags::EvaluationFlags integration_flag)
value_type get_value(int q_point) const
void evaluate(const EvaluationFlags::EvaluationFlags evaluate_flag)
void distribute_local_to_global(Number *dst) const
gradient_type get_gradient(int q_point) const
void read_dof_values(const Number *src)
void submit_value(const value_type &val_in, int q_point)
void submit_gradient(const gradient_type &grad_in, int q_point)
void apply_for_each_quad_point(const Functor &func)
void compute_diagonal(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, VectorType &diagonal_global, const std::function< void(FEEvaluation< dim, fe_degree, n_q_points_1d, n_components, Number, VectorizedArrayType > &)> &cell_operation, const unsigned int dof_no=0, const unsigned int quad_no=0, const unsigned int first_selected_component=0, const unsigned int first_vector_component=0)
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
Definition mpi.cc:107
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)
unsigned int local_q_point_id(const unsigned int cell, const unsigned int n_q_points, const unsigned int q_point) const
Portable::MatrixFree< dim, Number >::point_type & get_quadrature_point(const unsigned int cell, const unsigned int q_point) const