From f0f137b6a03eaa69a98df9eea0a89970123075e8 Mon Sep 17 00:00:00 2001 From: "Thomas C. Clevenger" Date: Wed, 23 Oct 2019 12:14:09 +0200 Subject: [PATCH] step-50 cleanup fix several issues listed in #8482 --- examples/step-50/step-50.cc | 220 +++++++++++++++++------------------- 1 file changed, 102 insertions(+), 118 deletions(-) diff --git a/examples/step-50/step-50.cc b/examples/step-50/step-50.cc index 4618477061..0c7218a698 100644 --- a/examples/step-50/step-50.cc +++ b/examples/step-50/step-50.cc @@ -107,6 +107,8 @@ namespace Step50 public: LaplaceProblem(const unsigned int deg); void run(); + using MatrixType = LA::MPI::SparseMatrix; + using VectorType = LA::MPI::Vector; private: void setup_system(); @@ -116,31 +118,30 @@ namespace Step50 void refine_grid(); void output_results(const unsigned int cycle) const; + MPI_Comm mpi_communicator; ConditionalOStream pcout; parallel::distributed::Triangulation triangulation; FE_Q fe; - DoFHandler mg_dof_handler; + DoFHandler dof_handler; - using matrix_t = LA::MPI::SparseMatrix; - using vector_t = LA::MPI::Vector; - matrix_t system_matrix; + MatrixType system_matrix; IndexSet locally_relevant_set; AffineConstraints constraints; - vector_t solution; - vector_t system_rhs; + VectorType solution; + VectorType system_rhs; const unsigned int degree; // Finally we are storing the various parallel multigrid matrices. Our // problem is self-adjoint, so the interface matrices are the transpose // of each other, so we only need to compute/store them once. - MGLevelObject mg_matrices; - MGLevelObject mg_interface_matrices; + MGLevelObject mg_matrices; + MGLevelObject mg_interface_matrices; // MGConstrainedDoFs mg_constrained_dofs; }; @@ -225,13 +226,15 @@ namespace Step50 // triangulation class. template LaplaceProblem::LaplaceProblem(const unsigned int degree) - : pcout(std::cout, (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)) - , triangulation(MPI_COMM_WORLD, + : mpi_communicator(MPI_COMM_WORLD) + , pcout(std::cout, + (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)) + , triangulation(mpi_communicator, Triangulation::limit_level_difference_at_vertices, parallel::distributed::Triangulation< dim>::construct_multigrid_hierarchy) , fe(degree) - , mg_dof_handler(triangulation) + , dof_handler(triangulation) , degree(degree) {} @@ -245,14 +248,14 @@ namespace Step50 template void LaplaceProblem::setup_system() { - mg_dof_handler.distribute_dofs(fe); - mg_dof_handler.distribute_mg_dofs(); + dof_handler.clear(); + dof_handler.distribute_dofs(fe); + dof_handler.distribute_mg_dofs(); - DoFTools::extract_locally_relevant_dofs(mg_dof_handler, - locally_relevant_set); + DoFTools::extract_locally_relevant_dofs(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); + solution.reinit(dof_handler.locally_owned_dofs(), mpi_communicator); + system_rhs.reinit(dof_handler.locally_owned_dofs(), mpi_communicator); // But it starts to be a wee bit different // here, although this still doesn't have @@ -274,25 +277,26 @@ namespace Step50 // correctly into the global linear system // right away, without the need for a later // clean-up stage: + constraints.clear(); constraints.reinit(locally_relevant_set); - DoFTools::make_hanging_node_constraints(mg_dof_handler, constraints); - - std::set dirichlet_boundary_ids; - std::map *> dirichlet_boundary; - Functions::ConstantFunction homogeneous_dirichlet_bc(1.0); - dirichlet_boundary_ids.insert(0); - dirichlet_boundary[0] = &homogeneous_dirichlet_bc; - VectorTools::interpolate_boundary_values(mg_dof_handler, - dirichlet_boundary, + DoFTools::make_hanging_node_constraints(dof_handler, constraints); + + const types::boundary_id boundary = 0; + std::set dirichlet_boundary_ids{boundary}; + Functions::ConstantFunction boundary_values(1.0); + + VectorTools::interpolate_boundary_values(dof_handler, + boundary, + boundary_values, constraints); 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(), + system_matrix.clear(); + DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs()); + DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints); + system_matrix.reinit(dof_handler.locally_owned_dofs(), dsp, - MPI_COMM_WORLD, + mpi_communicator, true); @@ -302,8 +306,8 @@ namespace Step50 // pass the dirichlet_boundary // here as well. mg_constrained_dofs.clear(); - mg_constrained_dofs.initialize(mg_dof_handler); - mg_constrained_dofs.make_zero_boundary_constraints(mg_dof_handler, + mg_constrained_dofs.initialize(dof_handler); + mg_constrained_dofs.make_zero_boundary_constraints(dof_handler, dirichlet_boundary_ids); @@ -357,22 +361,33 @@ namespace Step50 // matrices. for (unsigned int level = 0; level < n_levels; ++level) { - DynamicSparsityPattern dsp(mg_dof_handler.n_dofs(level), - mg_dof_handler.n_dofs(level)); - MGTools::make_sparsity_pattern(mg_dof_handler, dsp, level); - - mg_matrices[level].reinit(mg_dof_handler.locally_owned_mg_dofs(level), - mg_dof_handler.locally_owned_mg_dofs(level), - dsp, - MPI_COMM_WORLD, - true); - - mg_interface_matrices[level].reinit( - mg_dof_handler.locally_owned_mg_dofs(level), - mg_dof_handler.locally_owned_mg_dofs(level), - dsp, - MPI_COMM_WORLD, - true); + { + DynamicSparsityPattern dsp(dof_handler.n_dofs(level), + dof_handler.n_dofs(level)); + MGTools::make_sparsity_pattern(dof_handler, dsp, level); + + mg_matrices[level].reinit(dof_handler.locally_owned_mg_dofs(level), + dof_handler.locally_owned_mg_dofs(level), + dsp, + mpi_communicator, + true); + } + + { + DynamicSparsityPattern dsp(dof_handler.n_dofs(level), + dof_handler.n_dofs(level)); + MGTools::make_interface_sparsity_pattern(dof_handler, + mg_constrained_dofs, + dsp, + level); + + mg_interface_matrices[level].reinit( + dof_handler.locally_owned_mg_dofs(level), + dof_handler.locally_owned_mg_dofs(level), + dsp, + mpi_communicator, + true); + } } } @@ -416,7 +431,7 @@ namespace Step50 const Coefficient coefficient; std::vector coefficient_values(n_q_points); - for (const auto &cell : mg_dof_handler.active_cell_iterators()) + for (const auto &cell : dof_handler.active_cell_iterators()) if (cell->is_locally_owned()) { cell_matrix = 0; @@ -526,12 +541,11 @@ namespace Step50 // AffineConstraints::add_lines(): std::vector> boundary_constraints( triangulation.n_global_levels()); - AffineConstraints empty_constraints; for (unsigned int level = 0; level < triangulation.n_global_levels(); ++level) { IndexSet dofset; - DoFTools::extract_locally_relevant_level_dofs(mg_dof_handler, + DoFTools::extract_locally_relevant_level_dofs(dof_handler, level, dofset); boundary_constraints[level].reinit(dofset); @@ -549,7 +563,7 @@ namespace Step50 // need a right hand side, and more significantly (ii) we don't // just loop over all active cells, but in fact all cells, active // or not. - for (const auto &cell : mg_dof_handler.cell_iterators()) + for (const auto &cell : dof_handler.cell_iterators()) if (cell->level_subdomain_id() == triangulation.locally_owned_subdomain()) { cell_matrix = 0; @@ -619,41 +633,13 @@ namespace Step50 // the solve() function) be able to just pass // the transpose matrix where necessary. - 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; i < dofs_per_cell; ++i) for (unsigned int j = 0; j < dofs_per_cell; ++j) - if (interface_dofs_on_level.is_element( - local_dof_indices[i]) // at_refinement_edge(i) - && !interface_dofs_on_level.is_element( - local_dof_indices[j]) // !at_refinement_edge(j) - && - ((!mg_constrained_dofs.is_boundary_index( - lvl, local_dof_indices[i]) && - !mg_constrained_dofs.is_boundary_index( - lvl, - local_dof_indices[j])) // ( !boundary(i) && !boundary(j) ) - || (mg_constrained_dofs.is_boundary_index( - lvl, local_dof_indices[i]) && - local_dof_indices[i] == - local_dof_indices[j]) // ( boundary(i) && boundary(j) - // && i==j ) - )) - { - // do nothing, so add entries to interface matrix - } - else - { - cell_matrix(i, j) = 0; - } - - - empty_constraints.distribute_local_to_global( - cell_matrix, - local_dof_indices, - mg_interface_matrices[cell->level()]); + if (mg_constrained_dofs.is_interface_matrix_entry( + cell->level(), local_dof_indices[i], local_dof_indices[j])) + mg_interface_matrices[cell->level()].add(local_dof_indices[i], + local_dof_indices[j], + cell_matrix(i, j)); } for (unsigned int i = 0; i < triangulation.n_global_levels(); ++i) @@ -693,20 +679,25 @@ namespace Step50 { // Create the object that deals with the transfer between // different refinement levels. - MGTransferPrebuilt mg_transfer(mg_constrained_dofs); + MGTransferPrebuilt mg_transfer(mg_constrained_dofs); // Now the prolongation matrix has to be built. - mg_transfer.build_matrices(mg_dof_handler); + mg_transfer.build_matrices(dof_handler); - matrix_t &coarse_matrix = mg_matrices[0]; + MatrixType &coarse_matrix = mg_matrices[0]; SolverControl coarse_solver_control(1000, 1e-10, false, false); - SolverCG coarse_solver(coarse_solver_control); - PreconditionIdentity id; - MGCoarseGridIterativeSolver, - matrix_t, + SolverCG coarse_solver(coarse_solver_control); + + PreconditionIdentity prec; + // TrilinosWrappers::PreconditionAMG prec; + // prec.initialize(coarse_matrix); + + MGCoarseGridIterativeSolver, + MatrixType, + // TrilinosWrappers::PreconditionAMG PreconditionIdentity> - coarse_grid_solver(coarse_solver, coarse_matrix, id); + coarse_grid_solver(coarse_solver, coarse_matrix, prec); // The next component of a multilevel solver or preconditioner is // that we need a smoother on each level. A common choice for this @@ -739,7 +730,7 @@ namespace Step50 // symmetric operator even for nonsymmetric // smoothers: using Smoother = LA::MPI::PreconditionJacobi; - MGSmootherPrecondition mg_smoother; + MGSmootherPrecondition mg_smoother; mg_smoother.initialize(mg_matrices, Smoother::AdditionalData(0.5)); mg_smoother.set_steps(2); // mg_smoother.set_symmetric(false); @@ -757,27 +748,27 @@ namespace Step50 // both up and down versions of the // operator with the matrices we already // built: - mg::Matrix mg_matrix(mg_matrices); - mg::Matrix mg_interface_up(mg_interface_matrices); - mg::Matrix mg_interface_down(mg_interface_matrices); + mg::Matrix mg_matrix(mg_matrices); + mg::Matrix mg_interface_up(mg_interface_matrices); + mg::Matrix mg_interface_down(mg_interface_matrices); // Now, we are ready to set up the // V-cycle operator and the // multilevel preconditioner. - Multigrid mg( + Multigrid mg( mg_matrix, coarse_grid_solver, mg_transfer, mg_smoother, mg_smoother); // mg.set_debug(6); mg.set_edge_matrices(mg_interface_down, mg_interface_up); - PreconditionMG> preconditioner( - mg_dof_handler, mg, mg_transfer); + PreconditionMG> + preconditioner(dof_handler, mg, mg_transfer); // With all this together, we can finally // get about solving the linear system in // the usual way: - SolverControl solver_control(500, 1e-8 * system_rhs.l2_norm(), false); - SolverCG solver(solver_control); + SolverControl solver_control(500, 1e-8 * system_rhs.l2_norm(), false); + SolverCG solver(solver_control); if (false) { @@ -826,17 +817,17 @@ namespace Step50 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.reinit(locally_relevant_set, mpi_communicator); temp_solution = solution; KellyErrorEstimator::estimate( - mg_dof_handler, + dof_handler, QGauss(degree + 1), std::map *>(), temp_solution, estimated_error_per_cell); - parallel::distributed::GridRefinement::refine_and_coarsen_fixed_fraction( + parallel::distributed::GridRefinement::refine_and_coarsen_fixed_number( triangulation, estimated_error_per_cell, 0.3, 0.0); triangulation.execute_coarsening_and_refinement(); @@ -850,18 +841,11 @@ namespace Step50 DataOut data_out; LA::MPI::Vector temp_solution; - temp_solution.reinit(locally_relevant_set, MPI_COMM_WORLD); + temp_solution.reinit(locally_relevant_set, mpi_communicator); temp_solution = solution; - - LA::MPI::Vector temp = solution; - system_matrix.residual(temp, solution, system_rhs); - LA::MPI::Vector res_ghosted = temp_solution; - res_ghosted = temp; - - data_out.attach_dof_handler(mg_dof_handler); + data_out.attach_dof_handler(dof_handler); data_out.add_data_vector(temp_solution, "solution"); - data_out.add_data_vector(res_ghosted, "res"); Vector subdomain(triangulation.n_active_cells()); for (unsigned int i = 0; i < subdomain.size(); ++i) subdomain(i) = triangulation.locally_owned_subdomain(); @@ -931,11 +915,11 @@ namespace Step50 setup_system(); - pcout << " Number of degrees of freedom: " << mg_dof_handler.n_dofs() + pcout << " Number of degrees of freedom: " << dof_handler.n_dofs() << " (by level: "; for (unsigned int level = 0; level < triangulation.n_global_levels(); ++level) - pcout << mg_dof_handler.n_dofs(level) + pcout << dof_handler.n_dofs(level) << (level == triangulation.n_global_levels() - 1 ? ")" : ", "); pcout << std::endl; -- 2.39.5