From 70e86f56c168a091b80f7d302919bd09fb617ba3 Mon Sep 17 00:00:00 2001 From: kronbichler Date: Sat, 8 Oct 2011 15:42:30 +0000 Subject: [PATCH] Remove unfinished tutorial program. git-svn-id: https://svn.dealii.org/branches/releases/Branch-7-1@24577 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-37/Makefile | 143 -- deal.II/examples/step-37/doc/builds-on | 1 - deal.II/examples/step-37/doc/intro.dox | 406 ------ deal.II/examples/step-37/doc/kind | 1 - deal.II/examples/step-37/doc/results.dox | 311 ---- deal.II/examples/step-37/doc/tooltip | 1 - deal.II/examples/step-37/step-37.cc | 1639 ---------------------- 7 files changed, 2502 deletions(-) delete mode 100644 deal.II/examples/step-37/Makefile delete mode 100644 deal.II/examples/step-37/doc/builds-on delete mode 100644 deal.II/examples/step-37/doc/intro.dox delete mode 100644 deal.II/examples/step-37/doc/kind delete mode 100644 deal.II/examples/step-37/doc/results.dox delete mode 100644 deal.II/examples/step-37/doc/tooltip delete mode 100644 deal.II/examples/step-37/step-37.cc diff --git a/deal.II/examples/step-37/Makefile b/deal.II/examples/step-37/Makefile deleted file mode 100644 index 2c39a273a8..0000000000 --- a/deal.II/examples/step-37/Makefile +++ /dev/null @@ -1,143 +0,0 @@ -# $Id$ - - -# For the small projects Makefile, you basically need to fill in only -# four fields. -# -# The first is the name of the application. It is assumed that the -# application name is the same as the base file name of the single C++ -# file from which the application is generated. -target = step-37 - -# The second field determines whether you want to run your program in -# debug or optimized mode. The latter is significantly faster, but no -# run-time checking of parameters and internal states is performed, so -# you should set this value to `on' while you develop your program, -# and to `off' when running production computations. -debug-mode = on - - -# As third field, we need to give the path to the top-level deal.II -# directory. You need to adjust this to your needs. Since this path is -# probably the most often needed one in the Makefile internals, it is -# designated by a single-character variable, since that can be -# reference using $D only, i.e. without the parentheses that are -# required for most other parameters, as e.g. in $(target). -D = ../../ - - -# The last field specifies the names of data and other files that -# shall be deleted when calling `make clean'. Object and backup files, -# executables and the like are removed anyway. Here, we give a list of -# files in the various output formats that deal.II supports. -clean-up-files = *gmv *gnuplot *gpl *eps *pov *vtk *ucd *.d2 - - - - -# -# -# Usually, you will not need to change anything beyond this point. -# -# -# The next statement tells the `make' program where to find the -# deal.II top level directory and to include the file with the global -# settings -include $D/common/Make.global_options - - -# Since the whole project consists of only one file, we need not -# consider difficult dependencies. We only have to declare the -# libraries which we want to link to the object file. deal.II has two -# libraries: one for the debug mode version of the -# application and one for optimized mode. -libs.g := $(lib-deal2.g) -libs.o := $(lib-deal2.o) - - -# We now use the variable defined above to switch between debug and -# optimized mode to select the set of libraries to link with. Included -# in the list of libraries is the name of the object file which we -# will produce from the single C++ file. Note that by default we use -# the extension .g.o for object files compiled in debug mode and .o for -# object files in optimized mode (or whatever local default on your -# system is instead of .o) -ifeq ($(debug-mode),on) - libraries = $(target).g.$(OBJEXT) $(libs.g) -else - libraries = $(target).$(OBJEXT) $(libs.o) -endif - - -# Now comes the first production rule: how to link the single object -# file produced from the single C++ file into the executable. Since -# this is the first rule in the Makefile, it is the one `make' selects -# if you call it without arguments. -$(target)$(EXEEXT) : $(libraries) - @echo ============================ Linking $@ - @$(CXX) -o $@ $^ $(LIBS) $(LDFLAGS) - - -# To make running the application somewhat independent of the actual -# program name, we usually declare a rule `run' which simply runs the -# program. You can then run it by typing `make run'. This is also -# useful if you want to call the executable with arguments which do -# not change frequently. You may then want to add them to the -# following rule: -run: $(target)$(EXEEXT) - @echo ============================ Running $< - @./$(target)$(EXEEXT) - - -# As a last rule to the `make' program, we define what to do when -# cleaning up a directory. This usually involves deleting object files -# and other automatically created files such as the executable itself, -# backup files, and data files. Since the latter are not usually quite -# diverse, you needed to declare them at the top of this file. -clean: - -rm -f *.$(OBJEXT) *~ Makefile.dep $(target)$(EXEEXT) $(clean-up-files) - - -# Since we have not yet stated how to make an object file from a C++ -# file, we should do so now. Since the many flags passed to the -# compiler are usually not of much interest, we suppress the actual -# command line using the `at' sign in the first column of the rules -# and write the string indicating what we do instead. -./%.g.$(OBJEXT) : - @echo "==============debug========= $( $@" - @$(CXX) $(CXXFLAGS.g) -c $< -o $@ -./%.$(OBJEXT) : - @echo "==============optimized===== $( $@" - @$(CXX) $(CXXFLAGS.o) -c $< -o $@ - - -# The following statement tells make that the rules `run' and `clean' -# are not expected to produce files of the same name as Makefile rules -# usually do. -.PHONY: run clean - - -# Finally there is a rule which you normally need not care much about: -# since the executable depends on some include files from the library, -# besides the C++ application file of course, it is necessary to -# re-generate the executable when one of the files it depends on has -# changed. The following rule creates a dependency file -# `Makefile.dep', which `make' uses to determine when to regenerate -# the executable. This file is automagically remade whenever needed, -# i.e. whenever one of the cc-/h-files changed. Make detects whether -# to remake this file upon inclusion at the bottom of this file. -# -# If the creation of Makefile.dep fails, blow it away and fail -Makefile.dep: $(target).cc Makefile \ - $(shell echo $D/include/deal.II/*/*.h) - @echo ============================ Remaking $@ - @$D/common/scripts/make_dependencies $(INCLUDE) -B. $(target).cc \ - > $@ \ - || (rm -f $@ ; false) - @if test -s $@ ; then : else rm $@ ; fi - -# To make the dependencies known to `make', we finally have to include -# them: -include Makefile.dep - - diff --git a/deal.II/examples/step-37/doc/builds-on b/deal.II/examples/step-37/doc/builds-on deleted file mode 100644 index 42c2846921..0000000000 --- a/deal.II/examples/step-37/doc/builds-on +++ /dev/null @@ -1 +0,0 @@ -step-16 diff --git a/deal.II/examples/step-37/doc/intro.dox b/deal.II/examples/step-37/doc/intro.dox deleted file mode 100644 index 9d488ece7f..0000000000 --- a/deal.II/examples/step-37/doc/intro.dox +++ /dev/null @@ -1,406 +0,0 @@ -
- - -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. - - - - -

Introduction

- -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. - -

Matrix-vector product implementation

- -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: -@f{eqnarray*} -A = \sum_{\mathrm{cell}=1}^{\mathrm{n\_cells}} P_\mathrm{cell,{loc-glob}}^T A_\mathrm{cell} -P_\mathrm{cell,{loc-glob}}. -@f} -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 -@f{eqnarray*} -y &=& A\cdot x = \left(\sum_{\text{cell}=1}^{\mathrm{n\_cells}} P_\mathrm{cell,{loc-glob}}^T -A_\mathrm{cell} P_\mathrm{cell,{loc-glob}}\right) \cdot x -\\ -&=& \sum_{\mathrm{cell}=1}^{\mathrm{n\_cells}} P_\mathrm{cell,{loc-glob}}^T -A_\mathrm{cell} x_\mathrm{cell} -\\ -&=& \sum_{\mathrm{cell}=1}^{\mathrm{n\_cells}} P_\mathrm{cell,{loc-glob}}^T -y_\mathrm{cell}, -@f} -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: -@code -MatrixFree::vmult (Vector &dst, - const Vector &src) const -{ - dst = 0; - - QGauss quadrature_formula(fe.degree+1); - FEValues 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 cell_matrix (dofs_per_cell, dofs_per_cell); - Vector cell_src (dofs_per_cell), - cell_dst (dofs_per_cell); - - std::vector local_dof_indices (dofs_per_cell); - - typename DoFHandler::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_pointget_dof_indices (local_dof_indices); - - for (unsigned int i=0; isrc 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, -@f{eqnarray*} -A_\mathrm{cell} = B_\mathrm{cell}^T D_\mathrm{cell} B_\mathrm{cell}, -@f} -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 @p dofs_per_cell columns). The matrix -Dcell is diagonal and contains the values -fe_values.JxW(q) (or, rather, @p dim copies of it). - -Every numerical analyst learns in one of her first classes that for -forming a product of the form -@f{eqnarray*} -A_\mathrm{cell}\cdot x_\mathrm{cell} = B_\mathrm{cell} D_\mathrm{cell} - B_\mathrm{cell}^T \cdot x_\mathrm{cell}, -@f} -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: -@code -... - for (; cell!=endc; ++cell) - { - fe_values.reinit (cell); - - cell->get_dof_indices (local_dof_indices); - - for (unsigned int i=0; id 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 $\mathcal {O}(\mathrm{dofs\_per\_cell}^3)$ to $\mathcal -{O}(\mathrm{dofs\_per\_cell}^2)$. 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 -@f{eqnarray*} -B_\mathrm{cell} = J_\mathrm{cell} B_\mathrm{ref\_cell}, -@f} -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 -@f{eqnarray*} -A_\mathrm{cell} = B_\mathrm{cell}^T D B_\mathrm{cell} - = B_\mathrm{ref\_cell}^T J_\mathrm{cell}^T - D_\mathrm{cell} - J_\mathrm{cell} B_\mathrm{ref\_cell}, -@f} -so we calculate the product (starting the local product from the right) -@f{eqnarray*} -y_\mathrm{cell} = B_\mathrm{ref\_cell}^T J_\mathrm{cell}^T D J_\mathrm{cell} -B_\mathrm{ref\_cell} x_\mathrm{cell}, \quad -y = \sum_{\mathrm{cell}=1}^{\mathrm{n\_cells}} P_\mathrm{cell,{loc-glob}}^T -y_\mathrm{cell}. -@f} -@code -... - FEValues fe_values_reference (fe, quadrature_formula, - update_gradients); - Triangulation reference_cell; - GridGenerator::hyper_cube(reference_cell, 0., 1.); - fe_values_reference.reinit (reference_cell.begin()); - - FEValues 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 temp; - for (unsigned int d=0; ddim 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: -
    -
  • We pre-compute the inverse of the Jacobian of the transformation and - store it in an extra array. This allows us to fuse the three - operations JcellT Dcell - Jcell (apply Jacobian, multiply by - weights, apply transposed - Jacobian) into one second-rank tensor that is also symmetric (so we - only need to store half the tensor). -
  • We work on several cells at once when we apply the gradients of the - unit cell (it is always the same matrix with the reference cell - data). This allows us to replace the matrix-vector product by a - matrix-matrix product (several vectors of cell-data form a matrix), - which enables a faster implementation. Obviously, we need some adapted - data structures for that, but it isn't too hard to provide that. What - is nice is that dense matrix-matrix products are close to today's - processors' peak performance if the matrices are neither too small nor - too large — and these operations are the most expensive part in - the implementation shown here. -
- -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. - - -

Combination with multigrid

- -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. - - -

The test case

- -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. diff --git a/deal.II/examples/step-37/doc/kind b/deal.II/examples/step-37/doc/kind deleted file mode 100644 index c1d9154931..0000000000 --- a/deal.II/examples/step-37/doc/kind +++ /dev/null @@ -1 +0,0 @@ -techniques diff --git a/deal.II/examples/step-37/doc/results.dox b/deal.II/examples/step-37/doc/results.dox deleted file mode 100644 index 6e00f0b1a1..0000000000 --- a/deal.II/examples/step-37/doc/results.dox +++ /dev/null @@ -1,311 +0,0 @@ -

Results

- -

Program output

- -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 ($Q_2$) elements, we get the -following output: -@code -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. -@endcode - -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): - -@code -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. -@endcode - - -

Comparison with a sparse matrix

- -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 @p 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 (@p setup and @p 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 consumptionTime assemblyTime 50 MV, 1 CPUTime 50 MV, 4 CPUs
element orderSpMM-FSpMM-FSpMM-FSpMM-F
1299 MiB394 MiB8.09 s3.43 s5.50 s22.4 s4.30 s11.0 s
2698 MiB177 MiB12.43 s1.32 s12.0 s18.6 s9.10 s6.31 s
31295 MiB124 MiB41.1 s1.31 s21.2 s23.7 s16.0 s7.43 s
42282 MiB107 MiB117 s1.97 s40.8 s36.3 s19.7 s10.9 s
53597 MiB96.4 MiB510 s5.52 s75.7 s53.9 s29.3 s15.7 s
65679 MiB96.3 MiB2389 s26.1 s135 s79.1 s45.8 s24.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 (= @p 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 consumptionTime assemblyTime 50 MV, 4 CPUs
element orderSpMM-FSpMM-FSpMM-F
1659 MiB661 MiB18.8 s6.45 s11.0 s28.8 s
21119 MiB391 MiB15.6 s2.46 s17.1 s16.2 s
31711 MiB318 MiB17.4 s1.82 s23.1 s13.7 s
42434 MiB285 MiB24.2 s1.34 s31.1 s14.6 s
53289 MiB266 MiB35.9 s1.26 s29.6 s16.7 s
64274 MiB254 MiB58.0 s1.12 s35.9 s19.4 s
diff --git a/deal.II/examples/step-37/doc/tooltip b/deal.II/examples/step-37/doc/tooltip deleted file mode 100644 index d469776e90..0000000000 --- a/deal.II/examples/step-37/doc/tooltip +++ /dev/null @@ -1 +0,0 @@ -Matrix-free methods. Multigrid. diff --git a/deal.II/examples/step-37/step-37.cc b/deal.II/examples/step-37/step-37.cc deleted file mode 100644 index bbc6b85f0e..0000000000 --- a/deal.II/examples/step-37/step-37.cc +++ /dev/null @@ -1,1639 +0,0 @@ -/* $Id$ */ -/* Author: Katharina Kormann, Martin Kronbichler, Uppsala University, 2009 */ - -/* $Id$ */ -/* */ -/* Copyright (C) 2009, 2010 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. */ - - - // To start with the include files are more - // or less the same as in step-16: -#include -#include -#include -#include - -#include -#include -#include -#include - -#include -#include - -#include -#include -#include -#include -#include - -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include - -#include -#include - -using namespace dealii; - - - - // @sect3{Equation data} - - // We define a variable coefficient function - // for the Poisson problem. It is similar to - // the function in step-5 but we use the form - // $a(\mathbf x)=\frac{1}{0.1 + \|\bf x\|^2}$ - // instead of a discontinuous one. It is - // merely to demonstrate the possibilities of - // this implementation, rather than making - // much sense physically. -template -class Coefficient : public Function -{ - public: - Coefficient () : Function() {} - - virtual double value (const Point &p, - const unsigned int component = 0) const; - - virtual void value_list (const std::vector > &points, - std::vector &values, - const unsigned int component = 0) const; -}; - - - -template -double Coefficient::value (const Point &p, - const unsigned int /*component*/) const -{ - return 1./(0.1+p.square()); -} - - - -template -void Coefficient::value_list (const std::vector > &points, - std::vector &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 - struct ScratchData - { - ScratchData (); - ScratchData (const ScratchData &scratch); - FullMatrix solutions; - }; - - template - ScratchData::ScratchData () - : - solutions () - {} - - template - ScratchData::ScratchData (const ScratchData &) - : - solutions () - {} - - template - struct CopyData : public ScratchData - { - CopyData (); - CopyData (const CopyData &scratch); - unsigned int first_cell; - unsigned int n_dofs; - }; - - template - CopyData::CopyData () - : - ScratchData () - {} - - template - CopyData::CopyData (const CopyData &) - : - ScratchData () - {} - -} - - - - // 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 @p 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 -class MatrixFree : public Subscriptor -{ - public: - MatrixFree (); - - void reinit (const unsigned int n_dofs, - const unsigned int n_cells, - const FullMatrix &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 &local_dof_indices); - void set_derivative_data (const unsigned int cell_no, - const unsigned int quad_point, - const Transformation &trans_in); - - template - void vmult (Vector &dst, - const Vector &src) const; - template - void Tvmult (Vector &dst, - const Vector &src) const; - template - void vmult_add (Vector &dst, - const Vector &src) const; - template - void Tvmult_add (Vector &dst, - const Vector &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 - // @p 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 >::const_iterator - CellChunkIterator; - template - void local_vmult (CellChunkIterator cell_range, - WorkStreamData::ScratchData &scratch, - WorkStreamData::CopyData ©, - const Vector &src) const; - - template - void - copy_local_to_global (const WorkStreamData::CopyData ©, - Vector &dst) const; - - FullMatrix B_ref_cell; - Table<2,unsigned int> indices_local_to_global; - Table<2,Transformation> derivatives; - - ConstraintMatrix constraints; - - mutable Vector 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 > chunks; - } matrix_sizes; -}; - - - - // This is the constructor of the @p - // MatrixFree class. All it does is to - // subscribe to the general deal.II @p - // 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 -MatrixFree::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 -unsigned int -MatrixFree::m () const -{ - return matrix_sizes.n_dofs; -} - - - -template -unsigned int -MatrixFree::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 -ConstraintMatrix & -MatrixFree::get_constraints () -{ - return constraints; -} - - - - // The following function takes a vector of - // local dof indices on cell level and writes - // the data into the - // @p 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 @p false. -template -void MatrixFree:: -set_local_dof_indices (const unsigned int cell_no, - const std::vector &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 -void MatrixFree:: -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 @p 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 - // @f{eqnarray*} - // P^T_\mathrm{cell,local-global} A_\mathrm{cell} - // P_\mathrm{cell,local-global} x - // @f} - // where - // @f{eqnarray*} - // A_\mathrm{cell} = - // B_\mathrm{ref\_cell}^T J_\mathrm{cell}^T - // D_\mathrm{cell} - // J_\mathrm{cell} B_\mathrm{ref\_cell} - // @f} - // and Pcell,local-global - // is the transformation from local to global - // indices. - // - // To do this, we would have to do the - // following steps: - //
    - //
  1. Form $x_\mathrm{cell} = - // P_\mathrm{cell,local-global} x$. This is - // done by using the command - // ConstraintMatrix::get_dof_values. - //
  2. Form $x_1 = B_\mathrm{ref\_cell} - // x_\mathrm{cell}$. The vector - // x1 contains the - // reference cell gradient to the local - // cell vector. - //
  3. Form $x_2 = J_\mathrm{cell}^T - // D_\mathrm{cell} J_\mathrm{cell} - // x_1$. This is a block-diagonal - // operation, with the block size equal to - // @p 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 - // @p JxW values and the - // inhomogeneous coefficient. - //
  4. Form $y_\mathrm{cell} = - // B_\mathrm{ref\_cell}^T x_2$. This gives - // the local result of the matrix-vector - // product. - //
  5. Form $y \leftarrow y + - // P_\mathrm{cell,local-global}^T - // y_\mathrm{cell}$. 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 - // @p 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. - //
- // The steps 1 to 4 can be done in %parallel - // by multiple processes. - - // 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 @p - // dofs_per_cell rows and columns. This is an - // operation that can be much better - // optimized than matrix-vector products. The - // functions @p FullMatrix::mmult and - // @p FullMatrix::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 @p - // scratch.solutions and @p 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 -template -void -MatrixFree:: -local_vmult (CellChunkIterator cell_range, - WorkStreamData::ScratchData &scratch, - WorkStreamData::CopyData ©, - const Vector &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 -template -void -MatrixFree:: -copy_local_to_global (const WorkStreamData::CopyData ©, - Vector &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 @p vmult function that is - // called externally: In addition to what we - // do in a @p vmult_add function, we set the - // destination to zero first. -template -template -void -MatrixFree::vmult (Vector &dst, - const Vector &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 @p vmult_add operation. -template -template -void -MatrixFree::Tvmult (Vector &dst, - const Vector &src) const -{ - dst = 0; - Tvmult_add (dst,src); -} - - - -template -template -void -MatrixFree::Tvmult_add (Vector &dst, - const Vector &src) const -{ - vmult_add (dst,src); -} - - - - // This is the @p vmult_add function that - // multiplies the matrix with vector @p src - // and adds the result to vector @p 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 - // @ref threads module). The subdivision into - // chunks will be performed in the reinit - // function and is stored in the field @p - // matrix_sizes.chunks. What the rather - // cryptic command to @p 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, std_cxx1x::_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 @p - // local_vmult function are bound: to @p this - // and a constant reference to the @p src in - // the first case, and @p this and a - // reference to the output vector in the - // second. Similarly, the placeholder - // @p _1 argument in the - // @p copy_local_to_global function - // sets the first explicit argument of that - // function, which is of class - // @p 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 -template -void -MatrixFree::vmult_add (Vector &dst, - const Vector &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:: - template local_vmult, - this, std_cxx1x::_1, std_cxx1x::_2, std_cxx1x::_3, boost::cref(src)), - std_cxx1x::bind(&MatrixFree:: - template copy_local_to_global, - this, std_cxx1x::_1, boost::ref(dst)), - WorkStreamData::ScratchData(), - WorkStreamData::CopyData(), - 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 -void MatrixFree:: -reinit (const unsigned int n_dofs_in, - const unsigned int n_cells_in, - const FullMatrix &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 @p - // 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 chunk; - for (unsigned int i=0; i 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 -void -MatrixFree::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 -number -MatrixFree::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 -void -MatrixFree::calculate_diagonal() const -{ - diagonal_values.reinit (matrix_sizes.n_dofs); - std::vector calculation (matrix_sizes.n); - for (unsigned int cell=0; cell -std::size_t MatrixFree::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; -} - - - - // @sect3{Laplace operator implementation} - - // 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 @p transform operation - // needed by the @p 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 - // @p double %numbers. Since we also want to - // use @p float %numbers for the multigrid - // preconditioner (in order to save memory - // and computing time), we manually implement - // this operator. Note that @p 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 -class LaplaceOperator -{ - public: - LaplaceOperator (); - - LaplaceOperator (const Tensor<2,dim> &tensor); - - void transform (number * result) const; - - LaplaceOperator& - operator = (const Tensor<2,dim> &tensor); - - unsigned int memory_consumption () const; - - private: - number transformation[dim*(dim+1)/2]; -}; - -template -LaplaceOperator::LaplaceOperator() -{} - - -template -LaplaceOperator::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 @p 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 @p if clause at this position - // in the code at the innermost loop, but - // note once again that @p dim is known when - // this piece of code is compiled, so the - // compiler can optimize away the @p if - // statement (and actually even inline these - // few lines of code into the @p MatrixFree - // class). -template -void LaplaceOperator::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 - // @p 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 -LaplaceOperator& -LaplaceOperator::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 -unsigned int -LaplaceOperator::memory_consumption () const -{ - return sizeof(*this); -} - - - - // @sect3{LaplaceProblem class} - - // This class is based on the same - // class in step-16. However, we - // replaced the SparseMatrix - // class by our matrix-free - // implementation, which means that - // we can also skip the sparsity - // patterns. -template -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 triangulation; - FE_Q fe; - MGDoFHandler mg_dof_handler; - - MatrixFree > system_matrix; - typedef MatrixFree > MatrixFreeType; - MGLevelObject mg_matrices; - FullMatrix coarse_matrix; - - Vector solution; - Vector system_rhs; -}; - - - -template -LaplaceProblem::LaplaceProblem (const unsigned int degree) - : - fe (degree), - mg_dof_handler (triangulation) -{} - - - - // @sect4{LaplaceProblem::setup_system} - - // 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 -void LaplaceProblem::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 quadrature_formula(fe.degree+1); - FEValues fe_values_reference (fe, quadrature_formula, - update_gradients); - Triangulation reference_cell; - GridGenerator::hyper_cube (reference_cell, 0, 1); - fe_values_reference.reinit (reference_cell.begin()); - FullMatrix ref_cell_gradients (fe.dofs_per_cell, - quadrature_formula.size()*dim); - for (unsigned int i=0; i(), - 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::type dirichlet_boundary; - ZeroFunction homogeneous_dirichlet_bc (1); - dirichlet_boundary[0] = &homogeneous_dirichlet_bc; - std::vector > boundary_indices(triangulation.n_levels()); - MGTools::make_boundary_list (mg_dof_handler, - dirichlet_boundary, - boundary_indices); - for (unsigned int level=0;level::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)); -} - - - - // @sect4{LaplaceProblem::assemble_system} - - // 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 @p - // 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 -void LaplaceProblem::assemble_system () -{ - QGauss quadrature_formula(fe.degree+1); - MappingQ mapping (fe.degree); - FEValues 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 local_dof_indices (dofs_per_cell); - const Coefficient coefficient; - std::vector coefficient_values (n_q_points); - - unsigned int cell_no = 0; - - typename DoFHandler::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 -void LaplaceProblem::assemble_multigrid () -{ - coarse_matrix = 0; - QGauss quadrature_formula(fe.degree+1); - MappingQ mapping (fe.degree); - FEValues 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 local_dof_indices (dofs_per_cell); - const Coefficient coefficient; - std::vector coefficient_values (n_q_points); - - std::vector cell_no(triangulation.n_levels()); - typename MGDoFHandler::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 -void LaplaceProblem::solve () -{ - GrowingVectorMemory<> vector_memory; - - MGTransferPrebuilt > mg_transfer; - mg_transfer.build_matrices(mg_dof_handler); - - MGCoarseGridHouseholder > mg_coarse; - mg_coarse.initialize(coarse_matrix); - - typedef PreconditionChebyshev > SMOOTHER; - MGSmootherPrecondition > - 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 - // $[\lambda_{\max}/10,\lambda_{\max}]$. 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 > - mg_matrix(&mg_matrices); - - Multigrid > mg(mg_dof_handler, - mg_matrix, - mg_coarse, - mg_transfer, - mg_smoother, - mg_smoother); - PreconditionMG, - MGTransferPrebuilt > > - 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; -} - - - - // @sect4{LaplaceProblem::output_results} - - // 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 -void LaplaceProblem::output_results (const unsigned int cycle) const -{ - DataOut 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); -} - - - - // @sect4{LaplaceProblem::run} - - // 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 -void LaplaceProblem::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 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; - }; -} - - - - // @sect3{The main function} - - // This 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; -} -- 2.39.5