From 7478b133616a6d839e5ef0d64e23004e2c55f287 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Mon, 31 Aug 2009 15:33:08 +0000 Subject: [PATCH] Add a first version of the step-37 tutorial program that implements a matrix-free method for the Laplace equation and combines it with multigrid. git-svn-id: https://svn.dealii.org/trunk@19355 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-37/Makefile | 156 ++++ deal.II/examples/step-37/doc/intro.dox | 14 + deal.II/examples/step-37/doc/results.dox | 47 ++ deal.II/examples/step-37/step-37.cc | 885 +++++++++++++++++++++++ 4 files changed, 1102 insertions(+) create mode 100644 deal.II/examples/step-37/Makefile create mode 100644 deal.II/examples/step-37/doc/intro.dox create mode 100644 deal.II/examples/step-37/doc/results.dox create 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 new file mode 100644 index 0000000000..feef6229a0 --- /dev/null +++ b/deal.II/examples/step-37/Makefile @@ -0,0 +1,156 @@ +# $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 = $(basename $(shell echo step-*.cc)) + +# 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 + + + + +# +# +# Usually, you will not need to change anything beyond this point. +# +# +# The next statement tell 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, and there need +# to be two sets of libraries: one for the debug mode version of the +# application and one for the optimized mode. Here we have selected +# the versions for 2d. Note that the order in which the libraries are +# given here is important and that your applications won't link +# properly if they are given in another order. +# +# You may need to augment the lists of libraries when compiling your +# program for other dimensions, or when using third party libraries +libs.g = $(lib-deal2-2d.g) \ + $(lib-lac.g) \ + $(lib-base.g) +libs.o = $(lib-deal2-2d.o) \ + $(lib-deal2-3d.o) \ + $(lib-lac.o) \ + $(lib-base.o) + + +# We now use the variable defined above which 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 the 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) : $(libraries) + @echo ============================ Linking $@ + @$(CXX) -o $@$(EXEEXT) $^ $(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) + @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========= $( $@ \ + || (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/intro.dox b/deal.II/examples/step-37/doc/intro.dox new file mode 100644 index 0000000000..0580cb04cd --- /dev/null +++ b/deal.II/examples/step-37/doc/intro.dox @@ -0,0 +1,14 @@ + +

Introduction

+ +This example shows how a matrix-free method can be implemented for a +second-order Poisson equation. + +

Implementation of a matrix-vector product

+TODO: Describe how one gets to the piece of code we use below. + +

Combination with multigrid

+TODO: Explain which kind of smoother we can use when we do not have +access to the matrix entries. + +

Parallization issues

diff --git a/deal.II/examples/step-37/doc/results.dox b/deal.II/examples/step-37/doc/results.dox new file mode 100644 index 0000000000..dfc70b7f4a --- /dev/null +++ b/deal.II/examples/step-37/doc/results.dox @@ -0,0 +1,47 @@ +

Results

+ +Since this example solves the same problem as @ref step_5 "step-5", we +refer to the graphical output there. The only difference between the two is +the solver and the implementation of the matrix-vector products. + +The output produced by this program is the following: +@code +Cycle 0 +Number of degrees of freedom: 337 +System matrix memory consumption: 0.02781 MBytes. +Multigrid objects memory consumption: 0.04799 MBytes. +Convergence in 8 CG iterations. + +Cycle 1 +Number of degrees of freedom: 1313 +System matrix memory consumption: 0.103 MBytes. +Multigrid objects memory consumption: 0.171 MBytes. +Convergence in 9 CG iterations. + +Cycle 2 +Number of degrees of freedom: 5185 +System matrix memory consumption: 0.4019 MBytes. +Multigrid objects memory consumption: 0.6524 MBytes. +Convergence in 9 CG iterations. + +Cycle 3 +Number of degrees of freedom: 20609 +System matrix memory consumption: 1.594 MBytes. +Multigrid objects memory consumption: 2.562 MBytes. +Convergence in 9 CG iterations. + +Cycle 4 +Number of degrees of freedom: 82177 +System matrix memory consumption: 6.357 MBytes. +Multigrid objects memory consumption: 10.18 MBytes. +Convergence in 9 CG iterations. + +Cycle 5 +Number of degrees of freedom: 328193 +System matrix memory consumption: 25.39 MBytes. +Multigrid objects memory consumption: 40.58 MBytes. +Convergence in 9 CG iterations. +@endcode + +As in step-16, we see that the number of CG iterations stays constant with +increasing number of degrees of freedom. \ No newline at end of file diff --git a/deal.II/examples/step-37/step-37.cc b/deal.II/examples/step-37/step-37.cc new file mode 100644 index 0000000000..83d9fb172e --- /dev/null +++ b/deal.II/examples/step-37/step-37.cc @@ -0,0 +1,885 @@ +/* $Id$ */ +/* Author: Martin Kronbichler, Uppsala University, 2009 */ + +/* $Id$ */ +/* */ +/* 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. */ + + + // The include files are mostly similar to + // the ones 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 +#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. As a + // difference, we use the formulation + // $\frac{1}{0.1 + \|\bf x\|^2}$ instead of + // a discontinuous one. It is merely to + // demonstrate the possibilities of this + // implemenation, rather than being + // physically reasonable. +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 +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_on_subrange (const unsigned int first_cell, + const unsigned int last_cell, + Vector &dst, + const Vector &src) const; + + 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; + + std::size_t memory_consumption () const; + +private: + FullMatrix small_matrix; + ConstraintMatrix constraints; + + Table<2,unsigned int> indices_local_to_global; + Table<2,Transformation> derivatives; + + mutable Vector diagonal_values; + mutable bool diagonal_is_calculated; + + struct SmallMatrixData + { + unsigned int m; + unsigned int n; + unsigned int n_points; + unsigned int n_comp; + }; + unsigned int n_dofs, n_cols, n_cells; + SmallMatrixData matrix_data; +}; + + + +template +MatrixFree::MatrixFree () +: + Subscriptor() +{} + + + +template +void MatrixFree:: +reinit (const unsigned int n_dofs_in, + const unsigned int n_cells_in, + const FullMatrix &small_matrix_in, + const unsigned int n_points_per_cell) +{ + n_dofs = n_dofs_in; + n_cells = n_cells_in; + small_matrix = small_matrix_in; + matrix_data.m = small_matrix.m(); + matrix_data.n = small_matrix.n(); + matrix_data.n_points = n_points_per_cell; + matrix_data.n_comp = small_matrix.n()/matrix_data.n_points; + + Assert(matrix_data.n_comp * n_points_per_cell == small_matrix.n(), + ExcInternalError()); + + derivatives.reinit (n_cells, n_points_per_cell); + indices_local_to_global.reinit (n_cells, small_matrix.m()); + diagonal_is_calculated = false; +} + + + +template +void +MatrixFree::clear () +{ + n_dofs = 0; + n_cells = 0; + small_matrix.reinit(0,0); + derivatives.reinit (0,0); + indices_local_to_global.reinit(0,0); + diagonal_values.reinit (0); + constraints.clear(); + diagonal_is_calculated = false; +} + + + +template +unsigned int +MatrixFree::m () const +{ + return n_dofs; +} + + + +template +unsigned int +MatrixFree::n () const +{ + return n_dofs; +} + + + +template +ConstraintMatrix & +MatrixFree::get_constraints () +{ + return constraints; +} + + + +template +void MatrixFree:: +set_local_dof_indices (const unsigned int cell_no, + const std::vector &local_dof_indices) +{ + Assert (local_dof_indices.size() == matrix_data.m, + ExcDimensionMismatch(local_dof_indices.size(), + matrix_data.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) +{ + derivatives(cell_no,quad_point) = trans_in; + diagonal_is_calculated = false; +} + + + +template +template +void +MatrixFree:: +vmult_on_subrange (const unsigned int first_cell, + const unsigned int last_cell, + Vector &dst, + const Vector &src) const +{ + FullMatrix solution_cells, solution_points; + + const unsigned int n_chunks = (last_cell-first_cell)/100 + 1; + const unsigned int chunk_size = + (last_cell-first_cell)/n_chunks + ((last_cell-first_cell)%n_chunks>0); + + for (unsigned int k=first_cell; klast_cell ? last_cell-k : chunk_size; + + solution_cells.reinit (current_chunk_size,matrix_data.m, true); + solution_points.reinit (current_chunk_size,matrix_data.n, true); + + for (unsigned int i=0; i +template +void +MatrixFree::vmult (Vector &dst, + const Vector &src) const +{ + dst = 0; + vmult_add (dst, src); +} + + + +template +template +void +MatrixFree::Tvmult (Vector &dst, + const Vector &src) const +{ + dst = 0; + Tvmult_add (dst,src); +} + + + +template +template +void +MatrixFree::vmult_add (Vector &dst, + const Vector &src) const +{ + Vector src_copy (src); + constraints.distribute(src_copy); + + vmult_on_subrange (0, n_cells, dst, src_copy); + constraints.condense (dst); + + // Need to do this in order to be + // consistent even at constrained + // dofs. Need to find a better solution in + // the future (e.g. by switching to smaller + // vectors that do not contain any + // constrained entries). + for (unsigned int i=0; i +template +void +MatrixFree::Tvmult_add (Vector &dst, + const Vector &src) const +{ + vmult_add (dst,src); +} + + + +template +number +MatrixFree::el (const unsigned int row, + const unsigned int col) const +{ + Assert (row == col, ExcNotImplemented()); + + if (diagonal_is_calculated == false) + { + diagonal_values.reinit (n_dofs); + std::vector calculation (matrix_data.n_comp); + 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() + + small_matrix.memory_consumption() + sizeof(*this); + return glob_size; +} + + + + // @sect3{Laplace operator.} + + // This implements the local action of a + // Laplace preconditioner. +template +class LaplaceOperator +{ +public: + LaplaceOperator (); + + LaplaceOperator (const Tensor<2,dim> &tensor); + + void transform (number * result) const; + + LaplaceOperator& + operator = (const Tensor<2,dim> &tensor); + + number transformation[dim][dim]; +}; + +template +LaplaceOperator::LaplaceOperator() +{} + +template +LaplaceOperator::LaplaceOperator(const Tensor<2,dim> &tensor) +{ + *this = tensor; +} + +template +void LaplaceOperator::transform (number* result) const +{ + number temp_result[dim]; + for (unsigned int d=0; d +LaplaceOperator& +LaplaceOperator::operator=(const Tensor<2,dim> &tensor) +{ + for (unsigned int d=0;d class by our + // matrix-free implementation, which means + // that we can 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) +{} + + + + // This is the function of step-16 with + // relevant changes due to the MatrixFree + // class. +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_values2 (fe, quadrature_formula, + update_gradients); + Triangulation tria; + GridGenerator::hyper_cube (tria, 0, 1); + fe_values2.reinit (tria.begin()); + FullMatrix data_matrix (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: " + << (double)system_matrix.memory_consumption()*std::pow(2.,-20.) << " MBytes." + << std::endl; + + solution.reinit (mg_dof_handler.n_dofs()); + system_rhs.reinit (mg_dof_handler.n_dofs()); + + // Initialize the matrices for the + // multigrid method on all the levels. + 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)); +} + + + +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_gradients | + 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); + + system_matrix.set_local_dof_indices (cell_no, local_dof_indices); + for (unsigned int q=0; q +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); + + // Initialize the smoother with our level + // matrices and the required, additional + // data for the Chebyshev smoother. Use a + // higher polynomial degree for higher + // order elements, since smoothing gets + // more difficult then. Smooth out a + // range of + // $[\lambda_{\max}/8,\lambda_{\max}]$. + typename SMOOTHER::AdditionalData smoother_data; + smoother_data.smoothing_range = 8.; + smoother_data.degree = fe.degree+1; + 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); + + double multigrid_memory = + (double)mg_matrices.memory_consumption() + + (double)mg_transfer.memory_consumption() + + (double)coarse_matrix.memory_consumption(); + + std::cout << "Multigrid objects memory consumption: " + << multigrid_memory*std::pow(2.,-20.) + << " MBytes." + << std::endl; + + // Finally, create the solver + // object and solve the system + 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 + // do a 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); +} + + + +template +void LaplaceProblem::run () +{ + for (unsigned int cycle=0; cycle<6; ++cycle) + { + std::cout << "Cycle " << cycle << std::endl; + + if (cycle == 0) + { + // Generate a simple hyperball grid. + GridGenerator::hyper_ball(triangulation); + static const HyperBallBoundary boundary; + triangulation.set_boundary (0, boundary); + triangulation.refine_global (0); + } + triangulation.refine_global (1); + setup_system (); + assemble_system (); + assemble_multigrid (); + solve (); + output_results (cycle); + std::cout << std::endl; + }; +} + + + +int main () +{ + deallog.depth_console (0); + LaplaceProblem<2> laplace_problem (2); + laplace_problem.run (); + + return 0; +} -- 2.39.5