From: Wolfgang Bangerth Date: Sat, 18 May 2019 12:08:39 +0000 (-0600) Subject: Rework the first part of step-64. X-Git-Tag: v9.1.0-rc3~7^2~5 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=0eb1aacb24aec65365feeee9cb820aaa4b403da9;p=dealii.git Rework the first part of step-64. --- diff --git a/examples/step-64/step-64.cu b/examples/step-64/step-64.cu index 5d256396fa..959d68514c 100644 --- a/examples/step-64/step-64.cu +++ b/examples/step-64/step-64.cu @@ -37,8 +37,8 @@ #include #include -// This includes the data structures for the implementation of matrix-free -// methods on GPU +// The following ones include the data structures for the +// implementation of matrix-free methods on GPU: #include #include @@ -47,17 +47,22 @@ #include + +// As usual, we enclose everything into a namespace of its own: namespace Step64 { using namespace dealii; - // 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 - // CUDAWrappers::MatrixFree object that expects the class to have - // an operator that fills the values provided in the constructor for a given - // cell. + // @sect3{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 CUDAWrappers::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 devide, so it needs to be marked as `__device__` + // for the compiler. template class VaryingCoefficientFunctor { @@ -85,6 +90,9 @@ namespace Step64 + // 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 __device__ void VaryingCoefficientFunctor::operator()( const unsigned int cell, @@ -100,16 +108,21 @@ namespace Step64 double p_square = 0.; for (unsigned int i = 0; i < dim; ++i) { - double coord = q_point[i]; + const double coord = q_point[i]; p_square += coord * coord; } coef[pos] = 10. / (0.05 + 2. * p_square); } - // 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. + // @sect3{Class HelmgholtzOperatorQuad} + + // 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. As before, the functions of this class need to run on + // the device, so need to be marked as `__device__` for the + // compiler. template class HelmholtzOperatorQuad { @@ -127,12 +140,13 @@ namespace Step64 }; - // The Helmholtz problem reads in weak form + // The Helmholtz problem we want to solve here reads in weak form as follows: // @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. + // 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 __device__ void HelmholtzOperatorQuad:: operator()(CUDAWrappers::FEEvaluation *fe_eval, @@ -143,6 +157,8 @@ namespace Step64 } + // @sect3{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. @@ -200,9 +216,14 @@ namespace Step64 } + // @sect3{Class HelmholtzOperator} - // The HelmholtzOperator class acts as wrapper for LocalHelmholtzOperator - // defining an interface that can be used with linear solvers like SolverCG. + // 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 product of + // the linear operator and a source vector. template class HelmholtzOperator { @@ -222,6 +243,19 @@ namespace Step64 + // 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 doing matrix-vector products. + // + // In the second half, we need to store the value of the coefficient + // for each quadrature point in every locally owned cells. In + // actuality, however, the code stores the coefficient on + // every cell of the local DoFHandler object, including the + // ghost and artificial cells. (See the glossary entries on + // @ref @ref GlossGhostCell "ghost cells" + // and + // @ref GlossArtificialCell "artificial cells".) template HelmholtzOperator::HelmholtzOperator( const DoFHandler & dof_handler, @@ -236,13 +270,13 @@ namespace Step64 const QGauss<1> quad(fe_degree + 1); mf_data.reinit(mapping, dof_handler, constraints, quad, additional_data); - // We need to store the value of the coefficient for each quadrature point - // in every locally owned cells. + const unsigned int n_owned_cells = - std::distance( + std::distance::active_cell_iterator>( dof_handler.begin_active(), dof_handler.end()); coef.reinit(Utilities::pow(fe_degree + 1, dim) * n_owned_cells); - VaryingCoefficientFunctor functor(coef.get_values()); + + const VaryingCoefficientFunctor functor(coef.get_values()); mf_data.evaluate_coefficients(functor); }