From 8eb95f84be6d2925111a53d7517797f15a34068b Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Tue, 11 Aug 2020 16:47:31 +0200 Subject: [PATCH] New test case --- tests/multigrid/step-16-08.cc | 650 ++++++++++++++++++++++++++++++ tests/multigrid/step-16-08.output | 57 +++ 2 files changed, 707 insertions(+) create mode 100644 tests/multigrid/step-16-08.cc create mode 100644 tests/multigrid/step-16-08.output diff --git a/tests/multigrid/step-16-08.cc b/tests/multigrid/step-16-08.cc new file mode 100644 index 0000000000..807f7e8add --- /dev/null +++ b/tests/multigrid/step-16-08.cc @@ -0,0 +1,650 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2003 - 2020 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +// Multigrid for continuous finite elements without MeshWorker. Similar to the +// step-16 test but directly applying the constraints + +#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 "../tests.h" + + +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; + + AffineConstraints constraints; + + Vector solution; + Vector system_rhs; + + const unsigned int degree; + + MGLevelObject mg_sparsity_patterns; + MGLevelObject> mg_matrices; + MGLevelObject mg_interface_sparsity_patterns; + 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 < n_points; ++i) + values[i] = Coefficient::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(); + deallog << "Number of degrees of freedom: " << mg_dof_handler.n_dofs(); + + for (unsigned int l = 0; l < triangulation.n_levels(); ++l) + deallog << " " << 'L' << l << ": " << mg_dof_handler.n_dofs(l); + deallog << std::endl; + + 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(); + DoFTools::make_hanging_node_constraints(mg_dof_handler, constraints); + std::map *> dirichlet_boundary; + Functions::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); + mg_constrained_dofs.make_zero_boundary_constraints(mg_dof_handler, {0}); +} + + +template +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_point < n_q_points; ++q_point) + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + for (unsigned int j = 0; j < dofs_per_cell; ++j) + cell_matrix(i, j) += + (coefficient_values[q_point] * + fe_values.shape_grad(i, q_point) * + fe_values.shape_grad(j, q_point) * fe_values.JxW(q_point)); + + cell_rhs(i) += (fe_values.shape_value(i, q_point) * 1.0 * + fe_values.JxW(q_point)); + } + + cell->get_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() +{ + 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); + mg_interface_sparsity_patterns.resize(0, n_levels - 1); + + std::vector> boundary_constraints( + triangulation.n_levels()); + AffineConstraints empty_constraints; + for (unsigned int level = 0; level < n_levels; ++level) + { + boundary_constraints[level].add_lines( + mg_constrained_dofs.get_refinement_edge_indices(level)); + boundary_constraints[level].add_lines( + mg_constrained_dofs.get_boundary_indices(level)); + boundary_constraints[level].close(); + + DynamicSparsityPattern csp; + csp.reinit(mg_dof_handler.n_dofs(level), mg_dof_handler.n_dofs(level)); + MGTools::make_sparsity_pattern( + mg_dof_handler, csp, level, boundary_constraints[level], false); + mg_sparsity_patterns[level].copy_from(csp); + + csp.reinit(mg_dof_handler.n_dofs(level), mg_dof_handler.n_dofs(level)); + MGTools::make_sparsity_pattern( + mg_dof_handler, csp, level, empty_constraints, true); + mg_interface_sparsity_patterns[level].copy_from(csp); + + mg_matrices[level].reinit(mg_sparsity_patterns[level]); + mg_interface_matrices[level].reinit( + mg_interface_sparsity_patterns[level]); + } + + 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); + + typename DoFHandler::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_point < n_q_points; ++q_point) + for (unsigned int i = 0; i < dofs_per_cell; ++i) + for (unsigned int j = 0; j < dofs_per_cell; ++j) + cell_matrix(i, j) += + (coefficient_values[q_point] * fe_values.shape_grad(i, q_point) * + fe_values.shape_grad(j, q_point) * fe_values.JxW(q_point)); + + cell->get_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; i < dofs_per_cell; ++i) + for (unsigned int j = 0; j < dofs_per_cell; ++j) + if (mg_constrained_dofs.at_refinement_edge(lvl, + local_dof_indices[i]) && + !mg_constrained_dofs.at_refinement_edge(lvl, + local_dof_indices[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; + std::cout << i << " " << j << "\n"; + } + + + empty_constraints.distribute_local_to_global( + cell_matrix, local_dof_indices, mg_interface_matrices[cell->level()]); + } +} + + + +// @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(mg_dof_handler); + + FullMatrix coarse_matrix; + coarse_matrix.copy_from(mg_matrices[0]); + MGCoarseGridHouseholder<> coarse_grid_solver; + coarse_grid_solver.initialize(coarse_matrix); + + typedef PreconditionSOR> Smoother; + GrowingVectorMemory<> vector_memory; + MGSmootherRelaxation, Smoother, Vector> + mg_smoother; + mg_smoother.initialize(mg_matrices); + 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_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( + mg_dof_handler, + QGauss(3), + std::map *>(), + 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 SphericalManifold boundary; + triangulation.set_manifold(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 < triangulation.n_levels(); ++level) + deallog << mg_dof_handler.n_dofs(level) + << (level == triangulation.n_levels() - 1 ? ")" : ", "); + deallog << std::endl; + + assemble_system(); + assemble_multigrid(); + + solve(); + output_results(cycle); + } +} + + +// @sect3{The main() function} +// +// This is again the same function as +// in step-6: +int +main() +{ + initlog(); + deallog << std::setprecision(4); + + try + { + LaplaceProblem<2> 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-08.output b/tests/multigrid/step-16-08.output new file mode 100644 index 0000000000..b84e79a293 --- /dev/null +++ b/tests/multigrid/step-16-08.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 7 value 0 +DEAL:: 7 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 8 value 0 +DEAL:: 8 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: 392 +DEAL::Number of degrees of freedom: 453 L0: 8 L1: 25 L2: 89 L3: 288 L4: 240 +DEAL:: Number of degrees of freedom: 453 (by level: 8, 25, 89, 288, 240) +DEAL:cg::Starting value 0.1854 +DEAL:cg::Convergence step 13 value 0 +DEAL:: 13 CG iterations needed to obtain convergence. +DEAL::Cycle 5: +DEAL:: Number of active cells: 752 +DEAL::Number of degrees of freedom: 865 L0: 8 L1: 25 L2: 89 L3: 288 L4: 760 L5: 60 +DEAL:: Number of degrees of freedom: 865 (by level: 8, 25, 89, 288, 760, 60) +DEAL:cg::Starting value 0.1471 +DEAL:cg::Convergence step 14 value 0 +DEAL:: 14 CG iterations needed to obtain convergence. +DEAL::Cycle 6: +DEAL:: Number of active cells: 1532 +DEAL::Number of degrees of freedom: 1741 L0: 8 L1: 25 L2: 89 L3: 304 L4: 1000 L5: 984 L6: 72 +DEAL:: Number of degrees of freedom: 1741 (by level: 8, 25, 89, 304, 1000, 984, 72) +DEAL:cg::Starting value 0.1188 +DEAL:cg::Convergence step 14 value 0 +DEAL:: 14 CG iterations needed to obtain convergence. +DEAL::Cycle 7: +DEAL:: Number of active cells: 3020 +DEAL::Number of degrees of freedom: 3417 L0: 8 L1: 25 L2: 89 L3: 328 L4: 1032 L5: 1976 L6: 1360 +DEAL:: Number of degrees of freedom: 3417 (by level: 8, 25, 89, 328, 1032, 1976, 1360) +DEAL:cg::Starting value 0.09316 +DEAL:cg::Convergence step 19 value 0 +DEAL:: 19 CG iterations needed to obtain convergence. -- 2.39.5