matrix that defines the index mapping from local degrees of freedom in the
current cell to the global degrees of freedom. The information from which this
operator can be built is usually encoded in the <code>local_dof_indices</code>
-variable is used in the assembly calls filling matrices in deal.II. Moreover,
+variable and is used in the assembly calls filling matrices in deal.II. Here,
<i>A</i><sub>cell</sub> denotes the cell matrix associated with <i>A</i>.
If we are to perform a matrix-vector product, we can hence use that
needs to be called on every cell. In this process, the transformation is
applied while computing the local matrix elements.
-In a matrix-free implementation, however, we will compute those integrals much
-more often than once. For example, iterative solvers will apply the matrix
-many times, and so we can think about whether we may be able to cache
-something between different operator applications, i.e., integral
+In a matrix-free implementation, however, we will compute those integrals very
+often because iterative solvers will apply the matrix many times during the
+solution process. Therefore, we need to think about whether we may be able to
+cache some data that gets reused in the operator applications, i.e., integral
computations. On the other hand, we realize that we must not cache too much
data since otherwise we get back to the situation where memory access becomes
-the dominating factor.
+the dominating factor. Therefore, we will not store the transformed gradients
+in the matrix <i>B</i>, as they would in general be different on every element
+for curved meshes.
The trick is to factor out the Jacobian transformation and first apply the
-gradient on the reference cell only. That transforms the vector of values on
-the local dofs to a vector of gradients on the quadrature points. There, we
-first apply the Jacobian that we factored out from the gradient, apply the
-weights of the quadrature, and finally apply the transposed Jacobian for
-preparing the third loop which again uses the gradients on the unit cell.
+gradient on the reference cell only. This operation interpolates the vector of
+values on the local dofs to a vector of (unit-coordinate) gradients on the
+quadrature points. There, we first apply the Jacobian that we factored out
+from the gradient, then apply the weights of the quadrature, and finally apply
+the transposed Jacobian for preparing the third loop which tests by the
+gradients on the unit cell and sums over quadrature points.
Let us again write this in terms of matrices. Let the matrix
<i>B</i><sub>cell</sub> denote the cell-related gradient matrix, with each row
gradients and how we initialize it to the reference cell. The actual
derivative data is then applied by the inverse, transposed Jacobians (deal.II
calls the Jacobian matrix from real to unit cell inverse_jacobian, as the
-forward transformation is from unit to real cell). Note that the factor
+forward transformation is from unit to real cell). The factor
$J_\mathrm{cell}^{-1} D_\mathrm{cell} J_\mathrm{cell}^{-\mathrm T}$ is
block-diagonal over quadrature. In this form, one realizes that variable
coefficients (possibly expressed through a tensor) and general grid topologies
-with Jacobian transformations have a similar effect on the coefficient that
-transforms the unit-cell derivatives.
-
-Finally, we are using tensor product basis functions and now that we have
-separated out the gradient on the reference cell <i>B</i><sub>ref_cell</sub>,
-we can exploit the tensor-product structure to further reduce the
-complexity. We illustrate this in two space dimensions, but the same technique
-can be used in higher dimensions. On the reference cell, the basis functions
-are of the tensor product form $\phi(x,y,z) = \varphi_i(x) \varphi_j(y)$. The
-part of the matrix <i>B</i><sub>ref_cell</sub> that computes the first
-component has the form $B_\mathrm{sub\_cell}^x = B_\mathrm{grad,x} \otimes
-B_\mathrm{val,y}$, where <i>B</i><sub>grad,x</sub> and
-<i>B</i><sub>val,y</sub> contain the evaluation of all the 1D basis functions
-on all the 1D quadrature points. Forming a matrix <i>U</i> with <i>U(j,i)</i>
-containing the coefficient belonging to basis function $\varphi_i(x)
-\varphi_j(y)$, we get $(B_\mathrm{grad,x} \otimes
+with Jacobian transformations have a similar effect on the coefficient
+transforming the unit-cell derivatives.
+
+At this point, one might wonder why we store the matrix
+$J_\mathrm{cell}^{-\mathrm T}$ and the coefficient separately, rather than
+only the complete factor $J_\mathrm{cell}^{-1} D_\mathrm{cell}
+J_\mathrm{cell}^{-\mathrm T}$. The latter would use less memory because the
+tensor is symmetric with six independent values in 3D, whereas for the former
+we would need nine entries for the inverse transposed Jacobian, one for the
+quadrature weight and Jacobian determinant, and one for the coefficient,
+totaling to 11 doubles. The reason is that the former approach allows for
+implementing generic differential operators through a common framework of
+cached data, whereas the latter specifically stores the coefficient for the
+Laplacian. In case applications demand for it, this specialization could pay
+off and would be worthwhile to consider. Note that the implementation in
+deal.II is smart enough to detenct Cartesian or affine geometries where the
+Jacobian is constant throughout the cell and needs not be stored for every
+cell (and indeed often is the same over different cells as well).
+
+The final optimization that is most crucial from an operation count point of
+view is to make use of the tensor product structure in the basis
+functions. This is possible because we have factored out the gradient from the
+reference cell operation described by <i>B</i><sub>ref_cell</sub>, i.e., an
+interpolation operation over the completely regular data fields of the
+reference cell. We illustrate the process of complexity reduction in two space
+dimensions, but the same technique can be used in higher dimensions. On the
+reference cell, the basis functions are of the tensor product form
+$\phi(x,y,z) = \varphi_i(x) \varphi_j(y)$. The part of the matrix
+<i>B</i><sub>ref_cell</sub> that computes the first component has the form
+$B_\mathrm{sub\_cell}^x = B_\mathrm{grad,x} \otimes B_\mathrm{val,y}$, where
+<i>B</i><sub>grad,x</sub> and <i>B</i><sub>val,y</sub> contain the evaluation
+of all the 1D basis functions on all the 1D quadrature points. Forming a
+matrix <i>U</i> with <i>U(j,i)</i> containing the coefficient belonging to
+basis function $\varphi_i(x) \varphi_j(y)$, we get $(B_\mathrm{grad,x} \otimes
B_\mathrm{val,y})u_\mathrm{cell} = B_\mathrm{val,y} U B_\mathrm{grad,x}$. This
reduces the complexity for computing this product from $p^4$ to $2 p^3$, where
<i>p</i>-1 is the degree of the finite element (i.e., equivalently, <i>p</i>
is the number of shape functions in each coordinate direction), or $p^{2d}$ to
-$d p^{d+1}$ in general. Note that the concentration on the polynomial degree
-in terms of complexity is relevant as the polynomial degree is increased
-rather than the grid resolution. Good algorithms for moderate degrees like the
-one used here are linear in the polynomial degree independent on the
-dimension, as opposed to matrix-based schemes or naive evaluation through
-FEValues. These techniques have been established in the spectral element
-community since the 1980s.
+$d p^{d+1}$ in general. The reason why we look at the complexity in terms of
+the polynomial degree is since we want to be able to go to high degrees and
+possibly increase the polynomial degree <i>p</i> instead of the grid
+resolution. Good algorithms for moderate degrees like the ones used here are
+linear in the polynomial degree independent on the dimension, as opposed to
+matrix-based schemes or naive evaluation through FEValues. The techniques used
+in the implementations of deal.II have been established in the spectral
+element community since the 1980s.
Implementing a matrix-free and cell-based finite element operator requires a
-somewhat different design compared to the usual matrix assembly codes shown in
-previous tutorial programs. The data structures for doing this are the
-MatrixFree class that collects all data and issues a (parallel) loop over all
-cells and the FEEvaluation class that evaluates finite element basis functions
-by making use of the tensor product structure.
-
+somewhat different program design as compared to the usual matrix assembly
+codes shown in previous tutorial programs. The data structures for doing this
+are the MatrixFree class that collects all data and issues a (parallel) loop
+over all cells and the FEEvaluation class that evaluates finite element basis
+functions by making use of the tensor product structure.
The implementation of the matrix-free matrix-vector product shown in this
tutorial is slower than a matrix-vector product using a sparse matrix for
transfer during computations. The impact of reduced memory transfer is
particularly beneficial when working on a multicore processor where several
processing units share access to memory. In that case, an algorithm which is
-computation bound will show almost perfect parallel speedup, whereas an
-algorithm that is bound by memory transfer might not achieve similar speedup
-(even when the work is perfectly parallel and one could expect perfect scaling
-like in sparse matrix-vector products). An additional gain with this
-implementation is that we do not have to build the sparse matrix itself, which
-can also be quite expensive depending on the underlying differential
-equation. Moreover, the above framework is simple to generalize to nonlinear
-operations, as we demonstrate in step-48.
+computation bound will show almost perfect parallel speedup (apart from
+possible changes of the processor's clock frequency through turbo modes
+depending on how many cores are active), whereas an algorithm that is bound by
+memory transfer might not achieve similar speedup (even when the work is
+perfectly parallel and one could expect perfect scaling like in sparse
+matrix-vector products). An additional gain with this implementation is that
+we do not have to build the sparse matrix itself, which can also be quite
+expensive depending on the underlying differential equation. Moreover, the
+above framework is simple to generalize to nonlinear operations, as we
+demonstrate in step-48.
<h3>Combination with multigrid</h3>
Above, we have gone to significant lengths to implement a matrix-vector
product that does not actually store the matrix elements. In many user codes,
-however, one wants more than just performing some number of matrix-vector
-products — one wants to do as little of these operations as possible
-when solving linear equation systems. In theory, we could use the CG method
-without preconditioning; however, that would not be very efficient. Rather,
-one uses preconditioners for improving speed. Unfortunately, most of the more
-frequently used preconditioners such as SSOR, ILU or algebraic multigrid (AMG)
-cannot be used here because their implementation requires knowledge of the
-elements of the system matrix.
+however, one wants more than just doing a few matrix-vector products —
+one wants to do as few of these operations as possible when solving linear
+systems. In theory, we could use the CG method without preconditioning;
+however, that would not be very efficient for the Laplacian. Rather,
+preconditioners are used for increasing the speed of
+convergence. Unfortunately, most of the more frequently used preconditioners
+such as SSOR, ILU or algebraic multigrid (AMG) cannot be used here because
+their implementation requires knowledge of the elements of the system matrix.
One solution is to use geometric multigrid methods as shown in step-16. They
-are known to be very fast, and they are suitable for our purpose since they
-can be designed based purely on matrix-vector products related to the a
-colletion of cells. All one needs to do is to find a smoother that works with
-matrix-vector products only (our choice requires knowledge of the diagonal
-entries of the matrix, though). One such candidate would be a damped Jacobi
-iteration, but that is often not sufficiently good in damping high-frequency
-errors. A Chebyshev iteration around the Jacobi method, eventually, is our
-choice in this tutorial program. It can be seen as an extension of the Jacobi
-method by using Chebyshev polynomials. With degree zero, the Jacobi method
-with optimal damping parameter is retrieved, whereas higher order corrections
-improve the smoothing properties if some parameters are suitably chosen. The
-effectiveness of Chebyshev smoothing in multigrid has been demonstrated, e.g.,
-in the article <a
-href="http://www.sciencedirect.com/science/article/pii/S0021999103001943">
+are known to be very fast, and they are suitable for our purpose since all
+ingredients, including the transfer between different grid levels, can be
+expressed in terms of matrix-vector products related to a colletion of
+cells. All one needs to do is to find a smoother that is based on
+matrix-vector products rather than all the matrix entries. One such candidate
+would be a damped Jacobi iteration that requires access to the matrix
+diagonal, but it is often not sufficiently good in damping all high-frequency
+errors. The properties of the Jacobi method can be improved by iterating it a
+few times with the so-called Chebyshev iteration. The Chebyshev iteration is
+described by a polynomial expression of the matrix-vector product where the
+coefficients can be chosen to achieve certain properties, in this case to
+smooth the high-frequency components of the error which are associated to the
+eigenvalues of the Jacobi-preconditioned matrix. At degree zero, the Jacobi
+method with optimal damping parameter is retrieved, whereas higher order
+corrections are used to improve the smoothing properties. The effectiveness of
+Chebyshev smoothing in multigrid has been demonstrated, e.g., in the article
+<a href="http://www.sciencedirect.com/science/article/pii/S0021999103001943">
<i>M. Adams, M. Brezina, J. Hu, R. Tuminaro. Parallel multigrid smoothers:
polynomial versus Gauss–Seidel, J. Comput. Phys. 188:593–610,
-2003</i> </a>. This publication also identifies one more advantage of
+2003</i></a>. This publication also identifies one more advantage of
Chebyshev smoothers that we exploit here, namely that they are easy to
parallelize, whereas SOR/Gauss–Seidel smoothing relies on substitutions,
for which a naive parallelization works on diagonal sub-blocks of the matrix,
data from several cells. Remember that it is not related to the spatial
dimension and the number of elements e.g. in a Tensor or Point.
-Note that vectorization depends on the CPU that is used for deal.II. In order
-to generate the fastest kernels of FEEvaluation for your computer, you should
-compile deal.II with the so-called <i>native</i> processor variant. When using
-the gcc compiler, it can be enabled by setting the variable
-<tt>CMAKE_CXX_FLAGS</tt> to <tt>"-march=native"</tt> in the cmake build
-settings (on the command line, specify
-<tt>-DCMAKE_CXX_FLAGS="-march=native"</tt>, see the deal.II README for more
-information). Similar options exist for other compilers.
+Note that vectorization depends on the CPU the code is running on and for
+which the code is compiled. In order to generate the fastest kernels of
+FEEvaluation for your computer, you should compile deal.II with the so-called
+<i>native</i> processor variant. When using the gcc compiler, it can be
+enabled by setting the variable <tt>CMAKE_CXX_FLAGS</tt> to
+<tt>"-march=native"</tt> in the cmake build settings (on the command line,
+specify <tt>-DCMAKE_CXX_FLAGS="-march=native"</tt>, see the deal.II README for
+more information). Similar options exist for other compilers.
<h3>Running multigrid on large-scale parallel computers</h3>
As mentioned above, all components in the matrix-free framework can easily be
-parallelized with MPI using domain decomposition. Thanks to the good support
-of large-scale parallel mesh framework provided by p4est (see step-40 for
-details) in deal.II, and the fact that cell-based loops with matrix-free
-evaluation <i>only</i> need a decomposition of the mesh into chunks of roughly
-equal size on each processor, there is relatively little to do for extending
-the program to MPI. While other tutorial programs using MPI have relied on
-either PETSc or Trilinos, this program uses deal.II's own parallel vector
-facilities.
+parallelized with MPI using domain decomposition. Thanks to the easy access to
+large-scale parallel meshes through p4est (see step-40 for details) in
+deal.II, and the fact that cell-based loops with matrix-free evaluation
+<i>only</i> need a decomposition of the mesh into chunks of roughly equal size
+on each processor, there is relatively little to do to write a parallel
+program working with distributed memory. While other tutorial programs using
+MPI have relied on either PETSc or Trilinos, this program uses deal.II's own
+parallel vector facilities.
The deal.II parallel vector class, LinearAlgebra::distributed::Vector, holds
-the processor-local part of the solution as well as information on and data
-fields for the ghost DoFs, i.e. DoFs that are owned by a remote processor but
-needed on cells that are treated by the present processor. Moreover, it holds
-the MPI-send information for DoFs that are owned locally but needed by other
-processors. A benefit of the design of that vector class is the way ghosted
-entries are accessed. In the storage scheme of the vector, the array of data
-does extend beyond the processor-local part of the solution with further
-vector entries available for the ghost degrees of freedom by a continguous
-index range. (Note that the index range depends on the exact configuration of
-the mesh.) Moreover, it holds the MPI-send information for DoFs that are owned
-locally but needed by other processors. Since matrix-free operations can be
-thought of doing linear algebra that is performance critical, and
-performance-critical code cannot waste time on doing MPI-global to MPI-local
-index translations, the availability of such access functionality in index
-spaces local to one MPI rank is fundamental. The way things are accessed here
-is a direct array access. A user can use the query
-LinearAlgebra::distributed::Vector::local_element for that purpose, but it is
-actually rarely used because all of this happens internally in FEEvaluation.
+the processor-local part of the solution as well as data fields for ghosted
+DoFs, i.e. DoFs that are owned by a remote processor but accessed by cells
+that are owned by the present processor. In the @ref GlossLocallyActiveDof
+"glossary" these degrees of freedom are referred to as locally active degrees
+of freedom. The function MatrixFree::initialize_dof_vector() prodides a method
+that sets this design. Note that hanging nodes can relate to additional
+ghosted degrees of freedom that must be included in the distributed vector but
+are not part of the locally active DoFs in the sense of the @ref
+GlossLocallyActiveDof "glossary". Moreover, the distributed vector holds the
+some MPI metadata for DoFs that are owned locally but needed by other
+processors. A benefit of the design of this vector class is the way ghosted
+entries are accessed. In the storage scheme of the vector, the data array
+extends beyond the processor-local part of the solution with further vector
+entries available for the ghosted degrees of freedom. This gives a continguous
+index range for all locally active degrees of freedom. (Note that the index
+range depends on the exact configuration of the mesh.) Since matrix-free
+operations can be thought of doing linear algebra that is performance
+critical, and performance-critical code cannot waste time on doing MPI-global
+to MPI-local index translations, the availability of an index spaces local to
+one MPI rank is fundamental. The way things are accessed here is a direct
+array access. This is provided through
+LinearAlgebra::distributed::Vector::local_element(), but it is actually rarely
+needed because all of this happens internally in FEEvaluation.
The design of LinearAlgebra::distributed::Vector is similar to the
PETScWrappers::MPI::Vector and TrilinosWrappers::MPI::Vector data types we
@endcode
As in step-16, we see that the number of CG iterations remains constant with
-increasing number of degrees of freedom. The constant number of iterations
+increasing number of degrees of freedom. A constant number of iterations
(together with optimal computational properties) means that the computing time
approximately quadruples as the problem size quadruples from one cycle to the
next. A look at the memory consumption reveals that the code is also very
cheaper than the setup, despite not building a matrix (approximately half of
which is spent in the DoFHandler::distribute_dofs() and
DoFHandler::distributed_mg_dofs() calls). This shows the high efficiency of
-this approach, but also the fact that most benefit will be obtained when
-solving multiple systems on the same mesh and DoFHandler structure.
+this approach, but also that the deal.II data structures are quite expensive
+to set up and the setup cost must be amortized over several system solves.
Not much changes if we run the program in three spatial dimensions, with the
exception that problem sizes grow by a factor eight when refining the mesh and
Since it is so easy, we look at what happens if we increase the polynomial
degree. When selecting the degree as four in 3D, i.e., on $\mathcal Q_4$
-elements, we get the following program output:
+elements, by changing the line <code>const unsigned int
+degree_finite_element=4;</code> at the top of the program, we get the
+following program output:
@code
Cycle 0
Time solve (7 iterations) (CPU/wall) 26.3705s/27.1077s
@endcode
-Since $\mathcal Q_4$ elements on a mesh correspond to $\mathcal Q_2$ elements
-on half the mesh size, we can start to compare the run time at cycle 4 with
-fourth degree polynomials with cycle 5 at quadratic polynomials, both at 2.1
-million degrees of freedom. The surprising effect is that the solver for
+Since $\mathcal Q_4$ elements on a certain mesh correspond to $\mathcal Q_2$
+elements on half the mesh size, we can compare the run time at cycle 4 with
+fourth degree polynomials with cycle 5 using quadratic polynomials, both at
+2.1 million degrees of freedom. The surprising effect is that the solver for
$\mathcal Q_4$ element is actually slightly faster than for the quadratic
case, despite using one more linear iteration. The effect that higher-degree
polynomials are similarly fast or even faster than lower degree ones is one of
factorization, see the <a
href="http://dx.doi.org/10.1016/j.compfluid.2012.04.012">matrix-free
paper</a>. This is fundamentally different to matrix-based methods that get
-more expensive as the degree increases and the coupling gets denser.
+more expensive per unknown as the polynomial degree increases and the coupling
+gets denser.
-The improved efficiency is also reflected in the initialization routines.
+In addition, also the setup gets a bit cheaper for higher order, which is
+because fewer elements need to be set up.
Finally, let us look at the timings with degree 8, which corresponds to
another round of mesh refinement in the lower order methods:
computes a 729 x 729 matrix on each cell and throws away everything but the
diagonal. The solver times, however, are again very close to the quartic case,
showing that the linear increase with the polynomial degree that is
-theoretically expected is offset by better computational characteristics and
-the fact that higher order methods have a smaller share of degrees of freedom
-living on several cells that slow the method down.
+theoretically expected is almost completely offset by better computational
+characteristics and the fact that higher order methods have a smaller share of
+degrees of freedom living on several cells that add to the evaluation
+complexity.
<h3>Comparison with a sparse matrix</h3>
In order to understand the capabilities of the matrix-free implementation, we
-compare the performance on the 3d example above with a SparseMatrix
-implementation and we measure the computation times for both initialization of
+compare the performance of the 3d example above with a SparseMatrix
+implementation by measuring both the computation times for the initialization of
the problem (distribute DoFs, setup and assemble matrices, setup multigrid
structures) and the actual solution for the matrix-free variant and the
variant based on sparse matrices. We base the preconditioner on float
<td align="right">125</td>
<td align="center">0.0048s</td>
<td align="center">0.00075s</td>
- <td align="center">0.0025s</td>
- <td align="center">0.00033s</td>
+ <td align="center">0.0023s</td>
+ <td align="center">0.00092s</td>
</tr>
<tr>
<td align="right">729</td>
<td align="center">0.014s</td>
<td align="center">0.0022s</td>
- <td align="center">0.0026s</td>
- <td align="center">0.0018s</td>
+ <td align="center">0.0029s</td>
+ <td align="center">0.0028s</td>
</tr>
<tr>
<td align="right">4,913</td>
<td align="center">0.10s</td>
<td align="center">0.012s</td>
- <td align="center">0.017s</td>
- <td align="center">0.013s</td>
+ <td align="center">0.014s</td>
+ <td align="center">0.011s</td>
</tr>
<tr>
<td align="right">35,937</td>
<td align="center">0.80s</td>
<td align="center">0.14s</td>
- <td align="center">0.11s</td>
- <td align="center">0.071s</td>
+ <td align="center">0.087s</td>
+ <td align="center">0.065s</td>
</tr>
<tr>
<td align="right">274,625</td>
<td align="center">5.93s</td>
<td align="center">1.05s</td>
- <td align="center">0.82s</td>
- <td align="center">0.51s</td>
+ <td align="center">0.60s</td>
+ <td align="center">0.43s</td>
</tr>
<tr>
<td align="right">2,146,689</td>
<td align="center">46.7s</td>
<td align="center">8.44s</td>
- <td align="center">6.39s</td>
- <td align="center">3.77s</td>
+ <td align="center">4.96s</td>
+ <td align="center">3.56s</td>
</tr>
</table>
that the times usually go up by a factor eight, too (as the solver iterations
are constant at six). The main deviation is in the sparse matrix between 5k
and 36k degrees of freedom, where the time increases by a factor 12. This is
-the threshold when the cache in the processor can no longer hold all data
-necessary for the matrix-vector products and all matrix elements must be
+the threshold where the (L3) cache in the processor can no longer hold all
+data necessary for the matrix-vector products and all matrix elements must be
fetched from main memory.
Of course, this picture does not necessarily translate to all cases, as there
// variable coefficient, we further implement a method that can fill the
// coefficiient values.
//
+ // Note that the file <code>include/deal.II/matrix_free/operators.h</code>
+ // already contains an implementation of the Laplacian through the class
+ // MatrixFreeOperators::LaplaceOperator. For educational purposes, the
+ // operator is re-implemented in this tutorial program, explaining the
+ // ingredients and concepts used there.
+ //
// This program makes use of the data cache for finite element operator
// application that is integrated in deal.II. This data cache class is
// called MatrixFree. It contains mapping information (Jacobians) and index
// 64-bit floating point) for the final matrix, but floats (single
// precision, 32-bit floating point numbers) for the multigrid level
// matrices (as that is only a preconditioner, and floats can be processed
- // twice as fast).
+ // twice as fast). The class FEEvaluation also takes a template argument for
+ // the number of quadrature points in one dimension. In the code below, we
+ // hard-code it to <code>fe_degree+1</code>. If we wanted to change it
+ // independently of the polynomial degree, we would need to add a template
+ // parameter as is done in the MatrixFreeOperators::LaplaceOperator class.
//
// As a sidenote, if we implemented several different operations on the same
// grid and degrees of freedom (like a mass matrix and a Laplace matrix), we
// would define two classes like the current one for each of the operators
- // (maybe with a common base class), and let both of them refer to the same
- // MatrixFree data cache from the general problem class. The interface
- // through MatrixFreeOperators::Base requires us to only provide a minimal
- // set of functions. This concept allows for writing complex application
- // codes with many matrix-free apply operations.
+ // (derived from the MatrixFreeOperators::Base class), and let both of them
+ // refer to the same MatrixFree data cache from the general problem
+ // class. The interface through MatrixFreeOperators::Base requires us to
+ // only provide a minimal set of functions. This concept allows for writing
+ // complex application codes with many matrix-free operations.
//
// @note Storing values of type <code>VectorizedArray<number></code>
// requires care: Here, we use the deal.II table class which is prepared to
// derivatives of order between zero and two. We only want gradients, no
// values and no second derivatives, so we set the function arguments to
// true in the gradient slot (second slot), and to false in the values slot
- // (first slot). There is also a third slot for the Hessian whose setting
- // defaults to false, so it needs not be given. Note that the FEEvaluation
+ // (first slot). There is also a third slot for the Hessian which is
+ // false by default, so it needs not be given. Note that the FEEvaluation
// class internally evaluates shape functions in an efficient way where one
// dimension is worked on at a time (using the tensor product form of shape
// functions and quadrature points as mentioned in the introduction). This
// apply_add() function, so we do not need to take further action here.
//
// When using the combination of MatrixFree and FEEvaluation in parallel
- // with MPI, there is one aspect to be careful about. This is the indexing
+ // with MPI, there is one aspect to be careful about — the indexing
// used for accessing the vector. For performance reasons, MatrixFree and
// FEEvaluation are designed to access vectors in MPI-local index space also
// when working with multiple processors. Working in local index space means
- // that no index translation needs to be performed at the vector access,
- // apart from the unavoidable indirect addressing. However, local index
- // spaces are ambiguous: While it is standard convention to access the
- // locally owned range of a vector with indices between 0 and the local
- // size, the case is involved for the ghosted entries. For the matrix-vector
- // product, only the indices appearing on locally owned cells (plus those
- // referenced via hanging node constraints) are necessary. However, in
- // deal.II we often set all degrees of freedom on ghosted elements as
- // ghosted vector entries. In that case, the MPI-local index of a ghosted
- // vector entry in general be different in the two possible ghost sets,
- // despite referring to the same global index. To avoid problems,
- // FEEvaluation checks that the partitioning of the vector used for the
- // matrix-vector product does indeed match with the partitioning of the
- // indices in MatrixFree by a check called
+ // that no index translation needs to be performed at the place the vector
+ // access happns, apart from the unavoidable indirect addressing. However,
+ // local index spaces are ambiguous: While it is standard convention to
+ // access the locally owned range of a vector with indices between 0 and the
+ // local size, the numbering is not so clear for the ghosted entries and
+ // somewhat arbitrary. For the matrix-vector product, only the indices
+ // appearing on locally owned cells (plus those referenced via hanging node
+ // constraints) are necessary. However, in deal.II we often set all the
+ // degrees of freedom on ghosted elements as ghosted vector entries, called
+ // the @ref GlossLocallyRelevantDof "locally relevant DoFs described in the
+ // glossary". In that case, the MPI-local index of a ghosted vector entry
+ // can in general be different in the two possible ghost sets, despite
+ // referring to the same global index. To avoid problems, FEEvaluation
+ // checks that the partitioning of the vector used for the matrix-vector
+ // product does indeed match with the partitioning of the indices in
+ // MatrixFree by a check called
// LinearAlgebra::distributed::Vector::partitioners_are_compatible. To
// facilitate things, the MatrixFreeOperators::Base class includes a
// mechanism to fit the ghost set to the correct layout. This happens in the
- // ghost region of the vector, so the ghost region might be different in
- // both the destination and source vector after a call to a vmult()
- // mehtod. This is legitimate because the ghost region of a distributed
- // deal.II vector is mutable. Vectors used in matrix-vector products must
- // not be ghosted upon entry of vmult() functions, so no information gets
- // lost.
+ // ghost region of the vector, so keep in mind that the ghost region might
+ // be modified in both the destination and source vector after a call to a
+ // vmult() method. This is legitimate because the ghost region of a
+ // distributed deal.II vector is a mutable section and filled on
+ // demand. Vectors used in matrix-vector products must not be ghosted upon
+ // entry of vmult() functions, so no information gets lost.
template <int dim, int fe_degree, typename number>
void
LaplaceOperator<dim,fe_degree,number>
LinearAlgebra::distributed::Vector<number> &inverse_diagonal =
this->inverse_diagonal_entries->get_vector();
this->data->initialize_dof_vector(inverse_diagonal);
- unsigned int dummy;
+ unsigned int dummy = 0;
this->data->cell_loop (&LaplaceOperator::local_compute_diagonal, this,
inverse_diagonal, dummy);
for (unsigned int i=0; i<inverse_diagonal.local_size(); ++i)
{
- Assert(inverse_diagonal.local_element(i) != 0.,
+ Assert(inverse_diagonal.local_element(i) > 0.,
ExcMessage("No diagonal entry in a positive definite operator "
"should be zero"));
inverse_diagonal.local_element(i) =
// an integral contribution of a constrained DoF to several other entries
// inside the distribute_local_to_global call, the vector interface used
// here does not exactly compute the diagonal entries, but lumps some
- // contributions located on the diagonal of the local matrix but would end
- // up on a off-diagonal position of the global matrix to the diagonal. The
- // result is correct up to discretization accuracy as
shown in <a
+ // contributions located on the diagonal of the local matrix that would end
+ // up in a off-diagonal position of the global matrix to the diagonal. The
+ // result is correct up to discretization accuracy as
explained in <a
// href="http://dx.doi.org/10.4208/cicp.101214.021015a">Kormann (2016),
- // section 5.3</a>, but not equal. In this tutorial program, no harm can
- // happen because the diagonal is only used for the multigrid level matrices
- // where no hanging node constraints appear.
+ // section 5.3</a>, but not mathematically equal. In this tutorial program,
+ // no harm can happen because the diagonal is only used for the multigrid
+ // level matrices where no hanging node constraints appear.
template <int dim, int fe_degree, typename number>
void
LaplaceOperator<dim,fe_degree,number>
- // @sect4{MGInterfaceMatrix class}
-
- // The adaptive multigrid realization in deal.II implements an approach
- // called local smoothing. This means that the smoothing on the finest level
- // only covers the local part of the mesh defined by the fixed (finest) grid
- // level and ignores parts of the computational domain where the terminal
- // cells are coarser than this level. As the method progresses to coarser
- // levels, more and more of the global mesh will be covered. At some coarser
- // level, the whole mesh will be covered. Since all level matrices in the
- // multigrid method cover a single level in the mesh, no hanging nodes
- // appear on the level matrices. At the interface between multigrid levels,
- // homogeneous Dirichlet boundary conditions are set while smoothing. When
- // the residual is transferred to the next coarser level, however, the
- // coupling over the multigrid interface needs to be taken into account.
- // This is done by the so-called interface (or edge) matrices that compute
- // the part of the residual that is missed by the level matrix with
- // homogeneous Dirichlet conditions. We refer to the @ref mg_paper
- // "Multigrid paper by Janssen and Kanschat" for more details.
- //
- // For the implementation of those interface matrices, most infrastructure
- // is already in place and provided by MatrixFreeOperators::Base through the
- // two multiplication routines vmult_interface_down() and
- // vmult_interface_up(). The only thing we have to do in this tutorial
- // program is to wrap those multiplication routines into a separate class
- // where they are accessed via the @p vmult() and @p Tvmult interface as
- // expected by the multigrid routines (that were originally written for
- // matrices, hence expecting those names). Note that the
- // vmult_interface_down is used during the restriction phase of the
- // multigrid V-cycle, whereas vmult_interface_up is used during the
- // prolongation phase.
- template <typename OperatorType>
- class MGInterfaceMatrix : public Subscriptor
- {
- public:
- typedef typename OperatorType::value_type value_type;
-
- void initialize (const OperatorType &operator_in)
- {
- this->mf_base_operator = &operator_in;
- }
-
- void vmult (LinearAlgebra::distributed::Vector<value_type> &dst,
- const LinearAlgebra::distributed::Vector<value_type> &src) const
- {
- mf_base_operator->vmult_interface_down(dst, src);
- }
-
- void Tvmult (LinearAlgebra::distributed::Vector<value_type> &dst,
- const LinearAlgebra::distributed::Vector<value_type> &src) const
- {
- mf_base_operator->vmult_interface_up(dst, src);
- }
-
- private:
- SmartPointer<const OperatorType> mf_base_operator;
- };
-
-
-
// @sect3{LaplaceProblem class}
// This class is based on the one in step-16. However, we replaced the
SystemMatrixType system_matrix;
MGConstrainedDoFs mg_constrained_dofs;
- MGLevelObject<ConstraintMatrix> mg_constraints;
MGLevelObject<MatrixFree<dim,float> > mg_mf_storage;
typedef LaplaceOperator<dim,degree_finite_element,float> LevelMatrixType;
MGLevelObject<LevelMatrixType> mg_matrices;
fe (degree_finite_element),
dof_handler (triangulation),
pcout (std::cout, Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0),
+ // The LaplaceProblem class holds an additional output stream that
+ // collects detailed timings about the setup phase. This stream, called
+ // time_details, is disabled by default through the @p false argument
+ // specified here. For detailed timings, removing the @p false argument
+ // prints all the details.
time_details (std::cout, false &&
Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
{}
// step-40.
// Once we have created the multigrid dof_handler and the constraints, we
- // can call the reinit function for each level of the multigrid routine (and
- // the active cells). The main purpose of the reinit function is to setup
- // the <code> MatrixFree </code> instance for the problem. Also, the
- // coefficient is evaluated. For this, we need to activate the update flag
- // in the AdditionalData field of MatrixFree that enables the storage of
- // quadrature point coordinates in real space (by default, it only caches
- // data for gradients (inverse transposed Jacobians) and JxW values). Note
- // that if we call the reinit function without specifying the level (i.e.,
- // giving <code>level = numbers::invalid_unsigned_int</code>), MatrixFree
- // constructs a loop over the active cells.
+ // can call the reinit function for the global matrix operator as well as
+ // each level of the multigrid scheme. The main action is to set up the
+ // <code> MatrixFree </code> instance for the problem. For this, we need to
+ // activate the update flag in the AdditionalData field of MatrixFree that
+ // enables the storage of quadrature point coordinates in real space (by
+ // default, it only caches data for gradients (inverse transposed Jacobians)
+ // and JxW values). Note that if we call the reinit function without
+ // specifying the level (i.e., giving <code>level =
+ // numbers::invalid_unsigned_int</code>), MatrixFree constructs a loop over
+ // the active cells. In this tutorial, we do not use threads in addition to
+ // MPI, which is why we explicitly disable it by setting the
+ // MatrixFree::AdditionalData::tasks_parallel_scheme to
+ // MatrixFree::AdditionalData::none. Finally, the coefficient is evaluated
+ // and vectors are initialized as explained above.
template <int dim>
void LaplaceProblem<dim>::setup_system ()
{
system_matrix.clear();
mg_matrices.clear_elements();
- mg_constraints.clear_elements();
dof_handler.distribute_dofs (fe);
dof_handler.distribute_mg_dofs (fe);
const unsigned int nlevels = triangulation.n_global_levels();
mg_matrices.resize(0, nlevels-1);
mg_mf_storage.resize(0, nlevels-1);
- mg_constraints.resize (0, nlevels-1);
std::set<types::boundary_id> dirichlet_boundary;
dirichlet_boundary.insert(0);
IndexSet relevant_dofs;
DoFTools::extract_locally_relevant_level_dofs(dof_handler, level,
relevant_dofs);
- mg_constraints[level].reinit(relevant_dofs);
- mg_constraints[level].add_lines(mg_constrained_dofs.get_boundary_indices(level));
- mg_constraints[level].close();
+ ConstraintMatrix level_constraints;
+ level_constraints.reinit(relevant_dofs);
+ level_constraints.add_lines(mg_constrained_dofs.get_boundary_indices(level));
+ level_constraints.close();
typename MatrixFree<dim,float>::AdditionalData additional_data;
additional_data.tasks_parallel_scheme =
additional_data.mpi_communicator = MPI_COMM_WORLD;
additional_data.level_mg_handler = level;
- mg_mf_storage[level].reinit(dof_handler, mg_constraints[level],
+ mg_mf_storage[level].reinit(dof_handler, level_constraints,
QGauss<1>(fe.degree+1), additional_data);
mg_matrices[level].initialize(mg_mf_storage[level], mg_constrained_dofs, level);
// to use the Chebyshev iteration as a solver. PreconditionChebyshev
// allows the user to switch to solver mode where the number of iterations
// is internally chosen to the correct value. In the additional data
- // object, this setting is activated by choosen the polynomial degree to
+ // object, this setting is activated by choosing the polynomial degree to
// @p numbers::invalid_unsigned_int. The algorithm will then attack all
// eigenvalues between the smallest and largest one in the coarse level
// matrix. The number of steps in the Chebyshev smoother are chosen such
// methods. The former involves only local communication between neighbors
// in the (coarse) mesh, whereas the latter requires global communication
// over all processors.
- typedef PreconditionChebyshev<LevelMatrixType,LinearAlgebra::distributed::Vector<float> > SMOOTHER;
- mg::SmootherRelaxation<SMOOTHER, LinearAlgebra::distributed::Vector<float> >
+ typedef PreconditionChebyshev<LevelMatrixType,LinearAlgebra::distributed::Vector<float> > SmootherType;
+ mg::SmootherRelaxation<SmootherType, LinearAlgebra::distributed::Vector<float> >
mg_smoother;
- MGLevelObject<typename SMOOTHER::AdditionalData> smoother_data;
+ MGLevelObject<typename SmootherType::AdditionalData> smoother_data;
smoother_data.resize(0, triangulation.n_global_levels()-1);
for (unsigned int level = 0; level<triangulation.n_global_levels(); ++level)
{
}
else
{
- smoother_data[level].smoothing_range = 1e-3;
+ smoother_data[0].smoothing_range = 1e-3;
smoother_data[0].degree = numbers::invalid_unsigned_int;
smoother_data[0].eig_cg_n_iterations = mg_matrices[0].m();
}
MGCoarseGridApplySmoother<LinearAlgebra::distributed::Vector<float> > mg_coarse;
mg_coarse.initialize(mg_smoother);
- // Next, we set up the interface matrices that are needed for the case
- // with hanging nodes. We have prepared the necessary class
- // MGInterfaceMatrix above for accessing the infrastructure of
- // MatrixFreeOperators::Base. All we have to do is to create a respective
- // MGLevelObject. Once that is done, we set up the remaining Multigrid
+ // The next step is to set up the interface matrices that are needed for the case
+ // with hanging nodes. The adaptive multigrid realization in deal.II implements an approach
+ // called local smoothing. This means that the smoothing on the finest level
+ // only covers the local part of the mesh defined by the fixed (finest) grid
+ // level and ignores parts of the computational domain where the terminal
+ // cells are coarser than this level. As the method progresses to coarser
+ // levels, more and more of the global mesh will be covered. At some coarser
+ // level, the whole mesh will be covered. Since all level matrices in the
+ // multigrid method cover a single level in the mesh, no hanging nodes
+ // appear on the level matrices. At the interface between multigrid levels,
+ // homogeneous Dirichlet boundary conditions are set while smoothing. When
+ // the residual is transferred to the next coarser level, however, the
+ // coupling over the multigrid interface needs to be taken into account.
+ // This is done by the so-called interface (or edge) matrices that compute
+ // the part of the residual that is missed by the level matrix with
+ // homogeneous Dirichlet conditions. We refer to the @ref mg_paper
+ // "Multigrid paper by Janssen and Kanschat" for more details.
+ //
+ // For the implementation of those interface matrices, there is already a
+ // pre-defined class MatrixFreeOperators::MGInterfaceOperator that wraps
+ // the routines MatrixFreeOperators::Base::vmult_interface_down() and
+ // MatrixFreeOperators::Base::vmult_interface_up() in a new class with @p
+ // vmult() and @p Tvmult() operations (that were originally written for
+ // matrices, hence expecting those names). Note that vmult_interface_down
+ // is used during the restriction phase of the multigrid V-cycle, whereas
+ // vmult_interface_up is used during the prolongation phase.
+ //
+ // Once the interface matrix created, we set up the remaining Multigrid
// preconditioner infrastructure in complete analogy to step-16 to obtain
// a @p preconditioner object that can be applied to a matrix.
mg::Matrix<LinearAlgebra::distributed::Vector<float> > mg_matrix(mg_matrices);
- MGLevelObject<MGInterfaceMatrix<LevelMatrixType> > mg_interface_matrices;
+ MGLevelObject<MatrixFreeOperators::MGInterfaceOperator<LevelMatrixType> > mg_interface_matrices;
mg_interface_matrices.resize(0, triangulation.n_global_levels()-1);
for (unsigned int level=0; level<triangulation.n_global_levels(); ++level)
mg_interface_matrices[level].initialize(mg_matrices[level]);
mg_transfer,
mg_smoother,
mg_smoother);
+ mg.set_edge_matrices(mg_interface, mg_interface);
+
PreconditionMG<dim, LinearAlgebra::distributed::Vector<float>,
MGTransferMatrixFree<dim,float> >
preconditioner(dof_handler, mg, mg_transfer);
// The setup of the multigrid routines was quite easy and one cannot see
- // any difference in the solve process as compared to step-16. The magic
- // is all hidden behind the implementation of the LaplaceOperator::vmult
+ // any difference in the solve process as compared to step-16. All the
+ // magic is hidden behind the implementation of the LaplaceOperator::vmult
// operation. Note that we print out the solve time and the accumulated
// setup time through standard out, i.e., in any case, whereas detailed
// times for the setup operations are only printed in case the flag for
projects / dealii.git / commitdiff