From: Denis Davydov Date: Tue, 4 Oct 2016 13:12:27 +0000 (+0200) Subject: make PreconditionJacobi work with MGSmootherPrecondition by adding a typedef for... X-Git-Tag: v8.5.0-rc1~579^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=13535d1bb26035fe5eca67b1b0f39dec1660174b;p=dealii.git make PreconditionJacobi work with MGSmootherPrecondition by adding a typedef for AdditionalData --- diff --git a/include/deal.II/lac/precondition.h b/include/deal.II/lac/precondition.h index d21c71a9e5..6998e0f332 100644 --- a/include/deal.II/lac/precondition.h +++ b/include/deal.II/lac/precondition.h @@ -452,7 +452,7 @@ protected: * @ref ConceptRelaxationType "relaxation concept". * * @code - * // Declare related objects + * // Declare related objects * * SparseMatrix A; * Vector x; @@ -461,10 +461,10 @@ protected: * * //...initialize and build A * - * // Define and initialize preconditioner + * // Define and initialize preconditioner: * * PreconditionJacobi > precondition; - * precondition.initialize (A, .6); + * precondition.initialize (A, PreconditionJacobi >::AdditionalData(.6)); * * solver.solve (A, x, b, precondition); * @endcode @@ -475,6 +475,11 @@ template > class PreconditionJacobi : public PreconditionRelaxation { public: + /** + * A typedef to the base class AdditionalData. + */ + typedef typename PreconditionRelaxation::AdditionalData AdditionalData; + /** * Apply preconditioner. */ @@ -529,7 +534,7 @@ public: * @ref ConceptRelaxationType "relaxation concept". * * @code - * // Declare related objects + * // Declare related objects * * SparseMatrix A; * Vector x; @@ -538,10 +543,10 @@ public: * * //...initialize and build A * - * // Define and initialize preconditioner + * // Define and initialize preconditioner * * PreconditionSOR > precondition; - * precondition.initialize (A, .6); + * precondition.initialize (A, PreconditionSOR >::AdditionalData(.6)); * * solver.solve (A, x, b, precondition); * @endcode @@ -552,6 +557,11 @@ template > class PreconditionSOR : public PreconditionRelaxation { public: + /** + * A typedef to the base class AdditionalData. + */ + typedef typename PreconditionRelaxation::AdditionalData AdditionalData; + /** * Apply preconditioner. */ @@ -596,10 +606,10 @@ public: * * //...initialize and build A * - * // Define and initialize preconditioner + * // Define and initialize preconditioner * * PreconditionSSOR > precondition; - * precondition.initialize (A, .6); + * precondition.initialize (A, PreconditionSSOR >::AdditionalData(.6)); * * solver.solve (A, x, b, precondition); * @endcode @@ -610,6 +620,11 @@ template > class PreconditionSSOR : public PreconditionRelaxation { public: + /** + * A typedef to the base class AdditionalData. + */ + typedef typename PreconditionRelaxation::AdditionalData AdditionalData; + /** * Declare type for container size. */ diff --git a/tests/multigrid/step-16-07.cc b/tests/multigrid/step-16-07.cc new file mode 100644 index 0000000000..b5b0014c74 --- /dev/null +++ b/tests/multigrid/step-16-07.cc @@ -0,0 +1,646 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +// same as step-16 but with PreconditionJacobi + MGSmootherPrecondition. + +#include "../tests.h" +#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 + +#include +#include + +using namespace dealii; + +template +class LaplaceProblem +{ +public: + LaplaceProblem (const unsigned int deg); + void run (); + +private: + void setup_system (); + void assemble_system (); + void assemble_multigrid (); + void solve (); + void refine_grid (); + void output_results (const unsigned int cycle) const; + + Triangulation triangulation; + FE_Q fe; + DoFHandler mg_dof_handler; + + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + + ConstraintMatrix constraints; + + Vector solution; + Vector system_rhs; + + const unsigned int degree; + + MGLevelObject mg_sparsity_patterns; + MGLevelObject > mg_matrices; + MGLevelObject > mg_interface_matrices; + MGConstrainedDoFs mg_constrained_dofs; +}; + + +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) const +{ + if (p.square() < 0.5*0.5) + return 20; + else + return 1; +} + + + +template +void Coefficient::value_list (const std::vector > &points, + std::vector &values, + const unsigned int component) const +{ + const unsigned int n_points = points.size(); + + Assert (values.size() == n_points, + ExcDimensionMismatch (values.size(), n_points)); + + Assert (component == 0, + ExcIndexRange (component, 0, 1)); + + for (unsigned int i=0; i::value (points[i]); +} + + +template +LaplaceProblem::LaplaceProblem (const unsigned int degree) + : + triangulation (Triangulation:: + limit_level_difference_at_vertices), + fe (degree), + mg_dof_handler (triangulation), + degree(degree) +{} + + +template +void LaplaceProblem::setup_system () +{ + mg_dof_handler.distribute_dofs(fe); + mg_dof_handler.distribute_mg_dofs (fe); + deallog << "Number of degrees of freedom: " + << mg_dof_handler.n_dofs(); + + for (unsigned int l=0; l&>(mg_dof_handler), + sparsity_pattern); + + solution.reinit (mg_dof_handler.n_dofs()); + system_rhs.reinit (mg_dof_handler.n_dofs()); + + constraints.clear (); + DoFTools::make_hanging_node_constraints (mg_dof_handler, constraints); + typename FunctionMap::type dirichlet_boundary; + ZeroFunction homogeneous_dirichlet_bc (1); + dirichlet_boundary[0] = &homogeneous_dirichlet_bc; + MappingQGeneric mapping(1); + VectorTools::interpolate_boundary_values (mapping, + mg_dof_handler, + dirichlet_boundary, + constraints); + constraints.close (); + constraints.condense (sparsity_pattern); + sparsity_pattern.compress(); + system_matrix.reinit (sparsity_pattern); + + mg_constrained_dofs.clear(); + mg_constrained_dofs.initialize(mg_dof_handler, dirichlet_boundary); + const unsigned int n_levels = triangulation.n_levels(); + + mg_interface_matrices.resize(0, n_levels-1); + mg_interface_matrices.clear_elements (); + mg_matrices.resize(0, n_levels-1); + mg_matrices.clear_elements (); + mg_sparsity_patterns.resize(0, n_levels-1); + + for (unsigned int level=0; level +void LaplaceProblem::assemble_system () +{ + const QGauss quadrature_formula(degree+1); + + FEValues fe_values (fe, quadrature_formula, + update_values | 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(); + + FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); + Vector cell_rhs (dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); + + const Coefficient coefficient; + std::vector coefficient_values (n_q_points); + + typename DoFHandler::active_cell_iterator + cell = mg_dof_handler.begin_active(), + endc = mg_dof_handler.end(); + for (; cell!=endc; ++cell) + { + cell_matrix = 0; + cell_rhs = 0; + + fe_values.reinit (cell); + + coefficient.value_list (fe_values.get_quadrature_points(), + coefficient_values); + + for (unsigned int q_point=0; q_pointget_dof_indices (local_dof_indices); + constraints.distribute_local_to_global (cell_matrix, cell_rhs, + local_dof_indices, + system_matrix, system_rhs); + } +} + + +template +void LaplaceProblem::assemble_multigrid () +{ + QGauss quadrature_formula(1+degree); + + FEValues fe_values (fe, quadrature_formula, + update_values | 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(); + + FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); + + const Coefficient coefficient; + std::vector coefficient_values (n_q_points); + + std::vector boundary_constraints (triangulation.n_levels()); + ConstraintMatrix empty_constraints; + for (unsigned int level=0; level::cell_iterator cell = mg_dof_handler.begin(), + endc = mg_dof_handler.end(); + + for (; cell!=endc; ++cell) + { + cell_matrix = 0; + fe_values.reinit (cell); + + coefficient.value_list (fe_values.get_quadrature_points(), + coefficient_values); + + for (unsigned int q_point=0; q_pointget_mg_dof_indices (local_dof_indices); + + boundary_constraints[cell->level()] + .distribute_local_to_global (cell_matrix, + local_dof_indices, + mg_matrices[cell->level()]); + + // The next step is again slightly more + // obscure (but explained in the @ref + // mg_paper): We need the remainder of + // the operator that we just copied + // into the mg_matrices + // object, namely the part on the + // interface between cells at the + // current level and cells one level + // coarser. This matrix exists in two + // directions: for interior DoFs (index + // $i$) of the current level to those + // sitting on the interface (index + // $j$), and the other way around. Of + // course, since we have a symmetric + // operator, one of these matrices is + // the transpose of the other. + // + // The way we assemble these matrices + // is as follows: since the are formed + // from parts of the local + // contributions, we first delete all + // those parts of the local + // contributions that we are not + // interested in, namely all those + // elements of the local matrix for + // which not $i$ is an interface DoF + // and $j$ is not. The result is one of + // the two matrices that we are + // interested in, and we then copy it + // into the + // mg_interface_matrices + // object. The + // boundary_interface_constraints + // object at the same time makes sure + // that we delete contributions from + // all degrees of freedom that are not + // only on the interface but also on + // the external boundary of the domain. + // + // The last part to remember is how to + // get the other matrix. Since it is + // only the transpose, we will later + // (in the solve() + // function) be able to just pass the + // transpose matrix where necessary. + const unsigned int lvl = cell->level(); + + for (unsigned int i=0; ilevel()]); + } +} + + + +// @sect4{LaplaceProblem::solve} + +// This is the other function that is +// significantly different in support of the +// multigrid solver (or, in fact, the +// preconditioner for which we use the +// multigrid method). +// +// Let us start out by setting up two of the +// components of multilevel methods: transfer +// operators between levels, and a solver on +// the coarsest level. In finite element +// methods, the transfer operators are +// derived from the finite element function +// spaces involved and can often be computed +// in a generic way independent of the +// problem under consideration. In that case, +// we can use the MGTransferPrebuilt class +// that, given the constraints on the global +// level and an DoFHandler object computes +// the matrices corresponding to these +// transfer operators. +// +// The second part of the following lines +// deals with the coarse grid solver. Since +// our coarse grid is very coarse indeed, we +// decide for a direct solver (a Householder +// decomposition of the coarsest level +// matrix), even if its implementation is not +// particularly sophisticated. If our coarse +// mesh had many more cells than the five we +// have here, something better suited would +// obviously be necessary here. +template +void LaplaceProblem::solve () +{ + MGTransferPrebuilt > mg_transfer(mg_constrained_dofs); + mg_transfer.build_matrices(mg_dof_handler); + + FullMatrix coarse_matrix; + coarse_matrix.copy_from (mg_matrices[0]); + MGCoarseGridHouseholder<> coarse_grid_solver; + coarse_grid_solver.initialize (coarse_matrix); + + typedef PreconditionJacobi > Smoother; + MGSmootherPrecondition, Smoother, Vector > + mg_smoother; + mg_smoother.initialize(mg_matrices, typename Smoother::AdditionalData(0.678)); + mg_smoother.set_steps(2); + mg_smoother.set_symmetric(true); + + mg::Matrix<> mg_matrix(mg_matrices); + mg::Matrix<> mg_interface_up(mg_interface_matrices); + mg::Matrix<> mg_interface_down(mg_interface_matrices); + + Multigrid > mg(mg_dof_handler, + mg_matrix, + coarse_grid_solver, + mg_transfer, + mg_smoother, + mg_smoother); + mg.set_edge_matrices(mg_interface_down, mg_interface_up); + + PreconditionMG, MGTransferPrebuilt > > + preconditioner(mg_dof_handler, mg, mg_transfer); + + SolverControl solver_control (1000, 1e-12); + SolverCG<> cg (solver_control); + + solution = 0; + + cg.solve (system_matrix, solution, system_rhs, + preconditioner); + constraints.distribute (solution); + + deallog << " " << solver_control.last_step() + << " CG iterations needed to obtain convergence." + << std::endl; +} + + + +// @sect4{Postprocessing} + +// The following two functions postprocess a +// solution once it is computed. In +// particular, the first one refines the mesh +// at the beginning of each cycle while the +// second one outputs results at the end of +// each such cycle. The functions are almost +// unchanged from those in step-6, with the +// exception of two minor differences: The +// KellyErrorEstimator::estimate function +// wants an argument of type DoFHandler, not +// DoFHandler, and so we have to cast from +// derived to base class; and we generate +// output in VTK format, to use the more +// modern visualization programs available +// today compared to those that were +// available when step-6 was written. +template +void LaplaceProblem::refine_grid () +{ + Vector estimated_error_per_cell (triangulation.n_active_cells()); + + KellyErrorEstimator::estimate (static_cast&>(mg_dof_handler), + QGauss(3), + typename FunctionMap::type(), + solution, + estimated_error_per_cell); + GridRefinement::refine_and_coarsen_fixed_number (triangulation, + estimated_error_per_cell, + 0.3, 0.03); + triangulation.execute_coarsening_and_refinement (); +} + + + +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} + +// Like several of the functions above, this +// is almost exactly a copy of of the +// corresponding function in step-6. The only +// difference is the call to +// assemble_multigrid that takes +// care of forming the matrices on every +// level that we need in the multigrid +// method. +template +void LaplaceProblem::run () +{ + for (unsigned int cycle=0; cycle<8; ++cycle) + { + deallog << "Cycle " << cycle << ':' << std::endl; + + if (cycle == 0) + { + GridGenerator::hyper_ball (triangulation); + + static const HyperBallBoundary boundary; + triangulation.set_boundary (0, boundary); + + triangulation.refine_global (1); + } + else + refine_grid (); + + + deallog << " Number of active cells: " + << triangulation.n_active_cells() + << std::endl; + + setup_system (); + + deallog << " Number of degrees of freedom: " + << mg_dof_handler.n_dofs() + << " (by level: "; + for (unsigned int level=0; level laplace_problem(1); + 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; +} diff --git a/tests/multigrid/step-16-07.output b/tests/multigrid/step-16-07.output new file mode 100644 index 0000000000..9847ced8bc --- /dev/null +++ b/tests/multigrid/step-16-07.output @@ -0,0 +1,57 @@ + +DEAL::Cycle 0: +DEAL:: Number of active cells: 20 +DEAL::Number of degrees of freedom: 25 L0: 8 L1: 25 +DEAL:: Number of degrees of freedom: 25 (by level: 8, 25) +DEAL:cg::Starting value 0.5107 +DEAL:cg::Convergence step 5 value 0 +DEAL:: 5 CG iterations needed to obtain convergence. +DEAL::Cycle 1: +DEAL:: Number of active cells: 44 +DEAL::Number of degrees of freedom: 57 L0: 8 L1: 25 L2: 48 +DEAL:: Number of degrees of freedom: 57 (by level: 8, 25, 48) +DEAL:cg::Starting value 0.4679 +DEAL:cg::Convergence step 7 value 0 +DEAL:: 7 CG iterations needed to obtain convergence. +DEAL::Cycle 2: +DEAL:: Number of active cells: 92 +DEAL::Number of degrees of freedom: 117 L0: 8 L1: 25 L2: 80 L3: 60 +DEAL:: Number of degrees of freedom: 117 (by level: 8, 25, 80, 60) +DEAL:cg::Starting value 0.3390 +DEAL:cg::Convergence step 9 value 0 +DEAL:: 9 CG iterations needed to obtain convergence. +DEAL::Cycle 3: +DEAL:: Number of active cells: 188 +DEAL::Number of degrees of freedom: 221 L0: 8 L1: 25 L2: 80 L3: 200 +DEAL:: Number of degrees of freedom: 221 (by level: 8, 25, 80, 200) +DEAL:cg::Starting value 0.2689 +DEAL:cg::Convergence step 12 value 0 +DEAL:: 12 CG iterations needed to obtain convergence. +DEAL::Cycle 4: +DEAL:: Number of active cells: 416 +DEAL::Number of degrees of freedom: 485 L0: 8 L1: 25 L2: 89 L3: 288 L4: 280 +DEAL:: Number of degrees of freedom: 485 (by level: 8, 25, 89, 288, 280) +DEAL:cg::Starting value 0.1841 +DEAL:cg::Convergence step 12 value 0 +DEAL:: 12 CG iterations needed to obtain convergence. +DEAL::Cycle 5: +DEAL:: Number of active cells: 800 +DEAL::Number of degrees of freedom: 925 L0: 8 L1: 25 L2: 89 L3: 288 L4: 784 L5: 132 +DEAL:: Number of degrees of freedom: 925 (by level: 8, 25, 89, 288, 784, 132) +DEAL:cg::Starting value 0.1440 +DEAL:cg::Convergence step 14 value 0 +DEAL:: 14 CG iterations needed to obtain convergence. +DEAL::Cycle 6: +DEAL:: Number of active cells: 1628 +DEAL::Number of degrees of freedom: 1865 L0: 8 L1: 25 L2: 89 L3: 304 L4: 1000 L5: 1164 L6: 72 +DEAL:: Number of degrees of freedom: 1865 (by level: 8, 25, 89, 304, 1000, 1164, 72) +DEAL:cg::Starting value 0.1174 +DEAL:cg::Convergence step 15 value 0 +DEAL:: 15 CG iterations needed to obtain convergence. +DEAL::Cycle 7: +DEAL:: Number of active cells: 3194 +DEAL::Number of degrees of freedom: 3603 L0: 8 L1: 25 L2: 89 L3: 328 L4: 1032 L5: 2200 L6: 1392 +DEAL:: Number of degrees of freedom: 3603 (by level: 8, 25, 89, 328, 1032, 2200, 1392) +DEAL:cg::Starting value 0.09098 +DEAL:cg::Convergence step 16 value 0 +DEAL:: 16 CG iterations needed to obtain convergence.