| Table of contents | |
|---|---|
This program was contributed by Katharina Kormann and Martin Kronbichler.
This program is currently under construction.
The algorithm for the matrix-vector product is built upon the report "MPI parallelization of a cell-based matrix-vector product for finite elements. An application from quantum dynamics" by Katharina Kormann, Uppsala University, June 2009.
This example shows how to implement a matrix-free method, that is, a method that does not explicitly store the matrix elements, for a second-order Poisson equation with variable coefficients on a fairly unstructured mesh representing a circle.
In order to find out how we can write a code that performs a matrix-vector product, but does not need to store the matrix elements, let us start at looking how some finite-element related matrix A is assembled:
In this formula, the matrix Pcell,loc-glob is a rectangular 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 local_dof_indices variable we have always used in the assembly of matrices.
If we are to perform a matrix-vector product, we can hence use that
where xcell are the values of x at the degrees of freedom of the respective cell, and xcell correspondingly for the result. A naive attempt to implement the local action of the Laplacian would hence be to use the following code:
MatrixFree::vmult (Vector<double> &dst, const Vector<double> &src) const { dst = 0; QGauss<dim> quadrature_formula(fe.degree+1); FEValues<dim> fe_values (fe, quadrature_formula, update_gradients | update_JxW_values); const unsigned int dofs_per_cell = fe.dofs_per_cell; const unsigned int n_q_points = quadrature_formula.size(); FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell); Vector<double> cell_src (dofs_per_cell), cell_dst (dofs_per_cell); std::vector<unsigned int> local_dof_indices (dofs_per_cell); typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end(); for (; cell!=endc; ++cell) { cell_matrix = 0; fe_values.reinit (cell); for (unsigned int q_point=0; q_point<n_q_points; ++q_point) for (unsigned int i=0; i<dofs_per_cell; ++i) for (unsigned int j=0; j<dofs_per_cell; ++j) cell_matrix(i,j) += (fe_values.shape_grad(i,q_point) * fe_values.shape_grad(j,q_point) * fe_values.JxW(q_point)); cell->get_dof_indices (local_dof_indices); for (unsigned int i=0; i<dofs_per_cell; ++i) cell_src(i) = src(local_dof_indices(i)); cell_matrix.vmult (cell_dst, cell_src); for (unsigned int i=0; i<dofs_per_cell; ++i) dst(local_dof_indices(i)) += cell_dst; } }
Here we neglected boundary conditions as well as any hanging nodes we may have, though neither would be very difficult to include using the ConstraintMatrix class. Note how we first generate the local matrix in the usual way. To form the actual product as expressed in the above formula, we read in the values of src of the cell-related degrees of freedom (the action of Pcell,loc-glob), multiply by the local matrix (the action of Acell), and finally add the result to the destination vector dst (the action of Pcell,loc-globT, added over all the elements). It is not more difficult than that, in principle.
While this code is completely correct, it is very slow. For every cell, we generate a local matrix, which takes three nested loops with as many elements as there are degrees of freedom on the actual cell to compute. The multiplication itself is then done by two nested loops, which means that it is much cheaper.
One way to improve this is to realize that conceptually the local matrix can be thought of as the product of three matrices,
where for the example of the Laplace operator the (q*dim+d,i)-th element of Bcell is given by fe_values.shape_grad(i,q)[d]. The matrix consists of dim*n_q_points rows and dofs_per_cell columns). The matrix Dcell is diagonal and contains the values fe_values.JxW(q) (or, rather, dim copies of it).
Every numerical analyst learns in one of her first classes that for forming a product of the form
one should never form the matrix-matrix products, but rather multiply with the vector from right to left so that only three successive matrix-vector products are formed. To put this into code, we can write:
... for (; cell!=endc; ++cell) { fe_values.reinit (cell); cell->get_dof_indices (local_dof_indices); for (unsigned int i=0; i<dofs_per_cell; ++i) cell_src(i) = src(local_dof_indices(i)); temp_vector = 0; for (unsigned int q_point=0; q_point<n_q_points; ++q_point) for (unsigned int d=0; d<dim; ++d) for (unsigned int i=0; i<dofs_per_cell; ++i) temp_vector(q_point*dim+d) += fe_values.shape_grad(i,q_point)[d] * cell_src(i); for (unsigned int q_point=0; q_point<n_q_points; ++q_point) for (unsigned int d=0; d<dim; ++d) temp_vector(q_point*dim+d) *= fe_values.JxW(q_point); cell_dst = 0; for (unsigned int i=0; i<dofs_per_cell; ++i) for (unsigned int q_point=0; q_point<n_q_points; ++q_point) for (unsigned int d=0; d<dim; ++d) cell_dst(i) += fe_values.shape_grad(i,q_point)[d] * temp_vector(q_point*dim+d); for (unsigned int i=0; i<dofs_per_cell; ++i) dst(local_dof_indices(i)) += cell_dst(i); } }
This removed the three nested loops in the calculation of the local matrix (here the loop over d is a not really a loop, rather two or three operations). What happens is as follows: We first transform the vector of values on the local dofs to a vector of gradients on the quadrature points. In the second loop, we multiply these gradients by the integration weight. The third loop applies the second gradient (in transposed form), so that we get back to a vector of (Laplacian) values on the cell dofs.
This improves the situation a lot and reduced the complexity of the product from something like
to
. In fact, all the remainder is just to make a slightly more clever use of data in order to gain some extra speed. It does not change the code structure, though.
The bottleneck in the above code is the operations done by the call fe_values.reinit(cell), which take about as much time as the other steps together (at least if the mesh is unstructured; deal.II can recognize that the gradients are often unchanged on structured meshes). That is certainly not ideal and we would like to do better than this. What the reinit function does is to calculate the gradient in real space by transforming the gradient on the reference cell using the Jacobian of the transformation from real to reference cell. This is done for each basis function on the cell, for each quadrature point. The Jacobian does not depend on the basis function, but it is different on different quadrature points in general. The trick is now to factor out the Jacobian transformation and first apply the operation that leads us to temp_vector only with the gradient on the reference cell. 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, then we apply the weights of the quadrature, and we apply with the transposed Jacobian for preparing the third loop which again uses the gradients on the unit cell.
Let us again write this in terms of matrices. Let the matrix Bcell denote the cell-related gradient matrix, with each row containing the values of the quadrature points. It is constructed by a matrix-matrix product as
where Bref_cell denotes the gradient on the reference cell and Jcell denotes the Jacobian transformation. Jcell is block-diagonal, and the blocks size is equal to the dimension of the problem. Each diagonal block is the Jacobian transformation that goes from the reference cell to the real cell.
Putting things together, we find that
so we calculate the product (starting the local product from the right)
... FEValues<dim> fe_values_reference (fe, quadrature_formula, update_gradients); Triangulation<dim> reference_cell; GridGenerator::hyper_cube(reference_cell, 0., 1.); fe_values_reference.reinit (reference_cell.begin()); FEValues<dim> fe_values (fe, quadrature_formula, update_inverse_jacobians | update_JxW_values); for (; cell!=endc; ++cell) { fe_values.reinit (cell); cell->get_dof_indices (local_dof_indices); for (unsigned int i=0; i<dofs_per_cell; ++i) cell_src(i) = src(local_dof_indices(i)); temp_vector = 0; for (unsigned int q_point=0; q_point<n_q_points; ++q_point) for (unsigned int d=0; d<dim; ++d) for (unsigned int i=0; i<dofs_per_cell; ++i) temp_vector(q_point*dim+d) += fe_values_reference.shape_grad(i,q_point)[d] * cell_src(i); for (unsigned int q_point=0; q_point<n_q_points; ++q_point) { // apply the Jacobian of the mapping from unit to real cell Tensor<1,dim> temp; for (unsigned int d=0; d<dim; ++d) temp[d] = temp_vector(q_point*dim+d); for (unsigned int d=0; d<dim; ++d) { double sum = 0; for (unsigned int e=0; e<dim; ++e) sum += fe_values.inverse_jacobian(q_point)[d][e] * temp[e]; temp_vector(q_point*dim+d) = sum; } // multiply with integration weight for (unsigned int d=0; d<dim; ++d) temp_vector(q_point*dim+d) *= fe_values.JxW(q_point); // apply the transpose of the Jacobian of the mapping from unit // to real cell for (unsigned int d=0; d<dim; ++d) temp[d] = temp_vector(q_point*dim+d); for (unsigned int d=0; d<dim; ++d) { double sum = 0; for (unsigned int e=0; e<dim; ++e) sum += fe_values.inverse_jacobian(q_point)[e][d] * temp[e]; temp_vector(q_point*dim+d) = sum; } } cell_dst = 0; for (unsigned int i=0; i<dofs_per_cell; ++i) for (unsigned int q_point=0; q_point<n_q_points; ++q_point) for (unsigned int d=0; d<dim; ++d) cell_dst(i) += fe_values_reference.shape_grad(i,q_point)[d] * temp_vector(q_point*dim+d); for (unsigned int i=0; i<dofs_per_cell; ++i) dst(local_dof_indices(i)) += cell_dst(i); } }
Note how we create an additional FEValues object for the reference cell gradients and how we initialize it to the reference cell. The actual derivative data is then applied by the Jacobians (deal.II calls the Jacobian matrix from unit to real cell inverse_jacobian, because the transformation direction in deal.II is from real to unit cell).
To sum up, we want to look at the additional costs introduced by this written-out matrix-matrix-vector product compared to a sparse matrix. To first approximation, we have increased the number of operations for the local matrix-vector product by a factor of 4 in 2D and 6 in 3D, i.e., two matrix-vector products with matrices that have dim as many columns than before (here we assume that the number of quadrature points is the same as the number of degrees of freedom per cell, which is usual for scalar problems). Then, we also need to keep in mind that we touch some degrees of freedom several times because they belong to several cells. This also increases computational costs. A realistic value compared to a sparse matrix is that we now have to perform about 10 times as many operations (a bit less in 2D, a bit more in 3D).
The above is, in essence, what happens in the code below and if you have difficulties in understanding the implementation, you should try to first understand what happens in the code above. In the actual implementation there are a few more points done to be even more efficient, namely:
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 linear and quadratic elements, but on par with third order elements and faster for even higher order elements. 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.
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 uncertain 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. On the other hand, most of the more frequently used preconditioners such as SSOR, ILU or algebraic multigrid (AMG) can now no longer be used here because their implementation requires knowledge of the elements of the system matrix.
One solution is to use 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. 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 preconditioner, eventually, is what we use here. 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 M. Adams, M. Brezina, J. Hu, R. Tuminaro. Parallel multigrid smoothers: polynomial versus Gauss–Seidel, J. Comput. Phys. 188:593–610, 2003. 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, which can often only be parallelized by working on diagonal sub-blocks of the matrix, which decreases efficiency.
The implementation into the multigrid framework is then straightforward. This program is based on an earlier version of step-16 that demonstrated multigrid on uniformly refined grids. However, the present matrix-free techniques would obviously also apply to the adaptive meshes the current step-16 uses.
In order to demonstrate the capabilities of the method, we work on a rather general Poisson problem, based on a more or less unstructured mesh (where the Jacobians are different from cell to cell), higher order mappings to a curved boundary, and a non-constant coefficient in the equation. If we worked on a constant-coefficient case with structured mesh, we could decrease the operation count by a factor of 4 in 2D and 6 in 3D by building a local matrix (which is then the same for all cells), and doing the products as in the first developing step of the above code pieces.
To start with the include files are more or less the same as in step-16:
#include <base/quadrature_lib.h> #include <base/function.h> #include <base/logstream.h> #include <base/work_stream.h> #include <lac/vector.h> #include <lac/full_matrix.h> #include <lac/solver_cg.h> #include <lac/precondition.h> #include <fe/fe_q.h> #include <fe/fe_values.h> #include <grid/tria.h> #include <grid/tria_accessor.h> #include <grid/tria_iterator.h> #include <grid/tria_boundary_lib.h> #include <grid/grid_generator.h> #include <multigrid/multigrid.h> #include <multigrid/mg_dof_handler.h> #include <multigrid/mg_dof_accessor.h> #include <multigrid/mg_transfer.h> #include <multigrid/mg_tools.h> #include <multigrid/mg_coarse.h> #include <multigrid/mg_smoother.h> #include <multigrid/mg_matrix.h> #include <numerics/data_out.h> #include <numerics/vectors.h> #include <fstream> #include <sstream> using namespace dealii;
We define a variable coefficient function for the Poisson problem. It is similar to the function in step-5 but we use the form
instead of a discontinuous one. It is merely to demonstrate the possibilities of this implementation, rather than making much sense physically.
template <int dim> class Coefficient : public Function<dim> { public: Coefficient () : Function<dim>() {} virtual double value (const Point<dim> &p, const unsigned int component = 0) const; virtual void value_list (const std::vector<Point<dim> > &points, std::vector<double> &values, const unsigned int component = 0) const; }; template <int dim> double Coefficient<dim>::value (const Point<dim> &p, const unsigned int / *component* /) const { return 1./(0.1+p.square()); } template <int dim> void Coefficient<dim>::value_list (const std::vector<Point<dim> > &points, std::vector<double> &values, const unsigned int component) const { Assert (values.size() == points.size(), ExcDimensionMismatch (values.size(), points.size())); Assert (component == 0, ExcIndexRange (component, 0, 1)); const unsigned int n_points = points.size(); for (unsigned int i=0; i<n_points; ++i) values[i] = 1./(0.1+points[i].square()); }
In this program, we want to make use of the ability of deal.II to runs things in parallel if compute resources are available. We will follow the general framework laid out in the Parallel computing with multiple processors accessing shared memory module and use the WorkStream class to do operations on the range of all cells.
To this end, we first have to have a few declarations that we use for defining the parallel layout of the vector multiplication function with the WorkStream concept in the Matrix-free class. These comprise so-called scratch data that we use for calculating cell-related information, and copy data that is eventually used in a separate function for writing local data into the global vector. The reason for this split-up definition is that many threads at a time can execute the local multiplications (and filling up the copy data), but than that copy data needs to be worked on by one process at a time.
namespace WorkStreamData { template <typename number> struct ScratchData { ScratchData (); ScratchData (const ScratchData &scratch); FullMatrix<number> solutions; }; template<typename number> ScratchData<number>::ScratchData () : solutions () {} template<typename number> ScratchData<number>::ScratchData (const ScratchData &) : solutions () {} template <typename number> struct CopyData : public ScratchData<number> { CopyData (); CopyData (const CopyData &scratch); unsigned int first_cell; unsigned int n_dofs; }; template <typename number> CopyData<number>::CopyData () : ScratchData<number> () {} template <typename number> CopyData<number>::CopyData (const CopyData &) : ScratchData<number> () {} }
Next comes the implementation of the matrix-free class. It provides some standard information we expect for matrices (like returning the dimensions of the matrix), it implements matrix-vector multiplications in several forms, and it provides functions for filling the matrix with data.
We choose to make this class generic, i.e., we do not implement the actual differential operator (here: Laplace operator) directly in this class. We instead let the actual transformation (which happens on the level of quadrature points, see the discussion in the introduction) be a template parameter that is implemented by another class. We then only have to store a list of these objects for each quadrature point on each cell in a big list – we choose a Table<2,Transformation> data format) – and call a transform command of the Transformation class. This template magic makes it easy to reuse this MatrixFree class for other problems that are based on a symmetric operation without the need for substantial changes.
template <typename number, class Transformation> class MatrixFree : public Subscriptor { public: MatrixFree (); void reinit (const unsigned int n_dofs, const unsigned int n_cells, const FullMatrix<double> &cell_matrix, const unsigned int n_points_per_cell); void clear(); unsigned int m () const; unsigned int n () const; ConstraintMatrix & get_constraints (); void set_local_dof_indices (const unsigned int cell_no, const std::vector<unsigned int> &local_dof_indices); void set_derivative_data (const unsigned int cell_no, const unsigned int quad_point, const Transformation &trans_in); template <typename number2> void vmult (Vector<number2> &dst, const Vector<number2> &src) const; template <typename number2> void Tvmult (Vector<number2> &dst, const Vector<number2> &src) const; template <typename number2> void vmult_add (Vector<number2> &dst, const Vector<number2> &src) const; template <typename number2> void Tvmult_add (Vector<number2> &dst, const Vector<number2> &src) const; number el (const unsigned int row, const unsigned int col) const; void calculate_diagonal () const; std::size_t memory_consumption () const;
The private member variables of the MatrixFree class are a small matrix that does the transformation from solution values to quadrature points, a list with the mapping between local degrees of freedom and global degrees of freedom for each cell (stored as a two-dimensional array, where each row corresponds to one cell, and the columns within individual cells are the local degrees of freedom), the transformation variable for implementing derivatives, a constraint matrix for handling boundary conditions as well as a few other variables that store matrix properties.
private: typedef std::vector<std::pair<unsigned int,unsigned int> >::const_iterator CellChunkIterator; template <typename number2> void local_vmult (CellChunkIterator cell_range, WorkStreamData::ScratchData<number> &scratch, WorkStreamData::CopyData<number> ©, const Vector<number2> &src) const; template <typename number2> void copy_local_to_global (const WorkStreamData::CopyData<number> ©, Vector<number2> &dst) const; FullMatrix<number> B_ref_cell; Table<2,unsigned int> indices_local_to_global; Table<2,Transformation> derivatives; ConstraintMatrix constraints; mutable Vector<number> diagonal_values; mutable bool diagonal_is_calculated; struct MatrixSizes { unsigned int n_dofs, n_cells; unsigned int m, n; unsigned int n_points, n_comp; std::vector<std::pair<unsigned int,unsigned int> > chunks; } matrix_sizes; };
This is the constructor of the MatrixFree class. All it does is to subscribe to the general deal.II Subscriptor scheme that makes sure that we do not delete an object of this class as long as it used somewhere else, e.g. in a preconditioner.
template <typename number, class Transformation> MatrixFree<number,Transformation>::MatrixFree () : Subscriptor() {}
The next functions return the number of rows and columns of the global matrix (i.e. the dimensions of the operator this class represents, the point of this tutorial program was, after all, that we don't actually store the elements of the rows and columns of this operator). Since the matrix is square, the returned numbers are the same.
template <typename number, class Transformation> unsigned int MatrixFree<number,Transformation>::m () const { return matrix_sizes.n_dofs; } template <typename number, class Transformation> unsigned int MatrixFree<number,Transformation>::n () const { return matrix_sizes.n_dofs; }
One more function that just returns an internal variable. Note that the user will need to change this variable, so it returns a non-constant reference to the ConstraintMatrix.
template <typename number, class Transformation> ConstraintMatrix & MatrixFree<number,Transformation>::get_constraints () { return constraints; }
The following function takes a vector of local dof indices on cell level and writes the data into the indices_local_to_global field in order to have fast access to it. It performs a few sanity checks like whether the sizes in the matrix are set correctly. One tiny thing: Whenever we enter this function, we probably make some modification to the matrix. This means that the diagonal of the matrix, which we might have computed to have fast access to those elements, is invalidated. We set the respective flag to false.
template <typename number, class Transformation> void MatrixFree<number,Transformation>:: set_local_dof_indices (const unsigned int cell_no, const std::vector<unsigned int> &local_dof_indices) { Assert (local_dof_indices.size() == matrix_sizes.m, ExcDimensionMismatch(local_dof_indices.size(), matrix_sizes.m)); for (unsigned int i=0; i<matrix_sizes.m; ++i) { Assert (local_dof_indices[i] < matrix_sizes.n_dofs, ExcInternalError()); indices_local_to_global(cell_no,i) = local_dof_indices[i]; } diagonal_is_calculated = false; }
Next a function that writes the derivative data on a certain cell and a certain quadrature point to the array that keeps the data around. Even though the array derivatives stands for the majority of the matrix memory consumption, it still pays off to have that data around since it would be quite expensive to manually compute it every time we make a matrix-vector product.
template <typename number, class Transformation> void MatrixFree<number,Transformation>:: set_derivative_data (const unsigned int cell_no, const unsigned int quad_point, const Transformation &trans_in) { Assert (quad_point < matrix_sizes.n_points, ExcInternalError()); derivatives(cell_no,quad_point) = trans_in; diagonal_is_calculated = false; }
Now finally to the central function of the matrix-free class, implementing the multiplication of the matrix with a vector. This function does not actually work on all cells of a mesh, but only the subset of cells specified by the first argument cell_range. Since this function operates similarly irrespective on which cell chunk we are sitting, we can call it simultaneously on many processors, but with different cell range data.
The goal of this function is to provide the multiplication of a vector with the local contributions of a set of cells. As mentioned in the introduction, if we were to deal with a single cell, this would amount to performing the product
where
and Pcell,local-global is the transformation from local to global indices.
To do this, we would have to do the following steps:
. This is done by using the command ConstraintMatrix::get_dof_values.
. The vector x1 contains the reference cell gradient to the local cell vector.
. This is a block-diagonal operation, with the block size equal to dim. The blocks just correspond to the individual quadrature points. The operation on each quadrature point is implemented by the Transformation class object that this class is equipped with. Compared to the introduction, the matrix Dcell now contains the JxW values and the inhomogeneous coefficient.
. This gives the local result of the matrix-vector product.
. This adds the local result to the global vector, which is realized using the method ConstraintMatrix::distribute_local_to_global. Note that we do this in an extra function called copy_local_to_global because that operation must not be done in parallel, in order to avoid two or more processes trying to add to the same positions in the result vector y.
Now, it turns out that the most expensive part of the above is the multiplication Bref_cell xcell in the second step and the transpose operation in step 4. Note that the matrix JT D J is block-diagonal, and hence, its application is cheaper. Since the matrix Bref_cell is the same for all cells, all that changes is the vector xcell. Hence, nothing prevents us from collecting several cell vectors to a (rectangular) matrix, and then perform a matrix-matrix product. These matrices are both full, but not very large, having of the order dofs_per_cell rows and columns. This is an operation that can be much better optimized than matrix-vector products. The functions FullMatrix<number>::mmult and FullMatrix<number>::mTmult use the BLAS dgemm function (as long as BLAS has been detected in deal.II configuration), which provides optimized kernels for doing this product. In our case, a matrix-matrix product is between three and five times faster than doing the matrix-vector product on one cell after the other. The variables that hold the solution on the respective cell's support points and the quadrature points are thus full matrices, which we set to the correct size as a first action in this function. The number of rows in the two matrices scratch.solutions and copy.solutions is given by the number of cells they work on, and the number of columns is the number of degrees of freedom per cell for the first and the number of quadrature points times the number of components per point for the latter.
template <typename number, class Transformation> template <typename number2> void MatrixFree<number,Transformation>:: local_vmult (CellChunkIterator cell_range, WorkStreamData::ScratchData<number> &scratch, WorkStreamData::CopyData<number> ©, const Vector<number2> &src) const { const unsigned int chunk_size = cell_range->second - cell_range->first; scratch.solutions.reinit (chunk_size, matrix_sizes.n, true); copy.solutions.reinit (chunk_size, matrix_sizes.m, true); copy.first_cell = cell_range->first; copy.n_dofs = chunk_size*matrix_sizes.m; constraints.get_dof_values(src, &indices_local_to_global(copy.first_cell,0), ©.solutions(0,0), ©.solutions(0,0)+copy.n_dofs); copy.solutions.mmult (scratch.solutions, B_ref_cell); for (unsigned int i=0, k = copy.first_cell; i<chunk_size; ++i, ++k) for (unsigned int j=0; j<matrix_sizes.n_points; ++j) derivatives(k,j).transform(&scratch.solutions(i, j*matrix_sizes.n_comp)); scratch.solutions.mTmult (copy.solutions, B_ref_cell); } template <typename number, class Transformation> template <typename number2> void MatrixFree<number,Transformation>:: copy_local_to_global (const WorkStreamData::CopyData<number> ©, Vector<number2> &dst) const { constraints.distribute_local_to_global (©.solutions(0,0), ©.solutions(0,0)+copy.n_dofs, &indices_local_to_global(copy.first_cell,0), dst); }
Now to the vmult function that is called externally: In addition to what we do in a vmult_add function, we set the destination to zero first.
template <typename number, class Transformation> template <typename number2> void MatrixFree<number,Transformation>::vmult (Vector<number2> &dst, const Vector<number2> &src) const { dst = 0; vmult_add (dst, src); }
Transposed matrix-vector products (needed for the multigrid operations to be well-defined): do the same. Since we implement a symmetric operation, we can refer to the vmult_add operation.
template <typename number, class Transformation> template <typename number2> void MatrixFree<number,Transformation>::Tvmult (Vector<number2> &dst, const Vector<number2> &src) const { dst = 0; Tvmult_add (dst,src); } template <typename number, class Transformation> template <typename number2> void MatrixFree<number,Transformation>::Tvmult_add (Vector<number2> &dst, const Vector<number2> &src) const { vmult_add (dst,src); }
This is the vmult_add function that multiplies the matrix with vector src and adds the result to vector dst. We include a few sanity checks to make sure that the size of the vectors is the same as the dimension of the matrix. We call a parallel function that applies the multiplication on a chunk of cells at once using the WorkStream module (cf. also the Parallel computing with multiple processors accessing shared memory module). The subdivision into chunks will be performed in the reinit function and is stored in the field matrix_sizes.chunks. What the rather cryptic command to std_cxx1x::bind does is to transform a function that has several arguments (source vector, chunk information) into a function which has three arguments (in the first case) or one argument (in the second), which is what the WorkStream::run function expects. The placeholders _1, _2, _3 in the local vmult specify variable input values, given by the chunk information, scratch data and copy data that the WorkStream::run function will provide, whereas the other arguments to the local_vmult function are bound: to this and a constant reference to the src in the first case, and this and a reference to the output vector in the second. Similarly, the placeholder _1 argument in the copy_local_to_global function sets the first explicit argument of that function, which is of class CopyData. We need to abstractly specify these arguments because the tasks defined by different cell chunks will be scheduled by the WorkStream class, and we will reuse available scratch and copy data.
template <typename number, class Transformation> template <typename number2> void MatrixFree<number,Transformation>::vmult_add (Vector<number2> &dst, const Vector<number2> &src) const { Assert (src.size() == n(), ExcDimensionMismatch(src.size(), n())); Assert (dst.size() == m(), ExcDimensionMismatch(dst.size(), m())); WorkStream::run (matrix_sizes.chunks.begin(), matrix_sizes.chunks.end(), std_cxx1x::bind(&MatrixFree<number,Transformation>:: template local_vmult<number2>, this, _1, _2, _3, boost::cref(src)), std_cxx1x::bind(&MatrixFree<number,Transformation>:: template copy_local_to_global<number2>, this, _1, boost::ref(dst)), WorkStreamData::ScratchData<number>(), WorkStreamData::CopyData<number>(), 2*multithread_info.n_default_threads,1);
One thing to be cautious about: The deal.II classes expect that the matrix still contains a diagonal entry for constrained dofs (otherwise, the matrix would be singular, which is not what we want). Since the distribute_local_to_global command of the constraint matrix which we used for adding the local elements into the global vector does not do anything with constrained elements, we have to circumvent that problem by artificially setting the diagonal to some non-zero value and adding the source values. We simply set it to one, which corresponds to copying the respective elements of the source vector into the matching entry of the destination vector.
for (unsigned int i=0; i<matrix_sizes.n_dofs; ++i) if (constraints.is_constrained(i) == true) dst(i) += 1.0 * src(i); }
The next function initializes the structures of the matrix. It writes the number of total degrees of freedom in the problem as well as the number of cells to the MatrixSizes struct and copies the small matrix that transforms the solution from support points to quadrature points. It uses the small matrix for determining the number of degrees of freedom per cell (number of rows in B_ref_cell). The number of quadrature points needs to be passed through the last variable n_points_per_cell, since the number of columns in the small matrix is dim*n_points_per_cell for the Laplace problem (the Laplacian is a tensor and has dim components). In this function, we also give the fields containing the derivative information and the local dof indices the correct sizes. They will be filled by calling the respective set function defined above.
template <typename number, class Transformation> void MatrixFree<number,Transformation>:: reinit (const unsigned int n_dofs_in, const unsigned int n_cells_in, const FullMatrix<double> &B_ref_cell_in, const unsigned int n_points_per_cell) { B_ref_cell = B_ref_cell_in; derivatives.reinit (n_cells_in, n_points_per_cell); indices_local_to_global.reinit (n_cells_in, B_ref_cell.m()); diagonal_is_calculated = false; matrix_sizes.n_dofs = n_dofs_in; matrix_sizes.n_cells = n_cells_in; matrix_sizes.m = B_ref_cell.m(); matrix_sizes.n = B_ref_cell.n(); matrix_sizes.n_points = n_points_per_cell; matrix_sizes.n_comp = B_ref_cell.n()/matrix_sizes.n_points; Assert(matrix_sizes.n_comp * n_points_per_cell == B_ref_cell.n(), ExcInternalError());
One thing to make the matrix-vector product with this class efficient is to decide how many cells should be combined to one chunk, which will determine the size of the full matrix that we work on. If we choose too few cells, then the gains from using the matrix-matrix product will not be fully utilized (dgemm tends to provide more efficiency the larger the matrix dimensions get), so we choose at least 60 cells for one chunk (except when there are very few cells, like on the coarse levels of the multigrid scheme). If we choose too many, we will degrade parallelization (we need to have sufficiently independent tasks). We need to also think about the fact that most high performance BLAS implementations internally work with square sub-matrices. Choosing as many cells in a chunk as there are degrees of freedom on each cell (coded in matrix_sizes.m) respects the BLAS GEMM design, whenever we exceed 60. Clearly, the chunk size is an architecture-dependent value and the interested user can squeeze out some extra performance by hand-tuning this parameter. Once we have chosen the number of cells we collect in one chunk, we determine how many chunks we have on the given cell range and recalculate the actual chunk size in order to evenly distribute the chunks.
const unsigned int divisor = std::max(60U, matrix_sizes.m); const unsigned int n_chunks = std::max (matrix_sizes.n_cells/divisor + 1, 2*multithread_info.n_default_threads); const unsigned int chunk_size = (matrix_sizes.n_cells/n_chunks + (matrix_sizes.n_cells%n_chunks>0)); std::pair<unsigned int, unsigned int> chunk; for (unsigned int i=0; i<n_chunks; ++i) { chunk.first = i*chunk_size; if ((i+1)*chunk_size > matrix_sizes.n_cells) chunk.second = matrix_sizes.n_cells; else chunk.second = (i+1)*chunk_size; if (chunk.second > chunk.first) matrix_sizes.chunks.push_back(chunk); else break; } }
Then we need a function if we want to delete the content of the matrix, e.g. when we are finished with one grid level and continue to the next one. Just set all the field sizes to 0.
template <typename number, class Transformation> void MatrixFree<number,Transformation>::clear () { B_ref_cell.reinit(0,0); derivatives.reinit (0,0); indices_local_to_global.reinit(0,0); constraints.clear(); diagonal_values.reinit (0); diagonal_is_calculated = false; matrix_sizes.n_dofs = 0; matrix_sizes.n_cells = 0; matrix_sizes.chunks.clear(); }
The next function returns the entries of the matrix. Since this class is intended not to store the matrix entries, it would make no sense to provide all those elements. However, diagonal entries are explicitly needed for the implementation of the Chebyshev smoother that we intend to use in the multigrid preconditioner. This matrix is equipped with a vector that stores the diagonal, and we compute it when this function is called for the first time.
template <typename number, class Transformation> number MatrixFree<number,Transformation>::el (const unsigned int row, const unsigned int col) const { Assert (row == col, ExcNotImplemented()); if (diagonal_is_calculated == false) calculate_diagonal(); return diagonal_values(row); }
Regarding the calculation of the diagonal, remember that this is as simple (or complicated) as assembling a right hand side in deal.II. Well, it is a bit easier to do this within this class since we have all the derivative information available. What we do is to go through all the cells (now in serial, since this function should not be called very often anyway), then all the degrees of freedom. At this place, we first copy the first basis functions in all the quadrature points to a temporary array, apply the derivatives from the Jacobian matrix, and finally multiply with the second basis function. This is exactly the value that would be written into the diagonal of a sparse matrix. Note that we need to condense hanging node constraints and set the constrained diagonals to one.
template <typename number, class Transformation> void MatrixFree<number,Transformation>::calculate_diagonal() const { diagonal_values.reinit (matrix_sizes.n_dofs); std::vector<number> calculation (matrix_sizes.n); for (unsigned int cell=0; cell<matrix_sizes.n_cells; ++cell) for (unsigned int dof=0; dof<matrix_sizes.m; ++dof) { memcpy (&calculation[0],&B_ref_cell(dof,0), matrix_sizes.n*sizeof(number)); for (unsigned int q=0; q<matrix_sizes.n_points; ++q) derivatives(cell,q).transform(&calculation[q*matrix_sizes.n_comp]); double diag_value = 0; for (unsigned int q=0; q<matrix_sizes.n; ++q) diag_value += calculation[q] * B_ref_cell(dof,q); diagonal_values(indices_local_to_global(cell,dof)) += diag_value; } constraints.condense (diagonal_values); for (unsigned int i=0; i<matrix_sizes.n_dofs; ++i) if (constraints.is_constrained(i) == true) diagonal_values(i) = 1.0; diagonal_is_calculated = true; }
Eventually, we provide a function that calculates how much memory this class uses. We just need to sum up the memory consumption of the arrays, the constraints, the small matrix and of the local variables. Just as a remark: In 2D and with data type double, about 80 per cent of the memory consumption is due to the derivatives array, while in 3D this number is even 85 per cent.
template <typename number, class Transformation> std::size_t MatrixFree<number,Transformation>::memory_consumption () const { std::size_t glob_size = derivatives.memory_consumption() + indices_local_to_global.memory_consumption() + constraints.memory_consumption() + B_ref_cell.memory_consumption() + diagonal_values.memory_consumption() + matrix_sizes.chunks.size()*2*sizeof(unsigned int) + sizeof(*this); return glob_size; }
This class implements the local action of a Laplace operator on a quadrature point. This is a very basic class implementation, providing functions for initialization with a Tensor of rank 2 and implementing the transform operation needed by the MatrixFree class. There is one point worth noting: The quadrature-point related action of the Laplace operator is a tensor of rank two. It is symmetric since it is the product of the inverse Jacobian transformation between unit and real cell with its transpose (times quadrature weights and a coefficient, which are scalar), so we can just save the diagonal and upper diagonal part. We could use the SymmetricTensor<2,dim> class for doing this, however, that class is only based on double numbers. Since we also want to use float numbers for the multigrid preconditioner (in order to save memory and computing time), we manually implement this operator. Note that dim is a template argument and hence known at compile-time, so the compiler knows that this symmetric rank-2 tensor has 3 entries if used in 2D and 6 entries if used in 3D.
template <int dim,typename number> class LaplaceOperator { public: LaplaceOperator (); LaplaceOperator (const Tensor<2,dim> &tensor); void transform (number * result) const; LaplaceOperator<dim,number>& operator = (const Tensor<2,dim> &tensor); private: number transformation[dim*(dim+1)/2]; }; template<int dim,typename number> LaplaceOperator<dim,number>::LaplaceOperator() {} template<int dim,typename number> LaplaceOperator<dim,number>::LaplaceOperator(const Tensor<2,dim> &tensor) { *this = tensor; }
Now implement the transformation, which is just a so-called contraction operation between a tensor of rank two and a tensor of rank one. Unfortunately, we need to implement this by hand, since we chose not to use the SymmetricTensor<2,dim> class (note that the resulting values are entries in a full matrix that consists of doubles or floats). It feels a bit unsafe to operate on a pointer to the data, but that is the only possibility if we do not want to copy data back and forth, which is expensive since this is the innermost position of the loop in the vmult operation of the MatrixFree class. We need to pay attention to the fact that we only saved half of the (symmetric) rank-two tensor.
At first sight, it seems inefficient that we have an if clause at this position in the code at the innermost loop, but note once again that dim is known when this piece of code is compiled, so the compiler can optimize away the if statement (and actually even inline these few lines of code into the MatrixFree class).
template <int dim, typename number> void LaplaceOperator<dim,number>::transform (number* result) const { if (dim == 2) { const number temp = result[0]; result[0] = transformation[0] * temp + transformation[1] * result[1]; result[1] = transformation[1] * temp + transformation[2] * result[1]; } else if (dim == 3) { const number temp1 = result[0]; const number temp2 = result[1]; result[0] = transformation[0] * temp1 + transformation[1] * temp2 + transformation[2] * result[2]; result[1] = transformation[1] * temp1 + transformation[3] * temp2 + transformation[4] * result[2]; result[2] = transformation[2] * temp1 + transformation[4] * temp2 + transformation[5] * result[2]; } else ExcNotImplemented(); }
The final function in this group takes the content of a rank-2 tensor and writes it to the field transformation of this class. We save the upper part of the symmetric tensor row-wise: we first take the (0,0)-entry, then the (0,1)-entry, and so on. We only implement this for dimensions two and three, which for the moment should do just fine:
template <int dim, typename number> LaplaceOperator<dim,number>& LaplaceOperator<dim,number>::operator=(const Tensor<2,dim> &tensor) { if (dim == 2) { transformation[0] = tensor[0][0]; transformation[1] = tensor[0][1]; transformation[2] = tensor[1][1]; Assert (std::fabs(tensor[1][0]-tensor[0][1])<1e-15, ExcInternalError()); } else if (dim == 3) { transformation[0] = tensor[0][0]; transformation[1] = tensor[0][1]; transformation[2] = tensor[0][2]; transformation[3] = tensor[1][1]; transformation[4] = tensor[1][2]; transformation[5] = tensor[2][2]; Assert (std::fabs(tensor[1][0]-tensor[0][1])<1e-15, ExcInternalError()); Assert (std::fabs(tensor[2][0]-tensor[0][2])<1e-15, ExcInternalError()); Assert (std::fabs(tensor[2][1]-tensor[1][2])<1e-15, ExcInternalError()); } else ExcNotImplemented(); return *this; }
This class is based on the same class in step-16. However, we replaced the SparseMatrix<double> class by our matrix-free implementation, which means that we can also skip the sparsity patterns.
template <int dim> class LaplaceProblem { public: LaplaceProblem (const unsigned int degree); void run (); private: void setup_system (); void assemble_system (); void assemble_multigrid (); void solve (); void output_results (const unsigned int cycle) const; Triangulation<dim> triangulation; FE_Q<dim> fe; MGDoFHandler<dim> mg_dof_handler; MatrixFree<double,LaplaceOperator<dim,double> > system_matrix; typedef MatrixFree<float,LaplaceOperator<dim,float> > MatrixFreeType; MGLevelObject<MatrixFreeType> mg_matrices; FullMatrix<float> coarse_matrix; Vector<double> solution; Vector<double> system_rhs; }; template <int dim> LaplaceProblem<dim>::LaplaceProblem (const unsigned int degree) : fe (degree), mg_dof_handler (triangulation) {}
This is the function of step-16 with relevant changes due to the MatrixFree class. What we need to do is to somehow create a local gradient matrix that does not contain any cell-related data (gradient on the reference cell). The way to get to this matrix is to create an FEValues object with gradient information on a cell that corresponds to the reference cell, which is a cube with side length 1. So we create a pseudo triangulation, initialize the FEValues to the only cell of that triangulation, and read off the gradients (which we put in a FullMatrix). That full matrix is then passed to the reinit function of the MatrixFree class used as a system matrix and, further down, as multigrid matrices on the individual levels. We need to implement Dirichlet boundary conditions here, which is done with the ConstraintMatrix function as shown, e.g., in step-22.
template <int dim> void LaplaceProblem<dim>::setup_system () { system_matrix.clear(); mg_matrices.clear(); mg_dof_handler.distribute_dofs (fe); std::cout << "Number of degrees of freedom: " << mg_dof_handler.n_dofs() << std::endl; const unsigned int nlevels = triangulation.n_levels(); mg_matrices.resize(0, nlevels-1); QGauss<dim> quadrature_formula(fe.degree+1); FEValues<dim> fe_values_reference (fe, quadrature_formula, update_gradients); Triangulation<dim> reference_cell; GridGenerator::hyper_cube (reference_cell, 0, 1); fe_values_reference.reinit (reference_cell.begin()); FullMatrix<double> ref_cell_gradients (fe.dofs_per_cell, quadrature_formula.size()*dim); for (unsigned int i=0; i<fe.dofs_per_cell; ++i) { for (unsigned int j=0; j<quadrature_formula.size(); ++j) { for (unsigned int d=0; d<dim; ++d) ref_cell_gradients(i,j*dim+d) = fe_values_reference.shape_grad(i,j)[d]; } } system_matrix.reinit (mg_dof_handler.n_dofs(), triangulation.n_active_cells(), ref_cell_gradients, quadrature_formula.size()); VectorTools::interpolate_boundary_values (mg_dof_handler, 0, ZeroFunction<dim>(), system_matrix.get_constraints()); system_matrix.get_constraints().close(); std::cout.precision(4); std::cout << "System matrix memory consumption: " << system_matrix.memory_consumption()/double(1<<20) << " MiB." << std::endl; solution.reinit (mg_dof_handler.n_dofs()); system_rhs.reinit (mg_dof_handler.n_dofs());
Next, initialize the matrices for the multigrid method on all the levels. Unfortunately, the function MGTools::make_boundary_list cannot write Dirichlet boundary conditions into a ConstraintMatrix object directly, so we first have to make the boundary list and then manually fill the boundary conditions using the command ConstraintMatrix::add_line. Once this is done, we close the ConstraintMatrix so it can be used for matrix-vector products.
typename FunctionMap<dim>::type dirichlet_boundary; ZeroFunction<dim> homogeneous_dirichlet_bc (1); dirichlet_boundary[0] = &homogeneous_dirichlet_bc; std::vector<std::set<unsigned int> > boundary_indices(triangulation.n_levels()); MGTools::make_boundary_list (mg_dof_handler, dirichlet_boundary, boundary_indices); for (unsigned int level=0;level<nlevels;++level) { mg_matrices[level].reinit(mg_dof_handler.n_dofs(level), triangulation.n_cells(level), ref_cell_gradients, quadrature_formula.size()); std::set<unsigned int>::iterator bc_it = boundary_indices[level].begin(); for ( ; bc_it != boundary_indices[level].end(); ++bc_it) mg_matrices[level].get_constraints().add_line(*bc_it); mg_matrices[level].get_constraints().close(); } coarse_matrix.reinit (mg_dof_handler.n_dofs(0), mg_dof_handler.n_dofs(0)); }
The assemble function is significantly reduced compared to step-16. All we need to do is to assemble the right hand side and to calculate the cell-dependent part of the Laplace operator. The first task is standard. The second is also not too hard given the discussion in the introduction: We need to take the inverse of the Jacobian of the transformation from unit to real cell, multiply it with its transpose and multiply the resulting rank-2 tensor with the quadrature weights and the coefficient values at the quadrature points. To make this work, we add the update flag update_inverse_jacobians to the FEValues constructor, and query the inverse of the Jacobian in a loop over the quadrature points (note that the Jacobian is not related to any kind of degrees of freedom directly). In the end, we condense the constraints from Dirichlet boundary conditions away from the right hand side.
template <int dim> void LaplaceProblem<dim>::assemble_system () { QGauss<dim> quadrature_formula(fe.degree+1); MappingQ<dim> mapping (fe.degree); FEValues<dim> fe_values (mapping, fe, quadrature_formula, update_values | update_inverse_jacobians | update_quadrature_points | update_JxW_values); const unsigned int dofs_per_cell = fe.dofs_per_cell; const unsigned int n_q_points = quadrature_formula.size(); std::vector<unsigned int> local_dof_indices (dofs_per_cell); const Coefficient<dim> coefficient; std::vector<double> coefficient_values (n_q_points); unsigned int cell_no = 0; typename DoFHandler<dim>::active_cell_iterator cell = mg_dof_handler.begin_active(), endc = mg_dof_handler.end(); for (; cell!=endc; ++cell, ++cell_no) { cell->get_dof_indices (local_dof_indices); fe_values.reinit (cell); coefficient.value_list (fe_values.get_quadrature_points(), coefficient_values); for (unsigned int i=0; i<dofs_per_cell; ++i) { double rhs_val = 0; for (unsigned int q=0; q<n_q_points; ++q) rhs_val += (fe_values.shape_value(i,q) * 1.0 * fe_values.JxW(q)); system_rhs(local_dof_indices[i]) += rhs_val; } system_matrix.set_local_dof_indices (cell_no, local_dof_indices); for (unsigned int q=0; q<n_q_points; ++q) system_matrix.set_derivative_data (cell_no, q, (transpose (fe_values.inverse_jacobian(q)) * fe_values.inverse_jacobian(q)) * fe_values.JxW(q) * coefficient_values[q]); } system_matrix.get_constraints().condense(system_rhs); }
Here is another assemble function. The integration core is the same as above. Only the loop goes over all existing cells now and the results must be entered into the correct matrix.
Since we only do multilevel preconditioning, no right-hand side is assembled here. Compared to step-16, there is one new thing here: we manually calculate the matrix on the coarsest level. In step-16, we could simply copy the entries from the respective sparse matrix, but this is obviously not possible here. We could have integrated this into the MatrixFree class as well, but it is simple enough, so calculate it here instead.
template <int dim> void LaplaceProblem<dim>::assemble_multigrid () { coarse_matrix = 0; QGauss<dim> quadrature_formula(fe.degree+1); MappingQ<dim> mapping (fe.degree); FEValues<dim> fe_values (mapping, fe, quadrature_formula, update_gradients | update_inverse_jacobians | update_quadrature_points | update_JxW_values); const unsigned int dofs_per_cell = fe.dofs_per_cell; const unsigned int n_q_points = quadrature_formula.size(); std::vector<unsigned int> local_dof_indices (dofs_per_cell); const Coefficient<dim> coefficient; std::vector<double> coefficient_values (n_q_points); std::vector<unsigned int> cell_no(triangulation.n_levels()); typename MGDoFHandler<dim>::cell_iterator cell = mg_dof_handler.begin(), endc = mg_dof_handler.end(); for (; cell!=endc; ++cell) { const unsigned int level = cell->level(); cell->get_mg_dof_indices (local_dof_indices); fe_values.reinit (cell); coefficient.value_list (fe_values.get_quadrature_points(), coefficient_values); mg_matrices[level].set_local_dof_indices (cell_no[level], local_dof_indices); for (unsigned int q=0; q<n_q_points; ++q) mg_matrices[level].set_derivative_data (cell_no[level], q, (transpose (fe_values.inverse_jacobian(q)) * fe_values.inverse_jacobian(q)) * fe_values.JxW(q) * coefficient_values[q]); ++cell_no[level]; if (level == 0) { for (unsigned int i=0; i<dofs_per_cell; ++i) for (unsigned int j=0; j<dofs_per_cell; ++j) { double add_value = 0; for (unsigned int q=0; q<n_q_points; ++q) add_value += (fe_values.shape_grad(i,q) * fe_values.shape_grad(j,q) * coefficient_values[q] * fe_values.JxW(q)); coarse_matrix(local_dof_indices[i], local_dof_indices[j]) += add_value; } } }
In a final step, we need to condense the boundary conditions on the coarse matrix. There is no built-in function for doing this on a full matrix, so manually delete the rows and columns of the matrix that are constrained.
for (unsigned int i=0; i<coarse_matrix.m(); ++i) if (mg_matrices[0].get_constraints().is_constrained(i)) for (unsigned int j=0; j<coarse_matrix.n(); ++j) if (i!=j) { coarse_matrix(i,j) = 0; coarse_matrix(j,i) = 0; } }
The solution process again looks like step-16. We now use a Chebyshev smoother instead of SOR (SOR would be very difficult to implement because we do not have the matrix elements available explicitly, and it is difficult to make it work efficiently in parallel). The multigrid classes provide a simple interface for using the Chebyshev smoother which is defined in a preconditioner class: MGSmootherPrecondition.
template <int dim> void LaplaceProblem<dim>::solve () { GrowingVectorMemory<> vector_memory; MGTransferPrebuilt<Vector<double> > mg_transfer; mg_transfer.build_matrices(mg_dof_handler); MGCoarseGridHouseholder<float, Vector<double> > mg_coarse; mg_coarse.initialize(coarse_matrix); typedef PreconditionChebyshev<MatrixFreeType,Vector<double> > SMOOTHER; MGSmootherPrecondition<MatrixFreeType, SMOOTHER, Vector<double> > mg_smoother(vector_memory);
Then, we initialize the smoother with our level matrices and the required, additional data for the Chebyshev smoother. In particular, we use a higher polynomial degree for higher order elements, since smoothing gets more difficult for these. Smooth out a range of
. In order to compute the maximum eigenvalue of the corresponding matrix, the Chebyshev initializations performs a few steps of a CG algorithm. Since all we need is a rough estimate, we choose some eight iterations (more if the finite element polynomial degree is larger, less if it is smaller than quadratic).
typename SMOOTHER::AdditionalData smoother_data; smoother_data.smoothing_range = 10.; smoother_data.degree = fe.degree; smoother_data.eig_cg_n_iterations = 4+2*fe.degree; mg_smoother.initialize(mg_matrices, smoother_data); MGMatrix<MatrixFreeType, Vector<double> > mg_matrix(&mg_matrices); Multigrid<Vector<double> > mg(mg_dof_handler, mg_matrix, mg_coarse, mg_transfer, mg_smoother, mg_smoother); PreconditionMG<dim, Vector<double>, MGTransferPrebuilt<Vector<double> > > preconditioner(mg_dof_handler, mg, mg_transfer);
Finally, write out the memory consumption of the Multigrid object (or rather, of its most significant components, since there is no built-in function for the total multigrid object), then create the solver object and solve the system. This is very easy, and we didn't even see any difference in the solve process compared to step-16. The magic is all hidden behind the implementation of the MatrixFree::vmult operation.
const unsigned int multigrid_memory = (mg_matrices.memory_consumption() + mg_transfer.memory_consumption() + coarse_matrix.memory_consumption()); std::cout << "Multigrid objects memory consumption: " << multigrid_memory/double(1<<20) << " MiB." << std::endl; SolverControl solver_control (1000, 1e-12); SolverCG<> cg (solver_control); cg.solve (system_matrix, solution, system_rhs, preconditioner); std::cout << "Convergence in " << solver_control.last_step() << " CG iterations." << std::endl; }
Here is the data output, which is a simplified version of step-5. We use the standard VTK output for each grid produced in the refinement process.
template <int dim> void LaplaceProblem<dim>::output_results (const unsigned int cycle) const { DataOut<dim> data_out; data_out.attach_dof_handler (mg_dof_handler); data_out.add_data_vector (solution, "solution"); data_out.build_patches (); std::ostringstream filename; filename << "solution-" << cycle << ".vtk"; std::ofstream output (filename.str().c_str()); data_out.write_vtk (output); }
The function that runs the program is very similar to the one in step-16. We make less refinement steps in 3D compared to 2D, but that's it.
template <int dim> void LaplaceProblem<dim>::run () { for (unsigned int cycle=0; cycle<8-dim; ++cycle) { std::cout << "Cycle " << cycle << std::endl; if (cycle == 0) { GridGenerator::hyper_ball(triangulation); static const HyperBallBoundary<dim> boundary; triangulation.set_boundary (0, boundary); triangulation.refine_global (3-dim); } triangulation.refine_global (1); setup_system (); assemble_system (); assemble_multigrid (); solve (); output_results (cycle); std::cout << std::endl; }; }
main functionThis is as in all other programs:
int main () { try { deallog.depth_console (0); LaplaceProblem<2> laplace_problem (2); laplace_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; }
Since this example solves the same problem as step-5 (except for a different coefficient), we refer to the graphical output there. Here, we evaluate some aspects of the multigrid solver.
When we run this program in 2D for quadratic (
) elements, we get the following output:
Cycle 0 Number of degrees of freedom: 337 System matrix memory consumption: 0.02573 MiB. Multigrid objects memory consumption: 0.05083 MiB. Convergence in 10 CG iterations. Cycle 1 Number of degrees of freedom: 1313 System matrix memory consumption: 0.09257 MiB. Multigrid objects memory consumption: 0.1794 MiB. Convergence in 10 CG iterations. Cycle 2 Number of degrees of freedom: 5185 System matrix memory consumption: 0.3553 MiB. Multigrid objects memory consumption: 0.6779 MiB. Convergence in 10 CG iterations. Cycle 3 Number of degrees of freedom: 20609 System matrix memory consumption: 1.397 MiB. Multigrid objects memory consumption: 2.645 MiB. Convergence in 10 CG iterations. Cycle 4 Number of degrees of freedom: 82177 System matrix memory consumption: 5.546 MiB. Multigrid objects memory consumption: 10.46 MiB. Convergence in 10 CG iterations. Cycle 5 Number of degrees of freedom: 328193 System matrix memory consumption: 22.11 MiB. Multigrid objects memory consumption: 41.65 MiB. Convergence in 10 CG iterations.
As in step-16, we see that the number of CG iterations remains constant with increasing number of degrees of freedom. We can also see that the various objects we have to store for the multigrid method on the individual levels of our mesh together make up about twice as much as the matrix on the finest level.
Not much changes if we run the program in three spatial dimensions, with the exception that the multilevel objects now take up comparatively less space (because in 3d, each level has only one eighth the number of cells of the next finer one, whereas in 2d this factor if one quarter):
Cycle 0 Number of degrees of freedom: 517 System matrix memory consumption: 0.1001 MiB. Multigrid objects memory consumption: 0.1463 MiB. Convergence in 9 CG iterations. Cycle 1 Number of degrees of freedom: 3817 System matrix memory consumption: 0.6613 MiB. Multigrid objects memory consumption: 0.8896 MiB. Convergence in 10 CG iterations. Cycle 2 Number of degrees of freedom: 29521 System matrix memory consumption: 5.1 MiB. Multigrid objects memory consumption: 6.653 MiB. Convergence in 10 CG iterations. Cycle 3 Number of degrees of freedom: 232609 System matrix memory consumption: 40.4 MiB. Multigrid objects memory consumption: 52.24 MiB. Convergence in 11 CG iterations. Cycle 4 Number of degrees of freedom: 1847617 System matrix memory consumption: 322 MiB. Multigrid objects memory consumption: 415.1 MiB. Convergence in 11 CG iterations.
In order to understand the capabilities of this class, we compare the memory consumption and execution (wallclock) time for assembly and 50 matrix-vector products (MV) on a 3D problem with one million unknowns the classical sparse matrix implementation (SpM) and the MatrixFree implementation shown here (M-F). Both matrices are based on double numbers. The program is run on a 2.8 GHz Opteron processor with the ACML BLAS. We present results running on one core core and four cores, respectively. Moreover, we measure the time it takes to construct the individual matrices and filling them with data (setup and assemble functions). The sparse matrix is initialized using a CompressedSimpleSparsityPattern for calling the DoFTools::make_sparsity_pattern function, and then copied to a SparsityPattern object. The boundary nodes are eliminated using the ConstraintMatrix class, so that only elements that are actually nonzero are stored in the matrix.
| Memory consumption | Time assembly | Time 50 MV, 1 CPU | Time 50 MV, 4 CPUs | |||||
|---|---|---|---|---|---|---|---|---|
| element order | SpM | M-F | SpM | M-F | SpM | M-F | SpM | M-F |
| 1 | 299 MiB | 394 MiB | 8.09 s | 3.43 s | 5.50 s | 22.4 s | 4.30 s | 11.0 s |
| 2 | 698 MiB | 177 MiB | 12.43 s | 1.32 s | 12.0 s | 18.6 s | 9.10 s | 6.31 s |
| 3 | 1295 MiB | 124 MiB | 41.1 s | 1.31 s | 21.2 s | 23.7 s | 16.0 s | 7.43 s |
| 4 | 2282 MiB | 107 MiB | 117 s | 1.97 s | 40.8 s | 36.3 s | 19.7 s | 10.9 s |
| 5 | 3597 MiB | 96.4 MiB | 510 s | 5.52 s | 75.7 s | 53.9 s | 29.3 s | 15.7 s |
| 6 | 5679 MiB | 96.3 MiB | 2389 s | 26.1 s | 135 s | 79.1 s | 45.8 s | 24.3 s |
There are a few interesting things with the numbers in this table.
Firstly, we see the disappointing fact that for linear elements the MatrixFree class does actually consume more memory than a SparseMatrix with its SparsityPattern, despite the efforts made in this program. As mentioned earlier, this is mostly because the Transformation data is stored for every quadrature point. For each quadrature point, the transformation consists of six doubles, and there are about eight times as many quadrature points as there are degrees of freedom. In first approximation, this means that the matrix consumes 384 (= sizeof(double) * 6 * 8) bytes for each degree of freedom. On the other hand, the sparse matrix has a bandwidth of 27 or less, so each dof gives rise to at most 324 (= 27 * 12) bytes. A more clever implementation would try to compress the Jacobian transformation data, by exploiting similarities between the mappings within the cells, as well as from one cell to the next. This could dramatically reduce the memory requirements, and hence, increase the speed for lower-order implementations.
Secondly, we observe that the memory requirements for a SparseMatrix grow quickly as the order of the elements increases. This is because there are increasingly many entries in each row, which exist due to more degrees of freedom that couple to each other. The matrix-free implementation does not suffer from this drawback. Here, the memory consumption decreases instead, since there are less DoFs that are shared among elements, which decreases the relative amount of quadrature points. Regarding the execution speed, we see that the matrix-free variant gets more competitive with higher order, and it does scale better when run on multiple processors (3.5 speedup with four processors compared to the serial case, compared to 2-2.5 speedup for the SparseMatrix). The advantage in parallel scaling was expected, because the matrix-free variant is less memory-bound for higher order implementations, so that the additional computing power from many cores can better be exploited.
A third thing, which is unrelated to this tutorial program, is the fact that standard matrix assembly gets really slow for high order elements. The numbers shown here are based on the usual routines that many other tutorial programs make use of. A closer analysis of this shows that the cell data does not fit into cache anymore. One could circumvent this problem by writing the assembly as a matrix-matrix product, and using (cache-aware) BLAS implementations.
For completeness, here comes a similar table for a 2D problem with 5.7 million unknowns. Since the excess in work for the matrix-free implementation is less compared to 3D, the implementation is more competitive for lower-order elements.
| Memory consumption | Time assembly | Time 50 MV, 4 CPUs | ||||
|---|---|---|---|---|---|---|
| element order | SpM | M-F | SpM | M-F | SpM | M-F |
| 1 | 659 MiB | 661 MiB | 18.8 s | 6.45 s | 11.0 s | 28.8 s |
| 2 | 1119 MiB | 391 MiB | 15.6 s | 2.46 s | 17.1 s | 16.2 s |
| 3 | 1711 MiB | 318 MiB | 17.4 s | 1.82 s | 23.1 s | 13.7 s |
| 4 | 2434 MiB | 285 MiB | 24.2 s | 1.34 s | 31.1 s | 14.6 s |
| 5 | 3289 MiB | 266 MiB | 35.9 s | 1.26 s | 29.6 s | 16.7 s |
| 6 | 4274 MiB | 254 MiB | 58.0 s | 1.12 s | 35.9 s | 19.4 s |
(If you are looking at a locally installed deal.II version, then the program can be found at /u /bangerth /tmp /homepage /deal.II /examples /step-37 /step-37.cc . Otherwise, this is only the path on some remote server.)
/ * Subversion Id: @ref step_37 "step-37".cc 20922 2010-03-29 18:15:05Z kronbichler * / / * Author: Katharina Kormann, Martin Kronbichler, Uppsala University, 2009 * / / * Subversion Id: @ref step_37 "step-37".cc 20922 2010-03-29 18:15:05Z kronbichler * / / * * / / * Copyright (C) 2009 by the deal.II authors * / / * * / / * This file is subject to QPL and may not be distributed * / / * without copyright and license information. Please refer * / / * to the file deal.II/doc/license.html for the text and * / / * further information on this license. * / #include <base/quadrature_lib.h> #include <base/function.h> #include <base/logstream.h> #include <base/work_stream.h> #include <lac/vector.h> #include <lac/full_matrix.h> #include <lac/solver_cg.h> #include <lac/precondition.h> #include <fe/fe_q.h> #include <fe/fe_values.h> #include <grid/tria.h> #include <grid/tria_accessor.h> #include <grid/tria_iterator.h> #include <grid/tria_boundary_lib.h> #include <grid/grid_generator.h> #include <multigrid/multigrid.h> #include <multigrid/mg_dof_handler.h> #include <multigrid/mg_dof_accessor.h> #include <multigrid/mg_transfer.h> #include <multigrid/mg_tools.h> #include <multigrid/mg_coarse.h> #include <multigrid/mg_smoother.h> #include <multigrid/mg_matrix.h> #include <numerics/data_out.h> #include <numerics/vectors.h> #include <fstream> #include <sstream> using namespace dealii;
template <int dim> class Coefficient : public Function<dim> { public: Coefficient () : Function<dim>() {} virtual double value (const Point<dim> &p, const unsigned int component = 0) const; virtual void value_list (const std::vector<Point<dim> > &points, std::vector<double> &values, const unsigned int component = 0) const; }; template <int dim> double Coefficient<dim>::value (const Point<dim> &p, const unsigned int / *component* /) const { return 1./(0.1+p.square()); } template <int dim> void Coefficient<dim>::value_list (const std::vector<Point<dim> > &points, std::vector<double> &values, const unsigned int component) const { Assert (values.size() == points.size(), ExcDimensionMismatch (values.size(), points.size())); Assert (component == 0, ExcIndexRange (component, 0, 1)); const unsigned int n_points = points.size(); for (unsigned int i=0; i<n_points; ++i) values[i] = 1./(0.1+points[i].square()); }
namespace WorkStreamData { template <typename number> struct ScratchData { ScratchData (); ScratchData (const ScratchData &scratch); FullMatrix<number> solutions; }; template<typename number> ScratchData<number>::ScratchData () : solutions () {} template<typename number> ScratchData<number>::ScratchData (const ScratchData &) : solutions () {} template <typename number> struct CopyData : public ScratchData<number> { CopyData (); CopyData (const CopyData &scratch); unsigned int first_cell; unsigned int n_dofs; }; template <typename number> CopyData<number>::CopyData () : ScratchData<number> () {} template <typename number> CopyData<number>::CopyData (const CopyData &) : ScratchData<number> () {} } template <typename number, class Transformation> class MatrixFree : public Subscriptor { public: MatrixFree (); void reinit (const unsigned int n_dofs, const unsigned int n_cells, const FullMatrix<double> &cell_matrix, const unsigned int n_points_per_cell); void clear(); unsigned int m () const; unsigned int n () const; ConstraintMatrix & get_constraints (); void set_local_dof_indices (const unsigned int cell_no, const std::vector<unsigned int> &local_dof_indices); void set_derivative_data (const unsigned int cell_no, const unsigned int quad_point, const Transformation &trans_in); template <typename number2> void vmult (Vector<number2> &dst, const Vector<number2> &src) const; template <typename number2> void Tvmult (Vector<number2> &dst, const Vector<number2> &src) const; template <typename number2> void vmult_add (Vector<number2> &dst, const Vector<number2> &src) const; template <typename number2> void Tvmult_add (Vector<number2> &dst, const Vector<number2> &src) const; number el (const unsigned int row, const unsigned int col) const; void calculate_diagonal () const; std::size_t memory_consumption () const; private: typedef std::vector<std::pair<unsigned int,unsigned int> >::const_iterator CellChunkIterator; template <typename number2> void local_vmult (CellChunkIterator cell_range, WorkStreamData::ScratchData<number> &scratch, WorkStreamData::CopyData<number> ©, const Vector<number2> &src) const; template <typename number2> void copy_local_to_global (const WorkStreamData::CopyData<number> ©, Vector<number2> &dst) const; FullMatrix<number> B_ref_cell; Table<2,unsigned int> indices_local_to_global; Table<2,Transformation> derivatives; ConstraintMatrix constraints; mutable Vector<number> diagonal_values; mutable bool diagonal_is_calculated; struct MatrixSizes { unsigned int n_dofs, n_cells; unsigned int m, n; unsigned int n_points, n_comp; std::vector<std::pair<unsigned int,unsigned int> > chunks; } matrix_sizes; }; template <typename number, class Transformation> MatrixFree<number,Transformation>::MatrixFree () : Subscriptor() {} template <typename number, class Transformation> unsigned int MatrixFree<number,Transformation>::m () const { return matrix_sizes.n_dofs; } template <typename number, class Transformation> unsigned int MatrixFree<number,Transformation>::n () const { return matrix_sizes.n_dofs; } template <typename number, class Transformation> ConstraintMatrix & MatrixFree<number,Transformation>::get_constraints () { return constraints; } template <typename number, class Transformation> void MatrixFree<number,Transformation>:: set_local_dof_indices (const unsigned int cell_no, const std::vector<unsigned int> &local_dof_indices) { Assert (local_dof_indices.size() == matrix_sizes.m, ExcDimensionMismatch(local_dof_indices.size(), matrix_sizes.m)); for (unsigned int i=0; i<matrix_sizes.m; ++i) { Assert (local_dof_indices[i] < matrix_sizes.n_dofs, ExcInternalError()); indices_local_to_global(cell_no,i) = local_dof_indices[i]; } diagonal_is_calculated = false; } template <typename number, class Transformation> void MatrixFree<number,Transformation>:: set_derivative_data (const unsigned int cell_no, const unsigned int quad_point, const Transformation &trans_in) { Assert (quad_point < matrix_sizes.n_points, ExcInternalError()); derivatives(cell_no,quad_point) = trans_in; diagonal_is_calculated = false; } template <typename number, class Transformation> template <typename number2> void MatrixFree<number,Transformation>:: local_vmult (CellChunkIterator cell_range, WorkStreamData::ScratchData<number> &scratch, WorkStreamData::CopyData<number> ©, const Vector<number2> &src) const { const unsigned int chunk_size = cell_range->second - cell_range->first; scratch.solutions.reinit (chunk_size, matrix_sizes.n, true); copy.solutions.reinit (chunk_size, matrix_sizes.m, true); copy.first_cell = cell_range->first; copy.n_dofs = chunk_size*matrix_sizes.m; constraints.get_dof_values(src, &indices_local_to_global(copy.first_cell,0), ©.solutions(0,0), ©.solutions(0,0)+copy.n_dofs); copy.solutions.mmult (scratch.solutions, B_ref_cell); for (unsigned int i=0, k = copy.first_cell; i<chunk_size; ++i, ++k) for (unsigned int j=0; j<matrix_sizes.n_points; ++j) derivatives(k,j).transform(&scratch.solutions(i, j*matrix_sizes.n_comp)); scratch.solutions.mTmult (copy.solutions, B_ref_cell); } template <typename number, class Transformation> template <typename number2> void MatrixFree<number,Transformation>:: copy_local_to_global (const WorkStreamData::CopyData<number> ©, Vector<number2> &dst) const { constraints.distribute_local_to_global (©.solutions(0,0), ©.solutions(0,0)+copy.n_dofs, &indices_local_to_global(copy.first_cell,0), dst); } template <typename number, class Transformation> template <typename number2> void MatrixFree<number,Transformation>::vmult (Vector<number2> &dst, const Vector<number2> &src) const { dst = 0; vmult_add (dst, src); } template <typename number, class Transformation> template <typename number2> void MatrixFree<number,Transformation>::Tvmult (Vector<number2> &dst, const Vector<number2> &src) const { dst = 0; Tvmult_add (dst,src); } template <typename number, class Transformation> template <typename number2> void MatrixFree<number,Transformation>::Tvmult_add (Vector<number2> &dst, const Vector<number2> &src) const { vmult_add (dst,src); } template <typename number, class Transformation> template <typename number2> void MatrixFree<number,Transformation>::vmult_add (Vector<number2> &dst, const Vector<number2> &src) const { Assert (src.size() == n(), ExcDimensionMismatch(src.size(), n())); Assert (dst.size() == m(), ExcDimensionMismatch(dst.size(), m())); WorkStream::run (matrix_sizes.chunks.begin(), matrix_sizes.chunks.end(), std_cxx1x::bind(&MatrixFree<number,Transformation>:: template local_vmult<number2>, this, _1, _2, _3, boost::cref(src)), std_cxx1x::bind(&MatrixFree<number,Transformation>:: template copy_local_to_global<number2>, this, _1, boost::ref(dst)), WorkStreamData::ScratchData<number>(), WorkStreamData::CopyData<number>(), 2*multithread_info.n_default_threads,1); for (unsigned int i=0; i<matrix_sizes.n_dofs; ++i) if (constraints.is_constrained(i) == true) dst(i) += 1.0 * src(i); } template <typename number, class Transformation> void MatrixFree<number,Transformation>:: reinit (const unsigned int n_dofs_in, const unsigned int n_cells_in, const FullMatrix<double> &B_ref_cell_in, const unsigned int n_points_per_cell) { B_ref_cell = B_ref_cell_in; derivatives.reinit (n_cells_in, n_points_per_cell); indices_local_to_global.reinit (n_cells_in, B_ref_cell.m()); diagonal_is_calculated = false; matrix_sizes.n_dofs = n_dofs_in; matrix_sizes.n_cells = n_cells_in; matrix_sizes.m = B_ref_cell.m(); matrix_sizes.n = B_ref_cell.n(); matrix_sizes.n_points = n_points_per_cell; matrix_sizes.n_comp = B_ref_cell.n()/matrix_sizes.n_points; Assert(matrix_sizes.n_comp * n_points_per_cell == B_ref_cell.n(), ExcInternalError()); const unsigned int divisor = std::max(60U, matrix_sizes.m); const unsigned int n_chunks = std::max (matrix_sizes.n_cells/divisor + 1, 2*multithread_info.n_default_threads); const unsigned int chunk_size = (matrix_sizes.n_cells/n_chunks + (matrix_sizes.n_cells%n_chunks>0)); std::pair<unsigned int, unsigned int> chunk; for (unsigned int i=0; i<n_chunks; ++i) { chunk.first = i*chunk_size; if ((i+1)*chunk_size > matrix_sizes.n_cells) chunk.second = matrix_sizes.n_cells; else chunk.second = (i+1)*chunk_size; if (chunk.second > chunk.first) matrix_sizes.chunks.push_back(chunk); else break; } } template <typename number, class Transformation> void MatrixFree<number,Transformation>::clear () { B_ref_cell.reinit(0,0); derivatives.reinit (0,0); indices_local_to_global.reinit(0,0); constraints.clear(); diagonal_values.reinit (0); diagonal_is_calculated = false; matrix_sizes.n_dofs = 0; matrix_sizes.n_cells = 0; matrix_sizes.chunks.clear(); } template <typename number, class Transformation> number MatrixFree<number,Transformation>::el (const unsigned int row, const unsigned int col) const { Assert (row == col, ExcNotImplemented()); if (diagonal_is_calculated == false) calculate_diagonal(); return diagonal_values(row); } template <typename number, class Transformation> void MatrixFree<number,Transformation>::calculate_diagonal() const { diagonal_values.reinit (matrix_sizes.n_dofs); std::vector<number> calculation (matrix_sizes.n); for (unsigned int cell=0; cell<matrix_sizes.n_cells; ++cell) for (unsigned int dof=0; dof<matrix_sizes.m; ++dof) { memcpy (&calculation[0],&B_ref_cell(dof,0), matrix_sizes.n*sizeof(number)); for (unsigned int q=0; q<matrix_sizes.n_points; ++q) derivatives(cell,q).transform(&calculation[q*matrix_sizes.n_comp]); double diag_value = 0; for (unsigned int q=0; q<matrix_sizes.n; ++q) diag_value += calculation[q] * B_ref_cell(dof,q); diagonal_values(indices_local_to_global(cell,dof)) += diag_value; } constraints.condense (diagonal_values); for (unsigned int i=0; i<matrix_sizes.n_dofs; ++i) if (constraints.is_constrained(i) == true) diagonal_values(i) = 1.0; diagonal_is_calculated = true; } template <typename number, class Transformation> std::size_t MatrixFree<number,Transformation>::memory_consumption () const { std::size_t glob_size = derivatives.memory_consumption() + indices_local_to_global.memory_consumption() + constraints.memory_consumption() + B_ref_cell.memory_consumption() + diagonal_values.memory_consumption() + matrix_sizes.chunks.size()*2*sizeof(unsigned int) + sizeof(*this); return glob_size; }
template <int dim,typename number> class LaplaceOperator { public: LaplaceOperator (); LaplaceOperator (const Tensor<2,dim> &tensor); void transform (number * result) const; LaplaceOperator<dim,number>& operator = (const Tensor<2,dim> &tensor); private: number transformation[dim*(dim+1)/2]; }; template<int dim,typename number> LaplaceOperator<dim,number>::LaplaceOperator() {} template<int dim,typename number> LaplaceOperator<dim,number>::LaplaceOperator(const Tensor<2,dim> &tensor) { *this = tensor; } template <int dim, typename number> void LaplaceOperator<dim,number>::transform (number* result) const { if (dim == 2) { const number temp = result[0]; result[0] = transformation[0] * temp + transformation[1] * result[1]; result[1] = transformation[1] * temp + transformation[2] * result[1]; } else if (dim == 3) { const number temp1 = result[0]; const number temp2 = result[1]; result[0] = transformation[0] * temp1 + transformation[1] * temp2 + transformation[2] * result[2]; result[1] = transformation[1] * temp1 + transformation[3] * temp2 + transformation[4] * result[2]; result[2] = transformation[2] * temp1 + transformation[4] * temp2 + transformation[5] * result[2]; } else ExcNotImplemented(); } template <int dim, typename number> LaplaceOperator<dim,number>& LaplaceOperator<dim,number>::operator=(const Tensor<2,dim> &tensor) { if (dim == 2) { transformation[0] = tensor[0][0]; transformation[1] = tensor[0][1]; transformation[2] = tensor[1][1]; Assert (std::fabs(tensor[1][0]-tensor[0][1])<1e-15, ExcInternalError()); } else if (dim == 3) { transformation[0] = tensor[0][0]; transformation[1] = tensor[0][1]; transformation[2] = tensor[0][2]; transformation[3] = tensor[1][1]; transformation[4] = tensor[1][2]; transformation[5] = tensor[2][2]; Assert (std::fabs(tensor[1][0]-tensor[0][1])<1e-15, ExcInternalError()); Assert (std::fabs(tensor[2][0]-tensor[0][2])<1e-15, ExcInternalError()); Assert (std::fabs(tensor[2][1]-tensor[1][2])<1e-15, ExcInternalError()); } else ExcNotImplemented(); return *this; }
template <int dim> class LaplaceProblem { public: LaplaceProblem (const unsigned int degree); void run (); private: void setup_system (); void assemble_system (); void assemble_multigrid (); void solve (); void output_results (const unsigned int cycle) const; Triangulation<dim> triangulation; FE_Q<dim> fe; MGDoFHandler<dim> mg_dof_handler; MatrixFree<double,LaplaceOperator<dim,double> > system_matrix; typedef MatrixFree<float,LaplaceOperator<dim,float> > MatrixFreeType; MGLevelObject<MatrixFreeType> mg_matrices; FullMatrix<float> coarse_matrix; Vector<double> solution; Vector<double> system_rhs; }; template <int dim> LaplaceProblem<dim>::LaplaceProblem (const unsigned int degree) : fe (degree), mg_dof_handler (triangulation) {}
template <int dim> void LaplaceProblem<dim>::setup_system () { system_matrix.clear(); mg_matrices.clear(); mg_dof_handler.distribute_dofs (fe); std::cout << "Number of degrees of freedom: " << mg_dof_handler.n_dofs() << std::endl; const unsigned int nlevels = triangulation.n_levels(); mg_matrices.resize(0, nlevels-1); QGauss<dim> quadrature_formula(fe.degree+1); FEValues<dim> fe_values_reference (fe, quadrature_formula, update_gradients); Triangulation<dim> reference_cell; GridGenerator::hyper_cube (reference_cell, 0, 1); fe_values_reference.reinit (reference_cell.begin()); FullMatrix<double> ref_cell_gradients (fe.dofs_per_cell, quadrature_formula.size()*dim); for (unsigned int i=0; i<fe.dofs_per_cell; ++i) { for (unsigned int j=0; j<quadrature_formula.size(); ++j) { for (unsigned int d=0; d<dim; ++d) ref_cell_gradients(i,j*dim+d) = fe_values_reference.shape_grad(i,j)[d]; } } system_matrix.reinit (mg_dof_handler.n_dofs(), triangulation.n_active_cells(), ref_cell_gradients, quadrature_formula.size()); VectorTools::interpolate_boundary_values (mg_dof_handler, 0, ZeroFunction<dim>(), system_matrix.get_constraints()); system_matrix.get_constraints().close(); std::cout.precision(4); std::cout << "System matrix memory consumption: " << system_matrix.memory_consumption()/double(1<<20) << " MiB." << std::endl; solution.reinit (mg_dof_handler.n_dofs()); system_rhs.reinit (mg_dof_handler.n_dofs()); typename FunctionMap<dim>::type dirichlet_boundary; ZeroFunction<dim> homogeneous_dirichlet_bc (1); dirichlet_boundary[0] = &homogeneous_dirichlet_bc; std::vector<std::set<unsigned int> > boundary_indices(triangulation.n_levels()); MGTools::make_boundary_list (mg_dof_handler, dirichlet_boundary, boundary_indices); for (unsigned int level=0;level<nlevels;++level) { mg_matrices[level].reinit(mg_dof_handler.n_dofs(level), triangulation.n_cells(level), ref_cell_gradients, quadrature_formula.size()); std::set<unsigned int>::iterator bc_it = boundary_indices[level].begin(); for ( ; bc_it != boundary_indices[level].end(); ++bc_it) mg_matrices[level].get_constraints().add_line(*bc_it); mg_matrices[level].get_constraints().close(); } coarse_matrix.reinit (mg_dof_handler.n_dofs(0), mg_dof_handler.n_dofs(0)); }
template <int dim> void LaplaceProblem<dim>::assemble_system () { QGauss<dim> quadrature_formula(fe.degree+1); MappingQ<dim> mapping (fe.degree); FEValues<dim> fe_values (mapping, fe, quadrature_formula, update_values | update_inverse_jacobians | update_quadrature_points | update_JxW_values); const unsigned int dofs_per_cell = fe.dofs_per_cell; const unsigned int n_q_points = quadrature_formula.size(); std::vector<unsigned int> local_dof_indices (dofs_per_cell); const Coefficient<dim> coefficient; std::vector<double> coefficient_values (n_q_points); unsigned int cell_no = 0; typename DoFHandler<dim>::active_cell_iterator cell = mg_dof_handler.begin_active(), endc = mg_dof_handler.end(); for (; cell!=endc; ++cell, ++cell_no) { cell->get_dof_indices (local_dof_indices); fe_values.reinit (cell); coefficient.value_list (fe_values.get_quadrature_points(), coefficient_values); for (unsigned int i=0; i<dofs_per_cell; ++i) { double rhs_val = 0; for (unsigned int q=0; q<n_q_points; ++q) rhs_val += (fe_values.shape_value(i,q) * 1.0 * fe_values.JxW(q)); system_rhs(local_dof_indices[i]) += rhs_val; } system_matrix.set_local_dof_indices (cell_no, local_dof_indices); for (unsigned int q=0; q<n_q_points; ++q) system_matrix.set_derivative_data (cell_no, q, (transpose (fe_values.inverse_jacobian(q)) * fe_values.inverse_jacobian(q)) * fe_values.JxW(q) * coefficient_values[q]); } system_matrix.get_constraints().condense(system_rhs); }
template <int dim> void LaplaceProblem<dim>::assemble_multigrid () { coarse_matrix = 0; QGauss<dim> quadrature_formula(fe.degree+1); MappingQ<dim> mapping (fe.degree); FEValues<dim> fe_values (mapping, fe, quadrature_formula, update_gradients | update_inverse_jacobians | update_quadrature_points | update_JxW_values); const unsigned int dofs_per_cell = fe.dofs_per_cell; const unsigned int n_q_points = quadrature_formula.size(); std::vector<unsigned int> local_dof_indices (dofs_per_cell); const Coefficient<dim> coefficient; std::vector<double> coefficient_values (n_q_points); std::vector<unsigned int> cell_no(triangulation.n_levels()); typename MGDoFHandler<dim>::cell_iterator cell = mg_dof_handler.begin(), endc = mg_dof_handler.end(); for (; cell!=endc; ++cell) { const unsigned int level = cell->level(); cell->get_mg_dof_indices (local_dof_indices); fe_values.reinit (cell); coefficient.value_list (fe_values.get_quadrature_points(), coefficient_values); mg_matrices[level].set_local_dof_indices (cell_no[level], local_dof_indices); for (unsigned int q=0; q<n_q_points; ++q) mg_matrices[level].set_derivative_data (cell_no[level], q, (transpose (fe_values.inverse_jacobian(q)) * fe_values.inverse_jacobian(q)) * fe_values.JxW(q) * coefficient_values[q]); ++cell_no[level]; if (level == 0) { for (unsigned int i=0; i<dofs_per_cell; ++i) for (unsigned int j=0; j<dofs_per_cell; ++j) { double add_value = 0; for (unsigned int q=0; q<n_q_points; ++q) add_value += (fe_values.shape_grad(i,q) * fe_values.shape_grad(j,q) * coefficient_values[q] * fe_values.JxW(q)); coarse_matrix(local_dof_indices[i], local_dof_indices[j]) += add_value; } } } for (unsigned int i=0; i<coarse_matrix.m(); ++i) if (mg_matrices[0].get_constraints().is_constrained(i)) for (unsigned int j=0; j<coarse_matrix.n(); ++j) if (i!=j) { coarse_matrix(i,j) = 0; coarse_matrix(j,i) = 0; } }
template <int dim> void LaplaceProblem<dim>::solve () { GrowingVectorMemory<> vector_memory; MGTransferPrebuilt<Vector<double> > mg_transfer; mg_transfer.build_matrices(mg_dof_handler); MGCoarseGridHouseholder<float, Vector<double> > mg_coarse; mg_coarse.initialize(coarse_matrix); typedef PreconditionChebyshev<MatrixFreeType,Vector<double> > SMOOTHER; MGSmootherPrecondition<MatrixFreeType, SMOOTHER, Vector<double> > mg_smoother(vector_memory); typename SMOOTHER::AdditionalData smoother_data; smoother_data.smoothing_range = 10.; smoother_data.degree = fe.degree; smoother_data.eig_cg_n_iterations = 4+2*fe.degree; mg_smoother.initialize(mg_matrices, smoother_data); MGMatrix<MatrixFreeType, Vector<double> > mg_matrix(&mg_matrices); Multigrid<Vector<double> > mg(mg_dof_handler, mg_matrix, mg_coarse, mg_transfer, mg_smoother, mg_smoother); PreconditionMG<dim, Vector<double>, MGTransferPrebuilt<Vector<double> > > preconditioner(mg_dof_handler, mg, mg_transfer); const unsigned int multigrid_memory = (mg_matrices.memory_consumption() + mg_transfer.memory_consumption() + coarse_matrix.memory_consumption()); std::cout << "Multigrid objects memory consumption: " << multigrid_memory/double(1<<20) << " MiB." << std::endl; SolverControl solver_control (1000, 1e-12); SolverCG<> cg (solver_control); cg.solve (system_matrix, solution, system_rhs, preconditioner); std::cout << "Convergence in " << solver_control.last_step() << " CG iterations." << std::endl; }
template <int dim> void LaplaceProblem<dim>::output_results (const unsigned int cycle) const { DataOut<dim> data_out; data_out.attach_dof_handler (mg_dof_handler); data_out.add_data_vector (solution, "solution"); data_out.build_patches (); std::ostringstream filename; filename << "solution-" << cycle << ".vtk"; std::ofstream output (filename.str().c_str()); data_out.write_vtk (output); }
template <int dim> void LaplaceProblem<dim>::run () { for (unsigned int cycle=0; cycle<8-dim; ++cycle) { std::cout << "Cycle " << cycle << std::endl; if (cycle == 0) { GridGenerator::hyper_ball(triangulation); static const HyperBallBoundary<dim> boundary; triangulation.set_boundary (0, boundary); triangulation.refine_global (3-dim); } triangulation.refine_global (1); setup_system (); assemble_system (); assemble_multigrid (); solve (); output_results (cycle); std::cout << std::endl; }; }
int main () { try { deallog.depth_console (0); LaplaceProblem<2> laplace_problem (2); laplace_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; }
documentation generated on Wed Jul 28 23:06:04 2010 by
doxygen
1.5.6