This tutorial depends on step-7, step-37.
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:
- the data is moved from the host to the device,
- the computation is done on the device,
- 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;
LinearAlgebra::ReadWriteVector<double> owned_rw_vector(locally_owned_dofs);
...do something with the rw_vector...
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...
LinearAlgebra::ReadWriteVector<double>
relevant_rw_vector(locally_relevant_dofs);
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:
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)
{}
const unsigned
int q) const;
#define DEAL_II_HOST_DEVICE
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)
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 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);
}
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:
double *coef,
int cell)
: gpu_data(gpu_data)
, coef(coef)
, cell(cell)
{}
const
int q_point) const;
{
}
{
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:
double *coef;
int cell;
};
std::vector< index_type > data
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>
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 =
static constexpr unsigned int n_q_points =
LocalHelmholtzOperator(double *coefficient)
: coef(coefficient)
{}
operator()(const unsigned
int cell,
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>
const unsigned int cell,
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>
{
public:
void
&src) const;
void initialize_dof_vector(
const;
void compute_diagonal();
get_matrix_diagonal_inverse() const;
private:
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(
{
mf_data.reinit(mapping, dof_handler, constraints, quad, additional_data);
const unsigned int n_owned_cells =
&dof_handler.get_triangulation())
->n_locally_owned_active_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.
UpdateFlags mapping_update_flags
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) {
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>
HelmholtzOperator<dim, fe_degree>::get_matrix_diagonal_inverse() const
{
return inverse_diagonal_entries;
}
template <int dim, int fe_degree>
{
return mf_data.get_vector_partitioner()->size();
}
template <int dim, int fe_degree>
{
return mf_data.get_vector_partitioner()->size();
}
template <int dim, int fe_degree>
double
{
(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;
std::unique_ptr<HelmholtzOperator<dim, fe_degree>> system_matrix_dev;
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.
ghost_solution_host;
solution_dev;
system_rhs_dev;
};
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)
, fe(fe_degree)
, 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);
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)
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);
quadrature_formula,
const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
const unsigned int n_q_points = quadrature_formula.size();
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);
}
}
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();
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);
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;
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
{
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);
ghost_solution_host,
cellwise_norm,
const double global_norm =
cellwise_norm,
pcout << " solution norm: " << global_norm << std::endl;
}
void attach_dof_handler(const DoFHandler< dim, spacedim > &)
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)
setup_system();
pcout << " Number of active cells: "
<< " Number of degrees of freedom: " << dof_handler.n_dofs()
<< std::endl;
assemble_rhs();
solve();
output_results(cycle);
pcout << std::endl;
}
}
}
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;
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
#include <fstream>
namespace Step64
{
template <int dim, int fe_degree>
class VaryingCoefficientFunctor
{
public:
VaryingCoefficientFunctor(double *coefficient)
: coef(coefficient)
{}
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 unsigned int cell,
const unsigned int q) const
{
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:
double *coef,
int cell)
: gpu_data(gpu_data)
, coef(coef)
, cell(cell)
{}
const int q_point) const;
{
}
{
cell = new_cell;
}
static const unsigned int n_q_points =
static const unsigned int n_local_dofs = n_q_points;
private:
double *coef;
int cell;
};
template <int dim, int fe_degree>
const int q_point) const
{
const unsigned int pos =
}
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 =
static constexpr unsigned int n_q_points =
LocalHelmholtzOperator(double *coefficient)
: coef(coefficient)
{}
operator()(const unsigned int cell,
const double *src,
double *dst) const;
private:
double *coef;
};
template <int dim, int fe_degree>
const unsigned int cell,
const double *src,
double *dst) const
{
gpu_data, shared_data);
HelmholtzOperatorQuad<dim, fe_degree>(gpu_data, coef, cell));
}
template <int dim, int fe_degree>
{
public:
void
&src) const;
void initialize_dof_vector(
const;
get_matrix_diagonal_inverse() const;
private:
inverse_diagonal_entries;
};
template <int dim, int fe_degree>
HelmholtzOperator<dim, fe_degree>::HelmholtzOperator(
{
mf_data.reinit(mapping, dof_handler, constraints, quad, additional_data);
const unsigned int n_owned_cells =
->n_locally_owned_active_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(
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>
HelmholtzOperator<dim, fe_degree>::get_matrix_diagonal_inverse() const
{
return inverse_diagonal_entries;
}
template <int dim, int fe_degree>
{
return mf_data.get_vector_partitioner()->size();
}
template <int dim, int fe_degree>
{
return mf_data.get_vector_partitioner()->size();
}
template <int dim, int fe_degree>
double
{
(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();
private:
void setup_system();
void assemble_rhs();
void solve();
void output_results(const unsigned int cycle) const;
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)
, fe(fe_degree)
{}
template <int dim, int fe_degree>
void HelmholtzProblem<dim, fe_degree>::setup_system()
{
locally_relevant_dofs =
system_rhs_dev.reinit(locally_owned_dofs, mpi_communicator);
constraints.
reinit(locally_owned_dofs, locally_relevant_dofs);
0,
constraints);
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);
quadrature_formula,
const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
const unsigned int n_q_points = quadrature_formula.size();
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);
local_dof_indices,
system_rhs_host);
}
}
template <int dim, int fe_degree>
void HelmholtzProblem<dim, fe_degree>::solve()
{
system_matrix_dev->compute_diagonal();
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);
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;
ghost_solution_host.update_ghost_values();
}
template <int dim, int fe_degree>
void HelmholtzProblem<dim, fe_degree>::output_results(
const unsigned int cycle) const
{
"./", "solution", cycle, mpi_communicator, 2);
ghost_solution_host,
cellwise_norm,
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)
setup_system();
pcout << " Number of active cells: "
<<
" Number of degrees of freedom: " << dof_handler.
n_dofs()
<< std::endl;
assemble_rhs();
solve();
output_results(cycle);
pcout << std::endl;
}
}
}
int main(int argc, char *argv[])
{
try
{
using namespace Step64;
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)
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
Number * get_values() 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)
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
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
1.9.8