From df8572f9f29db0fa35455cd87274961ce90b5bbe Mon Sep 17 00:00:00 2001 From: Denis Davydov Date: Thu, 7 Apr 2016 09:49:05 +0200 Subject: [PATCH] add GMG unit tests which compare serial step-16 to parallel step-50 at the same conditions --- tests/multigrid/step-16-50-mpi.cc | 581 ++++++++++++++++++ ...linos=true.with_p4est=true.mpirun=1.output | 26 + ...linos=true.with_p4est=true.mpirun=3.output | 26 + ...linos=true.with_p4est=true.mpirun=7.output | 26 + tests/multigrid/step-16-50-serial.cc | 510 +++++++++++++++ tests/multigrid/step-16-50-serial.output | 26 + 6 files changed, 1195 insertions(+) create mode 100644 tests/multigrid/step-16-50-mpi.cc create mode 100644 tests/multigrid/step-16-50-mpi.with_mpi=true.with_trilinos=true.with_p4est=true.mpirun=1.output create mode 100644 tests/multigrid/step-16-50-mpi.with_mpi=true.with_trilinos=true.with_p4est=true.mpirun=3.output create mode 100644 tests/multigrid/step-16-50-mpi.with_mpi=true.with_trilinos=true.with_p4est=true.mpirun=7.output create mode 100644 tests/multigrid/step-16-50-serial.cc create mode 100644 tests/multigrid/step-16-50-serial.output diff --git a/tests/multigrid/step-16-50-mpi.cc b/tests/multigrid/step-16-50-mpi.cc new file mode 100644 index 0000000000..c1e2be4d2e --- /dev/null +++ b/tests/multigrid/step-16-50-mpi.cc @@ -0,0 +1,581 @@ +/* --------------------------------------------------------------------- + * + * Copyright (C) 2003 - 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. + * + * --------------------------------------------------------------------- + */ + +// compare serial version of step-16 to parallel version of step-50 +// using the same conditions: smoothers, solvers, marking strategy, +// RHS, coefficient, domain, etc + +#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 +#include +#include +#include +#include + + +#include + + +namespace LA +{ +#ifdef USE_PETSC_LA + using namespace dealii::LinearAlgebraPETSc; +#else + using namespace dealii::LinearAlgebraTrilinos; +#endif +} + +#include +#include +#include + +namespace Step50 +{ + 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 (); + + parallel::distributed::Triangulation triangulation; + FE_Q fe; + DoFHandler mg_dof_handler; + + typedef LA::MPI::SparseMatrix matrix_t; + typedef LA::MPI::Vector vector_t; + + matrix_t system_matrix; + + IndexSet locally_relevant_set; + + ConstraintMatrix hanging_node_constraints; + ConstraintMatrix constraints; + + vector_t solution; + vector_t system_rhs; + + const unsigned int degree; + + 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 5; + 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 (MPI_COMM_WORLD,Triangulation:: + limit_level_difference_at_vertices, + parallel::distributed::Triangulation::construct_multigrid_hierarchy), + 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); + + DoFTools::extract_locally_relevant_dofs (mg_dof_handler, + locally_relevant_set); + + solution.reinit(mg_dof_handler.locally_owned_dofs(), MPI_COMM_WORLD); + system_rhs.reinit(mg_dof_handler.locally_owned_dofs(), MPI_COMM_WORLD); + + constraints.reinit (locally_relevant_set); + hanging_node_constraints.reinit (locally_relevant_set); + DoFTools::make_hanging_node_constraints (mg_dof_handler, hanging_node_constraints); + DoFTools::make_hanging_node_constraints (mg_dof_handler, constraints); + + typename FunctionMap::type dirichlet_boundary; + ZeroFunction homogeneous_dirichlet_bc; + dirichlet_boundary[0] = &homogeneous_dirichlet_bc; + VectorTools::interpolate_boundary_values (mg_dof_handler, + dirichlet_boundary, + constraints); + constraints.close (); + hanging_node_constraints.close (); + + DynamicSparsityPattern dsp(mg_dof_handler.n_dofs(), mg_dof_handler.n_dofs()); + DoFTools::make_sparsity_pattern (mg_dof_handler, dsp, constraints); + system_matrix.reinit (mg_dof_handler.locally_owned_dofs(), dsp, MPI_COMM_WORLD, true); + + + mg_constrained_dofs.clear(); + mg_constrained_dofs.initialize(mg_dof_handler, dirichlet_boundary); + + + const unsigned int n_levels = triangulation.n_global_levels(); + + mg_interface_matrices.resize(0, n_levels-1); + mg_interface_matrices.clear (); + mg_matrices.resize(0, n_levels-1); + mg_matrices.clear (); + + 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) + if (cell->is_locally_owned()) + { + 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); + } + + system_matrix.compress(VectorOperation::add); + system_rhs.compress(VectorOperation::add); + } + + + + 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_global_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) + if (cell->level_subdomain_id()==triangulation.locally_owned_subdomain()) + { + 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()]); + + + const IndexSet &interface_dofs_on_level + = mg_constrained_dofs.get_refinement_edge_indices(cell->level()); + const unsigned int lvl = cell->level(); + + for (unsigned int i=0; ilevel()]); + } + + for (unsigned int i=0; i + void LaplaceProblem::solve () + { + MGTransferPrebuilt mg_transfer(hanging_node_constraints, mg_constrained_dofs); + mg_transfer.build_matrices(mg_dof_handler); + + matrix_t &coarse_matrix = mg_matrices[0]; + + SolverControl coarse_solver_control (1000, 1e-10, false, false); + SolverCG coarse_solver(coarse_solver_control); + PreconditionIdentity id; + MGCoarseGridLACIteration,vector_t> coarse_grid_solver(coarse_solver, + coarse_matrix, + id); + + typedef LA::MPI::PreconditionJacobi Smoother; + MGSmootherPrecondition mg_smoother; + mg_smoother.initialize(mg_matrices, Smoother::AdditionalData(0.5)); + mg_smoother.set_steps(2); + + 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 > + preconditioner(mg_dof_handler, mg, mg_transfer); + + + SolverControl solver_control (500, 1e-8*system_rhs.l2_norm(), false); + SolverCG solver (solver_control); + + solution = 0; + solver.solve (system_matrix, solution, system_rhs, + preconditioner); + constraints.distribute (solution); + } + + + + + template + void LaplaceProblem::refine_grid () + { + Vector estimated_error_per_cell (triangulation.n_active_cells()); + + LA::MPI::Vector temp_solution; + temp_solution.reinit(locally_relevant_set, MPI_COMM_WORLD); + temp_solution = solution; + + KellyErrorEstimator::estimate (static_cast&>(mg_dof_handler), + QGauss(degree+1), + typename FunctionMap::type(), + temp_solution, + estimated_error_per_cell); + + const double threshold = 0.6 * Utilities::MPI::max(estimated_error_per_cell.linfty_norm(),MPI_COMM_WORLD); + GridRefinement::refine (triangulation, estimated_error_per_cell, threshold); + + for (typename Triangulation::active_cell_iterator + cell = triangulation.begin_active(); + cell != triangulation.end(); ++cell) + if (cell->subdomain_id() != triangulation.locally_owned_subdomain()) + { + cell->clear_refine_flag (); + cell->clear_coarsen_flag (); + } + + triangulation.execute_coarsening_and_refinement (); + } + + template + void LaplaceProblem::run () + { + for (unsigned int cycle=0; cycle<5; ++cycle) + { + deallog << "Cycle " << cycle << ':' << std::endl; + + if (cycle == 0) + { + GridGenerator::hyper_cube (triangulation); + + triangulation.refine_global (4); + } + else + refine_grid (); + + deallog << " Number of active cells: " + << triangulation.n_global_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/*degree*/); + 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; + throw; + } + catch (...) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + throw; + } + + return 0; +} diff --git a/tests/multigrid/step-16-50-mpi.with_mpi=true.with_trilinos=true.with_p4est=true.mpirun=1.output b/tests/multigrid/step-16-50-mpi.with_mpi=true.with_trilinos=true.with_p4est=true.mpirun=1.output new file mode 100644 index 0000000000..3f3684cfb0 --- /dev/null +++ b/tests/multigrid/step-16-50-mpi.with_mpi=true.with_trilinos=true.with_p4est=true.mpirun=1.output @@ -0,0 +1,26 @@ + +DEAL::Cycle 0: +DEAL:: Number of active cells: 256 +DEAL:: Number of degrees of freedom: 289 (by level: 4, 9, 25, 81, 289) +DEAL:cg::Starting value 0.585938 +DEAL:cg::Convergence step 7 value 4.58200e-09 +DEAL::Cycle 1: +DEAL:: Number of active cells: 418 +DEAL:: Number of degrees of freedom: 495 (by level: 4, 9, 25, 81, 289, 305) +DEAL:cg::Starting value 0.533937 +DEAL:cg::Convergence step 10 value 7.40543e-10 +DEAL::Cycle 2: +DEAL:: Number of active cells: 865 +DEAL:: Number of degrees of freedom: 952 (by level: 4, 9, 25, 81, 289, 769, 153) +DEAL:cg::Starting value 0.378780 +DEAL:cg::Convergence step 11 value 1.39894e-09 +DEAL::Cycle 3: +DEAL:: Number of active cells: 1219 +DEAL:: Number of degrees of freedom: 1354 (by level: 4, 9, 25, 81, 289, 953, 277, 261) +DEAL:cg::Starting value 0.343729 +DEAL:cg::Convergence step 13 value 6.16552e-10 +DEAL::Cycle 4: +DEAL:: Number of active cells: 1297 +DEAL:: Number of degrees of freedom: 1446 (by level: 4, 9, 25, 81, 289, 961, 313, 349) +DEAL:cg::Starting value 0.339440 +DEAL:cg::Convergence step 13 value 7.00850e-10 diff --git a/tests/multigrid/step-16-50-mpi.with_mpi=true.with_trilinos=true.with_p4est=true.mpirun=3.output b/tests/multigrid/step-16-50-mpi.with_mpi=true.with_trilinos=true.with_p4est=true.mpirun=3.output new file mode 100644 index 0000000000..3f3684cfb0 --- /dev/null +++ b/tests/multigrid/step-16-50-mpi.with_mpi=true.with_trilinos=true.with_p4est=true.mpirun=3.output @@ -0,0 +1,26 @@ + +DEAL::Cycle 0: +DEAL:: Number of active cells: 256 +DEAL:: Number of degrees of freedom: 289 (by level: 4, 9, 25, 81, 289) +DEAL:cg::Starting value 0.585938 +DEAL:cg::Convergence step 7 value 4.58200e-09 +DEAL::Cycle 1: +DEAL:: Number of active cells: 418 +DEAL:: Number of degrees of freedom: 495 (by level: 4, 9, 25, 81, 289, 305) +DEAL:cg::Starting value 0.533937 +DEAL:cg::Convergence step 10 value 7.40543e-10 +DEAL::Cycle 2: +DEAL:: Number of active cells: 865 +DEAL:: Number of degrees of freedom: 952 (by level: 4, 9, 25, 81, 289, 769, 153) +DEAL:cg::Starting value 0.378780 +DEAL:cg::Convergence step 11 value 1.39894e-09 +DEAL::Cycle 3: +DEAL:: Number of active cells: 1219 +DEAL:: Number of degrees of freedom: 1354 (by level: 4, 9, 25, 81, 289, 953, 277, 261) +DEAL:cg::Starting value 0.343729 +DEAL:cg::Convergence step 13 value 6.16552e-10 +DEAL::Cycle 4: +DEAL:: Number of active cells: 1297 +DEAL:: Number of degrees of freedom: 1446 (by level: 4, 9, 25, 81, 289, 961, 313, 349) +DEAL:cg::Starting value 0.339440 +DEAL:cg::Convergence step 13 value 7.00850e-10 diff --git a/tests/multigrid/step-16-50-mpi.with_mpi=true.with_trilinos=true.with_p4est=true.mpirun=7.output b/tests/multigrid/step-16-50-mpi.with_mpi=true.with_trilinos=true.with_p4est=true.mpirun=7.output new file mode 100644 index 0000000000..3f3684cfb0 --- /dev/null +++ b/tests/multigrid/step-16-50-mpi.with_mpi=true.with_trilinos=true.with_p4est=true.mpirun=7.output @@ -0,0 +1,26 @@ + +DEAL::Cycle 0: +DEAL:: Number of active cells: 256 +DEAL:: Number of degrees of freedom: 289 (by level: 4, 9, 25, 81, 289) +DEAL:cg::Starting value 0.585938 +DEAL:cg::Convergence step 7 value 4.58200e-09 +DEAL::Cycle 1: +DEAL:: Number of active cells: 418 +DEAL:: Number of degrees of freedom: 495 (by level: 4, 9, 25, 81, 289, 305) +DEAL:cg::Starting value 0.533937 +DEAL:cg::Convergence step 10 value 7.40543e-10 +DEAL::Cycle 2: +DEAL:: Number of active cells: 865 +DEAL:: Number of degrees of freedom: 952 (by level: 4, 9, 25, 81, 289, 769, 153) +DEAL:cg::Starting value 0.378780 +DEAL:cg::Convergence step 11 value 1.39894e-09 +DEAL::Cycle 3: +DEAL:: Number of active cells: 1219 +DEAL:: Number of degrees of freedom: 1354 (by level: 4, 9, 25, 81, 289, 953, 277, 261) +DEAL:cg::Starting value 0.343729 +DEAL:cg::Convergence step 13 value 6.16552e-10 +DEAL::Cycle 4: +DEAL:: Number of active cells: 1297 +DEAL:: Number of degrees of freedom: 1446 (by level: 4, 9, 25, 81, 289, 961, 313, 349) +DEAL:cg::Starting value 0.339440 +DEAL:cg::Convergence step 13 value 7.00850e-10 diff --git a/tests/multigrid/step-16-50-serial.cc b/tests/multigrid/step-16-50-serial.cc new file mode 100644 index 0000000000..51091884f2 --- /dev/null +++ b/tests/multigrid/step-16-50-serial.cc @@ -0,0 +1,510 @@ +/* --------------------------------------------------------------------- + * + * Copyright (C) 2003 - 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. + * + * --------------------------------------------------------------------- + */ + +// compare serial version of step-16 to parallel version of step-50 +// using the same conditions: smoothers, solvers, marking strategy, +// RHS, coefficient, domain, etc + +#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 hanging_node_constraints; + 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 5; + 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); + + sparsity_pattern.reinit (mg_dof_handler.n_dofs(), + mg_dof_handler.n_dofs(), + mg_dof_handler.max_couplings_between_dofs()); + DoFTools::make_sparsity_pattern ( + static_cast&>(mg_dof_handler), + sparsity_pattern); + + solution.reinit (mg_dof_handler.n_dofs()); + system_rhs.reinit (mg_dof_handler.n_dofs()); + + constraints.clear (); + hanging_node_constraints.clear (); + DoFTools::make_hanging_node_constraints (mg_dof_handler, constraints); + DoFTools::make_hanging_node_constraints (mg_dof_handler, hanging_node_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 (); + hanging_node_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 (); + mg_matrices.resize(0, n_levels-1); + mg_matrices.clear (); + 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()]); + + const unsigned int lvl = cell->level(); + + for (unsigned int i=0; ilevel()]); + } +} + + + + +template +void LaplaceProblem::solve () +{ + typedef SparseMatrix matrix_t; + typedef Vector vector_t; + + MGTransferPrebuilt mg_transfer(hanging_node_constraints, mg_constrained_dofs); + mg_transfer.build_matrices(mg_dof_handler); + + matrix_t &coarse_matrix = mg_matrices[0]; + + SolverControl coarse_solver_control (1000, 1e-10, false, false); + SolverCG coarse_solver(coarse_solver_control); + PreconditionIdentity id; + MGCoarseGridLACIteration,vector_t> coarse_grid_solver(coarse_solver, + coarse_matrix, + id); + + typedef PreconditionJacobi Smoother; + MGSmootherPrecondition mg_smoother; + mg_smoother.initialize(mg_matrices, Smoother::AdditionalData(0.5)); + mg_smoother.set_steps(2); + + 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 > + preconditioner(mg_dof_handler, mg, mg_transfer); + + SolverControl solver_control (500, 1e-8*system_rhs.l2_norm(), false); + SolverCG solver (solver_control); + + solution = 0; + solver.solve (system_matrix, solution, system_rhs, + preconditioner); + constraints.distribute (solution); +} + + + + +template +void LaplaceProblem::refine_grid () +{ + Vector estimated_error_per_cell (triangulation.n_active_cells()); + + KellyErrorEstimator::estimate (static_cast&>(mg_dof_handler), + QGauss(degree+1), + typename FunctionMap::type(), + solution, + estimated_error_per_cell); + + const double threshold = 0.6 * estimated_error_per_cell.linfty_norm(); + GridRefinement::refine (triangulation, estimated_error_per_cell, threshold); + + triangulation.execute_coarsening_and_refinement (); +} + + +template +void LaplaceProblem::run () +{ + for (unsigned int cycle=0; cycle<5; ++cycle) + { + deallog << "Cycle " << cycle << ':' << std::endl; + + if (cycle == 0) + { + GridGenerator::hyper_cube (triangulation); + triangulation.refine_global (4); + } + 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-50-serial.output b/tests/multigrid/step-16-50-serial.output new file mode 100644 index 0000000000..3f3684cfb0 --- /dev/null +++ b/tests/multigrid/step-16-50-serial.output @@ -0,0 +1,26 @@ + +DEAL::Cycle 0: +DEAL:: Number of active cells: 256 +DEAL:: Number of degrees of freedom: 289 (by level: 4, 9, 25, 81, 289) +DEAL:cg::Starting value 0.585938 +DEAL:cg::Convergence step 7 value 4.58200e-09 +DEAL::Cycle 1: +DEAL:: Number of active cells: 418 +DEAL:: Number of degrees of freedom: 495 (by level: 4, 9, 25, 81, 289, 305) +DEAL:cg::Starting value 0.533937 +DEAL:cg::Convergence step 10 value 7.40543e-10 +DEAL::Cycle 2: +DEAL:: Number of active cells: 865 +DEAL:: Number of degrees of freedom: 952 (by level: 4, 9, 25, 81, 289, 769, 153) +DEAL:cg::Starting value 0.378780 +DEAL:cg::Convergence step 11 value 1.39894e-09 +DEAL::Cycle 3: +DEAL:: Number of active cells: 1219 +DEAL:: Number of degrees of freedom: 1354 (by level: 4, 9, 25, 81, 289, 953, 277, 261) +DEAL:cg::Starting value 0.343729 +DEAL:cg::Convergence step 13 value 6.16552e-10 +DEAL::Cycle 4: +DEAL:: Number of active cells: 1297 +DEAL:: Number of degrees of freedom: 1446 (by level: 4, 9, 25, 81, 289, 961, 313, 349) +DEAL:cg::Starting value 0.339440 +DEAL:cg::Convergence step 13 value 7.00850e-10 -- 2.39.5