From 9ebdcd4464c8370e3d7dcb36ae6fc75997665c6b Mon Sep 17 00:00:00 2001 From: David Wells Date: Sat, 1 Sep 2018 13:40:25 -0400 Subject: [PATCH] Add a test based on step-37. The discussion of step-37 modernizations lead to the conclusion that each proposed inhomogeneous boundary enforcement strategy should correspond to a distinct test. This commit includes a test based on the first such enforcement strategy. --- tests/matrix_free/step-37-inhomogeneous-1.cc | 657 ++++++++++++++++++ ...ogeneous-1.mpirun=4.with_p4est=true.output | 150 ++++ 2 files changed, 807 insertions(+) create mode 100644 tests/matrix_free/step-37-inhomogeneous-1.cc create mode 100644 tests/matrix_free/step-37-inhomogeneous-1.mpirun=4.with_p4est=true.output diff --git a/tests/matrix_free/step-37-inhomogeneous-1.cc b/tests/matrix_free/step-37-inhomogeneous-1.cc new file mode 100644 index 0000000000..838c1900e8 --- /dev/null +++ b/tests/matrix_free/step-37-inhomogeneous-1.cc @@ -0,0 +1,657 @@ +/* --------------------------------------------------------------------- + * + * Copyright (C) 2009 - 2018 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. + * + * --------------------------------------------------------------------- + * + * Authors: Katharina Kormann, Martin Kronbichler, Uppsala University, + * 2009-2012, updated to MPI version with parallel vectors in 2016 + */ + +// This test verifies that the first strategy for enforcing inhomogeneous +// boundary conditions (i.e., using read_dof_values_plain to compute the +// contribution of the constraints to the right hand side in assemble_rhs) +// given in step-37 works. + + +#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" + +namespace Step37 +{ + using namespace dealii; + + const unsigned int degree_finite_element = 2; + + template + class ManufacturedSolution : public Function + { + public: + virtual double + value(const Point &p, + const unsigned int /*component*/ = 0) const override + { + const double pi = numbers::PI; + return std::sin(pi * p[0]) + std::cos(pi * p[1]); + } + }; + + + + template + class ManufacturedForcing : public Function + { + public: + virtual double + value(const Point &p, + const unsigned int /*compononent*/ = 0) const override + { + return value(p); + } + + template + number + value(const Point &p) const + { + const double pi = numbers::PI; + const number d = 40.0 * p.square() + 1.0; + return 20 * pi * pi * std::cos(pi * p[1]) / d + + 20 * pi * pi * std::sin(pi * p[0]) / d + + 1600 * pi * p[0] * std::cos(pi * p[0]) / (d * d) - + 1600 * pi * p[1] * std::sin(pi * p[1]) / (d * d); + } + }; + + + template + class Coefficient : public Function + { + public: + Coefficient() + : Function() + {} + + virtual double + value(const Point &p, const unsigned int component = 0) const override; + + template + number + value(const Point &p, const unsigned int component = 0) const; + }; + + + + template + template + number + Coefficient::value(const Point &p, + const unsigned int /*component*/) const + { + return 1. / (0.05 + 2. * p.square()); + } + + + + template + double + Coefficient::value(const Point & p, + const unsigned int component) const + { + return value(p, component); + } + + + template + class LaplaceOperator + : public MatrixFreeOperators:: + Base> + { + public: + using value_type = number; + + LaplaceOperator(); + + void + clear() override; + + void + evaluate_coefficient(const Coefficient &coefficient_function); + + virtual void + compute_diagonal() override; + + const Table<2, VectorizedArray> & + get_coefficient() const + { + return coefficient; + } + + private: + virtual void + apply_add( + LinearAlgebra::distributed::Vector & dst, + const LinearAlgebra::distributed::Vector &src) const override; + + void + local_apply(const MatrixFree & data, + LinearAlgebra::distributed::Vector & dst, + const LinearAlgebra::distributed::Vector &src, + const std::pair &cell_range) const; + + void + local_compute_diagonal( + const MatrixFree & data, + LinearAlgebra::distributed::Vector & dst, + const unsigned int & dummy, + const std::pair &cell_range) const; + + Table<2, VectorizedArray> coefficient; + }; + + + + template + LaplaceOperator::LaplaceOperator() + : MatrixFreeOperators::Base>() + {} + + + + template + void + LaplaceOperator::clear() + { + coefficient.reinit(0, 0); + MatrixFreeOperators::Base>:: + clear(); + } + + + + template + void + LaplaceOperator::evaluate_coefficient( + const Coefficient &coefficient_function) + { + const unsigned int n_cells = this->data->n_macro_cells(); + FEEvaluation phi(*this->data); + + coefficient.reinit(n_cells, phi.n_q_points); + for (unsigned int cell = 0; cell < n_cells; ++cell) + { + phi.reinit(cell); + for (unsigned int q = 0; q < phi.n_q_points; ++q) + coefficient(cell, q) = + coefficient_function.value(phi.quadrature_point(q)); + } + } + + + + template + void + LaplaceOperator::local_apply( + const MatrixFree & data, + LinearAlgebra::distributed::Vector & dst, + const LinearAlgebra::distributed::Vector &src, + const std::pair & cell_range) const + { + FEEvaluation phi(data); + + for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) + { + AssertDimension(coefficient.size(0), data.n_macro_cells()); + AssertDimension(coefficient.size(1), phi.n_q_points); + + phi.reinit(cell); + phi.read_dof_values(src); + phi.evaluate(false, true); + for (unsigned int q = 0; q < phi.n_q_points; ++q) + phi.submit_gradient(coefficient(cell, q) * phi.get_gradient(q), q); + phi.integrate(false, true); + phi.distribute_local_to_global(dst); + } + } + + + + template + void + LaplaceOperator::apply_add( + LinearAlgebra::distributed::Vector & dst, + const LinearAlgebra::distributed::Vector &src) const + { + this->data->cell_loop(&LaplaceOperator::local_apply, this, dst, src); + } + + + + template + void + LaplaceOperator::compute_diagonal() + { + this->inverse_diagonal_entries.reset( + new DiagonalMatrix>()); + LinearAlgebra::distributed::Vector &inverse_diagonal = + this->inverse_diagonal_entries->get_vector(); + this->data->initialize_dof_vector(inverse_diagonal); + unsigned int dummy = 0; + this->data->cell_loop(&LaplaceOperator::local_compute_diagonal, + this, + inverse_diagonal, + dummy); + + this->set_constrained_entries_to_one(inverse_diagonal); + + for (unsigned int i = 0; i < inverse_diagonal.local_size(); ++i) + { + Assert(inverse_diagonal.local_element(i) > 0., + ExcMessage("No diagonal entry in a positive definite operator " + "should be zero")); + inverse_diagonal.local_element(i) = + 1. / inverse_diagonal.local_element(i); + } + } + + + + template + void + LaplaceOperator::local_compute_diagonal( + const MatrixFree & data, + LinearAlgebra::distributed::Vector &dst, + const unsigned int &, + const std::pair &cell_range) const + { + FEEvaluation phi(data); + + AlignedVector> diagonal(phi.dofs_per_cell); + + for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) + { + AssertDimension(coefficient.size(0), data.n_macro_cells()); + AssertDimension(coefficient.size(1), phi.n_q_points); + + phi.reinit(cell); + for (unsigned int i = 0; i < phi.dofs_per_cell; ++i) + { + for (unsigned int j = 0; j < phi.dofs_per_cell; ++j) + phi.submit_dof_value(VectorizedArray(), j); + phi.submit_dof_value(make_vectorized_array(1.), i); + + phi.evaluate(false, true); + for (unsigned int q = 0; q < phi.n_q_points; ++q) + phi.submit_gradient(coefficient(cell, q) * phi.get_gradient(q), + q); + phi.integrate(false, true); + diagonal[i] = phi.get_dof_value(i); + } + for (unsigned int i = 0; i < phi.dofs_per_cell; ++i) + phi.submit_dof_value(diagonal[i], i); + phi.distribute_local_to_global(dst); + } + } + + + + template + class LaplaceProblem + { + public: + LaplaceProblem(); + void + run(); + + private: + void + setup_system(); + void + assemble_rhs(); + void + solve(); + void + output_results(const unsigned int cycle) const; + + parallel::distributed::Triangulation triangulation; + + FE_Q fe; + DoFHandler dof_handler; + + AffineConstraints constraints; + using SystemMatrixType = + LaplaceOperator; + SystemMatrixType system_matrix; + + MGConstrainedDoFs mg_constrained_dofs; + using LevelMatrixType = LaplaceOperator; + MGLevelObject mg_matrices; + + LinearAlgebra::distributed::Vector solution; + LinearAlgebra::distributed::Vector system_rhs; + + ConditionalOStream pcout; + }; + + + + template + LaplaceProblem::LaplaceProblem() + : triangulation(MPI_COMM_WORLD, + Triangulation::limit_level_difference_at_vertices, + parallel::distributed::Triangulation< + dim>::construct_multigrid_hierarchy) + , fe(degree_finite_element) + , dof_handler(triangulation) + , pcout(deallog.get_file_stream(), + Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0) + {} + + + + template + void + LaplaceProblem::setup_system() + { + system_matrix.clear(); + mg_matrices.clear_elements(); + + dof_handler.distribute_dofs(fe); + dof_handler.distribute_mg_dofs(); + + pcout << "Number of degrees of freedom: " << dof_handler.n_dofs() + << std::endl; + + IndexSet locally_relevant_dofs; + DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs); + + constraints.clear(); + constraints.reinit(locally_relevant_dofs); + DoFTools::make_hanging_node_constraints(dof_handler, constraints); + VectorTools::interpolate_boundary_values(dof_handler, + 0, + ManufacturedSolution(), + constraints); + constraints.close(); + + { + typename MatrixFree::AdditionalData additional_data; + additional_data.tasks_parallel_scheme = + MatrixFree::AdditionalData::none; + additional_data.mapping_update_flags = + (update_gradients | update_JxW_values | update_quadrature_points); + std::shared_ptr> system_mf_storage( + new MatrixFree()); + system_mf_storage->reinit(dof_handler, + constraints, + QGauss<1>(fe.degree + 1), + additional_data); + system_matrix.initialize(system_mf_storage); + } + + system_matrix.evaluate_coefficient(Coefficient()); + + system_matrix.initialize_dof_vector(solution); + system_matrix.initialize_dof_vector(system_rhs); + + const unsigned int nlevels = triangulation.n_global_levels(); + mg_matrices.resize(0, nlevels - 1); + + std::set dirichlet_boundary; + dirichlet_boundary.insert(0); + mg_constrained_dofs.initialize(dof_handler); + mg_constrained_dofs.make_zero_boundary_constraints(dof_handler, + dirichlet_boundary); + + for (unsigned int level = 0; level < nlevels; ++level) + { + IndexSet relevant_dofs; + DoFTools::extract_locally_relevant_level_dofs(dof_handler, + level, + relevant_dofs); + AffineConstraints level_constraints; + level_constraints.reinit(relevant_dofs); + level_constraints.add_lines( + mg_constrained_dofs.get_boundary_indices(level)); + level_constraints.close(); + + typename MatrixFree::AdditionalData additional_data; + additional_data.tasks_parallel_scheme = + MatrixFree::AdditionalData::none; + additional_data.mapping_update_flags = + (update_gradients | update_JxW_values | update_quadrature_points); + additional_data.level_mg_handler = level; + std::shared_ptr> mg_mf_storage_level( + new MatrixFree()); + mg_mf_storage_level->reinit(dof_handler, + level_constraints, + QGauss<1>(fe.degree + 1), + additional_data); + + mg_matrices[level].initialize(mg_mf_storage_level, + mg_constrained_dofs, + level); + mg_matrices[level].evaluate_coefficient(Coefficient()); + } + } + + + + template + void + LaplaceProblem::assemble_rhs() + { + solution = 0.0; + constraints.distribute(solution); + solution.update_ghost_values(); + system_rhs = 0.0; + + const Table<2, VectorizedArray> &coefficient = + system_matrix.get_coefficient(); + ManufacturedForcing forcing; + + FEEvaluation phi( + *system_matrix.get_matrix_free()); + for (unsigned int cell = 0; + cell < system_matrix.get_matrix_free()->n_macro_cells(); + ++cell) + { + phi.reinit(cell); + phi.read_dof_values_plain(solution); + phi.evaluate(false, true); + for (unsigned int q = 0; q < phi.n_q_points; ++q) + { + phi.submit_gradient(-coefficient(cell, q) * phi.get_gradient(q), q); + phi.submit_value(forcing.value(phi.quadrature_point(q)), q); + } + phi.integrate(true, true); + phi.distribute_local_to_global(system_rhs); + } + system_rhs.compress(VectorOperation::add); + } + + + + template + void + LaplaceProblem::solve() + { + MGTransferMatrixFree mg_transfer(mg_constrained_dofs); + mg_transfer.build(dof_handler); + + using SmootherType = + PreconditionChebyshev>; + mg::SmootherRelaxation> + mg_smoother; + MGLevelObject smoother_data; + smoother_data.resize(0, triangulation.n_global_levels() - 1); + for (unsigned int level = 0; level < triangulation.n_global_levels(); + ++level) + { + if (level > 0) + { + smoother_data[level].smoothing_range = 15.; + smoother_data[level].degree = 4; + smoother_data[level].eig_cg_n_iterations = 10; + } + else + { + smoother_data[0].smoothing_range = 1e-3; + smoother_data[0].degree = numbers::invalid_unsigned_int; + smoother_data[0].eig_cg_n_iterations = mg_matrices[0].m(); + } + mg_matrices[level].compute_diagonal(); + smoother_data[level].preconditioner = + mg_matrices[level].get_matrix_diagonal_inverse(); + } + mg_smoother.initialize(mg_matrices, smoother_data); + + MGCoarseGridApplySmoother> + mg_coarse; + mg_coarse.initialize(mg_smoother); + + mg::Matrix> mg_matrix( + mg_matrices); + + MGLevelObject> + mg_interface_matrices; + mg_interface_matrices.resize(0, triangulation.n_global_levels() - 1); + for (unsigned int level = 0; level < triangulation.n_global_levels(); + ++level) + mg_interface_matrices[level].initialize(mg_matrices[level]); + mg::Matrix> mg_interface( + mg_interface_matrices); + + Multigrid> mg( + mg_matrix, mg_coarse, mg_transfer, mg_smoother, mg_smoother); + mg.set_edge_matrices(mg_interface, mg_interface); + + PreconditionMG, + MGTransferMatrixFree> + preconditioner(dof_handler, mg, mg_transfer); + + SolverControl solver_control(100, 1e-12 * system_rhs.l2_norm()); + SolverCG> cg(solver_control); + + constraints.set_zero(solution); + cg.solve(system_matrix, solution, system_rhs, preconditioner); + + constraints.distribute(solution); + } + + + + template + void + LaplaceProblem::output_results(const unsigned int cycle) const + { + if (triangulation.n_global_active_cells() > 1000000) + return; + + Vector errors; + errors.reinit(triangulation.n_active_cells()); + solution.update_ghost_values(); + VectorTools::integrate_difference(dof_handler, + solution, + ManufacturedSolution(), + errors, + QIterated(QTrapez<1>(), 4), + VectorTools::NormType::Linfty_norm); + double max_cell_error = 1.0; + if (errors.begin() != errors.end()) + max_cell_error = *std::max_element(errors.begin(), errors.end()); + max_cell_error = Utilities::MPI::max(max_cell_error, MPI_COMM_WORLD); + Assert(max_cell_error != 0.0, ExcInternalError()); + pcout << "max error: " << max_cell_error << '\n'; + static double error = max_cell_error; + pcout << "error ratio: " + << Utilities::MPI::max(error, MPI_COMM_WORLD) / max_cell_error + << '\n'; + error = max_cell_error; + } + + + + template + void + LaplaceProblem::run() + { + for (unsigned int cycle = 0; cycle < 8 - dim; ++cycle) + { + pcout << "Cycle " << cycle << std::endl; + + if (cycle == 0) + { + GridGenerator::hyper_cube(triangulation, 0., 1.); + triangulation.refine_global(3 - dim); + } + triangulation.refine_global(1); + setup_system(); + assemble_rhs(); + solve(); + output_results(cycle); + pcout << std::endl; + } + } +} // namespace Step37 + + + +int +main(int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, 1); + MPILogInitAll mpi_log; + + { + Step37::LaplaceProblem<2> laplace_problem; + laplace_problem.run(); + } + + { + Step37::LaplaceProblem<3> laplace_problem; + laplace_problem.run(); + } +} diff --git a/tests/matrix_free/step-37-inhomogeneous-1.mpirun=4.with_p4est=true.output b/tests/matrix_free/step-37-inhomogeneous-1.mpirun=4.with_p4est=true.output new file mode 100644 index 0000000000..6121bdd29e --- /dev/null +++ b/tests/matrix_free/step-37-inhomogeneous-1.mpirun=4.with_p4est=true.output @@ -0,0 +1,150 @@ + +Cycle 0 +Number of degrees of freedom: 81 +DEAL:0:cg::Starting value 45.0807 +DEAL:0:cg::Convergence step 6 value 2.23481e-13 +max error: 0.00772526 +error ratio: 1.00000 + +Cycle 1 +Number of degrees of freedom: 289 +DEAL:0:cg::Starting value 68.5863 +DEAL:0:cg::Convergence step 6 value 4.06929e-13 +max error: 0.000971239 +error ratio: 7.95402 + +Cycle 2 +Number of degrees of freedom: 1089 +DEAL:0:cg::Starting value 95.4381 +DEAL:0:cg::Convergence step 6 value 5.78305e-13 +max error: 0.000120278 +error ratio: 8.07497 + +Cycle 3 +Number of degrees of freedom: 4225 +DEAL:0:cg::Starting value 131.732 +DEAL:0:cg::Convergence step 6 value 1.31602e-12 +max error: 1.49228e-05 +error ratio: 8.05998 + +Cycle 4 +Number of degrees of freedom: 16641 +DEAL:0:cg::Starting value 183.473 +DEAL:0:cg::Convergence step 6 value 2.60867e-12 +max error: 1.85713e-06 +error ratio: 8.03544 + +Cycle 5 +Number of degrees of freedom: 66049 +DEAL:0:cg::Starting value 257.379 +DEAL:0:cg::Convergence step 6 value 4.17156e-12 +max error: 2.31590e-07 +error ratio: 8.01904 + +Cycle 0 +Number of degrees of freedom: 125 +DEAL:0:cg::Starting value 6.44073 +DEAL:0:cg::Convergence step 6 value 2.97002e-14 +max error: 0.0499110 +error ratio: 1.00000 + +Cycle 1 +Number of degrees of freedom: 729 +DEAL:0:cg::Starting value 8.68642 +DEAL:0:cg::Convergence step 6 value 6.80139e-14 +max error: 0.00756448 +error ratio: 6.59807 + +Cycle 2 +Number of degrees of freedom: 4913 +DEAL:0:cg::Starting value 9.47362 +DEAL:0:cg::Convergence step 6 value 6.56927e-14 +max error: 0.000965192 +error ratio: 7.83728 + +Cycle 3 +Number of degrees of freedom: 35937 +DEAL:0:cg::Starting value 9.30402 +DEAL:0:cg::Convergence step 6 value 7.06917e-14 +max error: 0.000120048 +error ratio: 8.04008 + +Cycle 4 +Number of degrees of freedom: 274625 +DEAL:0:cg::Starting value 8.99023 +DEAL:0:cg::Convergence step 6 value 8.15983e-14 +max error: 1.49140e-05 +error ratio: 8.04933 + + +DEAL:1:cg::Starting value 45.0807 +DEAL:1:cg::Convergence step 6 value 2.23481e-13 +DEAL:1:cg::Starting value 68.5863 +DEAL:1:cg::Convergence step 6 value 4.06929e-13 +DEAL:1:cg::Starting value 95.4381 +DEAL:1:cg::Convergence step 6 value 5.78305e-13 +DEAL:1:cg::Starting value 131.732 +DEAL:1:cg::Convergence step 6 value 1.31602e-12 +DEAL:1:cg::Starting value 183.473 +DEAL:1:cg::Convergence step 6 value 2.60867e-12 +DEAL:1:cg::Starting value 257.379 +DEAL:1:cg::Convergence step 6 value 4.17156e-12 +DEAL:1:cg::Starting value 6.44073 +DEAL:1:cg::Convergence step 6 value 2.97002e-14 +DEAL:1:cg::Starting value 8.68642 +DEAL:1:cg::Convergence step 6 value 6.80139e-14 +DEAL:1:cg::Starting value 9.47362 +DEAL:1:cg::Convergence step 6 value 6.56927e-14 +DEAL:1:cg::Starting value 9.30402 +DEAL:1:cg::Convergence step 6 value 7.06917e-14 +DEAL:1:cg::Starting value 8.99023 +DEAL:1:cg::Convergence step 6 value 8.15983e-14 + + +DEAL:2:cg::Starting value 45.0807 +DEAL:2:cg::Convergence step 6 value 2.23481e-13 +DEAL:2:cg::Starting value 68.5863 +DEAL:2:cg::Convergence step 6 value 4.06929e-13 +DEAL:2:cg::Starting value 95.4381 +DEAL:2:cg::Convergence step 6 value 5.78305e-13 +DEAL:2:cg::Starting value 131.732 +DEAL:2:cg::Convergence step 6 value 1.31602e-12 +DEAL:2:cg::Starting value 183.473 +DEAL:2:cg::Convergence step 6 value 2.60867e-12 +DEAL:2:cg::Starting value 257.379 +DEAL:2:cg::Convergence step 6 value 4.17156e-12 +DEAL:2:cg::Starting value 6.44073 +DEAL:2:cg::Convergence step 6 value 2.97002e-14 +DEAL:2:cg::Starting value 8.68642 +DEAL:2:cg::Convergence step 6 value 6.80139e-14 +DEAL:2:cg::Starting value 9.47362 +DEAL:2:cg::Convergence step 6 value 6.56927e-14 +DEAL:2:cg::Starting value 9.30402 +DEAL:2:cg::Convergence step 6 value 7.06917e-14 +DEAL:2:cg::Starting value 8.99023 +DEAL:2:cg::Convergence step 6 value 8.15983e-14 + + +DEAL:3:cg::Starting value 45.0807 +DEAL:3:cg::Convergence step 6 value 2.23481e-13 +DEAL:3:cg::Starting value 68.5863 +DEAL:3:cg::Convergence step 6 value 4.06929e-13 +DEAL:3:cg::Starting value 95.4381 +DEAL:3:cg::Convergence step 6 value 5.78305e-13 +DEAL:3:cg::Starting value 131.732 +DEAL:3:cg::Convergence step 6 value 1.31602e-12 +DEAL:3:cg::Starting value 183.473 +DEAL:3:cg::Convergence step 6 value 2.60867e-12 +DEAL:3:cg::Starting value 257.379 +DEAL:3:cg::Convergence step 6 value 4.17156e-12 +DEAL:3:cg::Starting value 6.44073 +DEAL:3:cg::Convergence step 6 value 2.97002e-14 +DEAL:3:cg::Starting value 8.68642 +DEAL:3:cg::Convergence step 6 value 6.80139e-14 +DEAL:3:cg::Starting value 9.47362 +DEAL:3:cg::Convergence step 6 value 6.56927e-14 +DEAL:3:cg::Starting value 9.30402 +DEAL:3:cg::Convergence step 6 value 7.06917e-14 +DEAL:3:cg::Starting value 8.99023 +DEAL:3:cg::Convergence step 6 value 8.15983e-14 + -- 2.39.5