From 7fe3d86ae8779c5dec4c9ecf6bf116813f56f331 Mon Sep 17 00:00:00 2001 From: Daniel Arndt Date: Fri, 10 May 2019 17:13:21 -0400 Subject: [PATCH] Improve documentation --- examples/step-64/doc/intro.dox | 51 +++++++++++++++------------- examples/step-64/step-64.cu | 61 +++++++++++++++------------------- 2 files changed, 55 insertions(+), 57 deletions(-) diff --git a/examples/step-64/doc/intro.dox b/examples/step-64/doc/intro.dox index 5350f84abd..a7d69b296f 100644 --- a/examples/step-64/doc/intro.dox +++ b/examples/step-64/doc/intro.dox @@ -29,11 +29,11 @@ amount is fairly small and the use of CUDA is limited to a few keywords.

The test case

In this example, we consider the Helmholtz problem @f{eqnarray*} - \nabla \cdot -\nabla u + a(\mathbf r) u &=&1,\\ u &=& 0 \quad \text{on} \partial \Omega @f} +\nabla u + a(\mathbf x) u &=&1,\\ u &=& 0 \quad \text{on } \partial \Omega @f} 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 r\|^2}$, Since the coefficient is symmetric around the origin but +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.

Moving data to and from the device

@@ -42,28 +42,32 @@ GPUs (we will use 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 is moved fromi 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, it allows a better control of data movement and more importantly it avoids to mistakenly run important kernels on the host instead of the device. -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 from the device or to send data to -the device. There are two types of vectors that can be used on the device: -LinearAlgebra::CUDAWrappers::Vector, which is similar to the more common -Vector, and LinearAlgebra::distributed::Vector, and +- LinearAlgebra::distributed::Vector, which is a regular -LinearAlgebra::distributed::Vector where we have specified which memory -space to use. If no memory space is specified, the default is MemorySpace::Host. +LinearAlgebra::distributed::Vector where we have specified which memory +space to use. + +If no memory space is specified, the default is MemorySpace::Host. Next, we show how to move data to/from the device using -LinearAlgebra::CUDAWrappers::Vector: - +LinearAlgebra::CUDAWrappers::Vector: +@code unsigned int size = 10; LinearAlgebra::CUDAWrappers::Vector vector_dev(10); LinearAlgebra::ReadWriteVector rw_vector(10); @@ -73,10 +77,11 @@ vector_dev.import(rw_vector, VectorOperations::insert); // Do some computations on the device // Move the data to the host rw_vector.import(vector_dev, VectorOperations::insert); - -Using LinearAlgebra::distributed::Vector is similar -but the import() stage may involve MPI communication: - +@endcode +Using LinearAlgebra::distributed::Vector is similar +but the `import()` stage may involve MPI communication: +@code IndexSet locally_owned_dofs, locally_relevant_dofs; // Fill the two IndexSet... LinearAlgebra::distributed::Vector @@ -91,15 +96,15 @@ distributed_vector_dev.import(owned_rw_vector, VectorOperations::insert); LinearAlgebra::ReadWriteVector relevant_rw_vector(locally_relevant_dofs); // Move the data to the host and do an MPI communication relevnt_rw_vector(distributed_vector_dev, VectorOperations::insert); - +@endcode -import() supports two kinds of strategies: VectorOperations::insert and -VectorOperations::add. +While importing a vector, values can either by inserted (using VectorOperation::insert) +or added (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 +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. diff --git a/examples/step-64/step-64.cu b/examples/step-64/step-64.cu index 8c727e5215..04453d555d 100644 --- a/examples/step-64/step-64.cu +++ b/examples/step-64/step-64.cu @@ -53,7 +53,7 @@ namespace Step64 // Define a class that implements the varying coefficients we want to use in - // the Helmholtzoperator. + // the Helmholtz operator. // Later, we want to pass an object of this type to a // CUDAWrappers::MatrixFree object that expects the class to have // an operator that fills the values provided in the constructor for a given @@ -92,7 +92,7 @@ namespace Step64 { const unsigned int pos = CUDAWrappers::local_q_point_id( cell, gpu_data, n_dofs_1d, n_q_points); - const auto q_point = + const Tensor<1, dim> q_point = CUDAWrappers::get_quadrature_point(cell, gpu_data, n_dofs_1d); @@ -107,9 +107,9 @@ namespace Step64 } - // The class HelmholtzOperatorQuad implements the evaluation of the Helmholtz - // operator in each quadrature point. It uses a similar mechanism as the - // MatrixFree framework introduced in step-37. + // 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. template class HelmholtzOperatorQuad { @@ -127,25 +127,25 @@ namespace Step64 }; - // The Helmholtz operator reads - // \begin{align*} - // -\Delta u + c\cdot u - // \end{align*} - // and consists of two parts that are correspond to the two function calls + // The Helmholtz problem reads in weak form + // @f{eqnarray*} + // (\nabla v, \nabla u)+ (v, a(\mathbf x) u) &=&(v,1) \quad \forall v. + // @f} + // The two terms on the left-hand side correspond to the two function calls // here. template __device__ void HelmholtzOperatorQuad:: operator()(CUDAWrappers::FEEvaluation *fe_eval, const unsigned int q) const { - fe_eval->submit_value(/*coef **/ fe_eval->get_value(q), q); + fe_eval->submit_value(coef * fe_eval->get_value(q), q); fe_eval->submit_gradient(fe_eval->get_gradient(q), q); } // Finally, we need to define a class that implements the whole operator - // evaluation that corresponds to matrix-vector product in matrix-based - // approaches. It corresponds + // evaluation that corresponds to a matrix-vector product in matrix-based + // approaches. template class LocalHelmholtzOperator { @@ -174,9 +174,10 @@ namespace Step64 // This is the call operator that performs the Helmholtz operator evaluation - // on a given cell similar to the MatrixFree framework. In particular, we need - // access to both values and gradients of the source vector and we write value - // and gradient information to the destincation vector. + // 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 __device__ void LocalHelmholtzOperator::operator()( const unsigned int cell, @@ -272,7 +273,6 @@ namespace Step64 { public: HelmholtzProblem(); - ~HelmholtzProblem(); void run(); @@ -289,8 +289,8 @@ namespace Step64 parallel::distributed::Triangulation triangulation; - DoFHandler dof_handler; FE_Q fe; + DoFHandler dof_handler; IndexSet locally_owned_dofs; IndexSet locally_relevant_dofs; @@ -326,21 +326,13 @@ namespace Step64 HelmholtzProblem::HelmholtzProblem() : mpi_communicator(MPI_COMM_WORLD) , triangulation(mpi_communicator) - , dof_handler(triangulation) , fe(fe_degree) + , dof_handler(triangulation) , pcout(std::cout, Utilities::MPI::this_mpi_process(mpi_communicator) == 0) {} - template - HelmholtzProblem::~HelmholtzProblem() - { - dof_handler.clear(); - } - - - template void HelmholtzProblem::setup_system() { @@ -429,11 +421,12 @@ namespace Step64 // This solve() function finally contains the calls to the new classes - // previously dicussed. Here we don't use any preconditioner, i.e. the - // identity, to focus just on the pecuiarities of the CUDA MatrixFree - // framework. Of course, in a real application the choice of a suitable - // preconditioner is crucial but we have at least the same restructions as in - // step-37 since matrix entries are computed on the fly and not stored. + // previously dicussed. Here we don't use any preconditioner, i.e. + // precondition by the identity, to focus just on the peculiarities of the + // CUDAWrappers::MatrixFree framework. Of course, in a real application the + // choice of a suitable preconditioner is crucial but we have at least the + // same restructions as in step-37 since matrix entries are computed on the + // fly and not stored. template void HelmholtzProblem::solve() { @@ -497,8 +490,8 @@ namespace Step64 } - // Nothing surprising in the run function as well. We simply compute the - // solution on a series of (globally) refined meshes. + // There is nothing surprising in the `run()` function either. We simply + // compute the solution on a series of (globally) refined meshes. template void HelmholtzProblem::run() { -- 2.39.5