From: Jean-Paul Pelteret Date: Sun, 12 Sep 2021 18:01:31 +0000 (+0200) Subject: Add tests for inverse_operator() using Trilinos payload X-Git-Tag: v9.4.0-rc1~984^2~1 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=0ee8ef001b1e617d8177971de37f9e94d62064b0;p=dealii.git Add tests for inverse_operator() using Trilinos payload --- diff --git a/tests/lac/step-40-linear_operator_01.cc b/tests/lac/step-40-linear_operator_01.cc new file mode 100644 index 0000000000..789713191c --- /dev/null +++ b/tests/lac/step-40-linear_operator_01.cc @@ -0,0 +1,373 @@ + +// --------------------------------------------------------------------- +// +// Copyright (C) 2009 - 2021 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. +// +// --------------------------------------------------------------------- + +// A lightly adapted version of the step-40 tutorial program. +// Tests different ways of settings up a Trilinos inverse_operator(), +// specifically targeting the correct configuration of the +// TrilinosPayload::inverse_payload() when using different exemplar matrices. + + +#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 "../tests.h" + +namespace Step40 +{ + using namespace dealii; + + template + class LaplaceProblem + { + public: + LaplaceProblem(); + + void + run(); + + private: + void + setup_system(); + void + assemble_system(); + void + solve(); + void + refine_grid(); + + MPI_Comm mpi_communicator; + + parallel::distributed::Triangulation triangulation; + + FE_Q fe; + DoFHandler dof_handler; + + IndexSet locally_owned_dofs; + IndexSet locally_relevant_dofs; + + AffineConstraints constraints; + + TrilinosWrappers::SparseMatrix system_matrix; + TrilinosWrappers::MPI::Vector locally_relevant_solution; + TrilinosWrappers::MPI::Vector system_rhs; + + ConditionalOStream pcout; + }; + + template + LaplaceProblem::LaplaceProblem() + : mpi_communicator(MPI_COMM_WORLD) + , triangulation(mpi_communicator, + typename Triangulation::MeshSmoothing( + Triangulation::smoothing_on_refinement | + Triangulation::smoothing_on_coarsening)) + , fe(2) + , dof_handler(triangulation) + , pcout(Utilities::MPI::this_mpi_process(mpi_communicator) == 0 ? + deallog.get_file_stream() : + std::cout, + (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)) + {} + + template + void + LaplaceProblem::setup_system() + { + dof_handler.distribute_dofs(fe); + + locally_owned_dofs = dof_handler.locally_owned_dofs(); + DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs); + + locally_relevant_solution.reinit(locally_owned_dofs, + locally_relevant_dofs, + mpi_communicator); + system_rhs.reinit(locally_owned_dofs, mpi_communicator); + + constraints.clear(); + constraints.reinit(locally_relevant_dofs); + DoFTools::make_hanging_node_constraints(dof_handler, constraints); + VectorTools::interpolate_boundary_values(dof_handler, + 0, + Functions::ZeroFunction(), + constraints); + constraints.close(); + + DynamicSparsityPattern dsp(locally_relevant_dofs); + + DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints, false); + SparsityTools::distribute_sparsity_pattern(dsp, + dof_handler.locally_owned_dofs(), + mpi_communicator, + locally_relevant_dofs); + + system_matrix.reinit(locally_owned_dofs, + locally_owned_dofs, + dsp, + mpi_communicator); + } + + template + void + LaplaceProblem::assemble_system() + { + const QGauss quadrature_formula(fe.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); + + for (const auto &cell : dof_handler.active_cell_iterators()) + if (cell->is_locally_owned()) + { + cell_matrix = 0.; + cell_rhs = 0.; + + fe_values.reinit(cell); + + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + { + const double rhs_value = + (fe_values.quadrature_point(q_point)[1] > + 0.5 + + 0.25 * std::sin(4.0 * numbers::PI * + fe_values.quadrature_point(q_point)[0]) ? + 1. : + -1.); + + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + for (unsigned int j = 0; j < dofs_per_cell; ++j) + cell_matrix(i, j) += fe_values.shape_grad(i, q_point) * + fe_values.shape_grad(j, q_point) * + fe_values.JxW(q_point); + + cell_rhs(i) += rhs_value * // + fe_values.shape_value(i, q_point) * // + 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); + } + + system_matrix.compress(VectorOperation::add); + system_rhs.compress(VectorOperation::add); + } + + template + void + LaplaceProblem::solve() + { + TrilinosWrappers::MPI::Vector completely_distributed_solution( + locally_owned_dofs, mpi_communicator); + + SolverControl solver_control(dof_handler.n_dofs(), 1e-12, false, false); + + TrilinosWrappers::SolverCG solver(solver_control); + + TrilinosWrappers::PreconditionAMG preconditioner; + + TrilinosWrappers::PreconditionAMG::AdditionalData data; + + preconditioner.initialize(system_matrix, data); + + auto K = linear_operator(system_matrix); + auto K_inv = inverse_operator(K, solver, preconditioner); + + check_solver_within_range(completely_distributed_solution = + K_inv * system_rhs, + solver_control.last_step(), + 20, + 25); + + pcout << " Solved in " << solver_control.last_step() << " iterations." + << std::endl; + + constraints.distribute(completely_distributed_solution); + + locally_relevant_solution = completely_distributed_solution; + } + + template + void + LaplaceProblem::refine_grid() + { + triangulation.refine_global(1); + } + + template + void + LaplaceProblem::run() + { + const unsigned int n_cycles = 2; + for (unsigned int cycle = 0; cycle < n_cycles; ++cycle) + { + pcout << "Cycle " << cycle << ':' << std::endl; + + if (cycle == 0) + { + GridGenerator::hyper_cube(triangulation); + triangulation.refine_global(5); + } + else + refine_grid(); + + setup_system(); + + pcout << " Number of active cells: " + << triangulation.n_global_active_cells() << std::endl + << " "; + const auto n_locally_owned_active_cells_per_processor = + Utilities::MPI::all_gather( + triangulation.get_communicator(), + triangulation.n_locally_owned_active_cells()); + for (unsigned int i = 0; + i < Utilities::MPI::n_mpi_processes(mpi_communicator); + ++i) + pcout << n_locally_owned_active_cells_per_processor[i] << '+'; + pcout << std::endl; + + pcout << " Number of degrees of freedom: " << dof_handler.n_dofs() + << std::endl + << " "; + for (unsigned int i = 0; + i < Utilities::MPI::n_mpi_processes(mpi_communicator); + ++i) + pcout << Utilities::MPI::all_gather( + MPI_COMM_WORLD, dof_handler.n_locally_owned_dofs())[i] + << '+'; + pcout << std::endl; + + assemble_system(); + solve(); + + pcout << std::endl; + } + } +} // namespace Step40 + + +int +test_mpi() +{ + try + { + using namespace Step40; + + + { + LaplaceProblem<2> laplace_problem_2d; + laplace_problem_2d.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; +} + + + +int +main(int argc, char *argv[]) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + + if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0) + { + initlog(); + + deallog.push("mpi"); + test_mpi(); + deallog.pop(); + } + else + test_mpi(); +} diff --git a/tests/lac/step-40-linear_operator_01.with_trilinos=true.mpirun=3.with_p4est=true.output b/tests/lac/step-40-linear_operator_01.with_trilinos=true.mpirun=3.with_p4est=true.output new file mode 100644 index 0000000000..f1fc100130 --- /dev/null +++ b/tests/lac/step-40-linear_operator_01.with_trilinos=true.mpirun=3.with_p4est=true.output @@ -0,0 +1,17 @@ + +Cycle 0: + Number of active cells: 1024 + 340+344+340+ + Number of degrees of freedom: 4225 + 1453+1412+1360+ +DEAL:mpi::Solver stopped within 20 - 25 iterations + Solved in 21 iterations. + +Cycle 1: + Number of active cells: 4096 + 1364+1368+1364+ + Number of degrees of freedom: 16641 + 5645+5540+5456+ +DEAL:mpi::Solver stopped within 20 - 25 iterations + Solved in 22 iterations. + diff --git a/tests/lac/step-40-linear_operator_02.cc b/tests/lac/step-40-linear_operator_02.cc new file mode 100644 index 0000000000..b8bec22f7c --- /dev/null +++ b/tests/lac/step-40-linear_operator_02.cc @@ -0,0 +1,376 @@ + +// --------------------------------------------------------------------- +// +// Copyright (C) 2009 - 2021 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. +// +// --------------------------------------------------------------------- + +// A lightly adapted version of the step-40 tutorial program. +// Tests different ways of settings up a Trilinos inverse_operator(), +// specifically targeting the correct configuration of the +// TrilinosPayload::inverse_payload() when using different exemplar matrices. + + +#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 "../tests.h" + +namespace Step40 +{ + using namespace dealii; + + template + class LaplaceProblem + { + public: + LaplaceProblem(); + + void + run(); + + private: + void + setup_system(); + void + assemble_system(); + void + solve(); + void + refine_grid(); + + MPI_Comm mpi_communicator; + + parallel::distributed::Triangulation triangulation; + + FE_Q fe; + DoFHandler dof_handler; + + IndexSet locally_owned_dofs; + IndexSet locally_relevant_dofs; + + AffineConstraints constraints; + + TrilinosWrappers::SparseMatrix system_matrix; + TrilinosWrappers::MPI::Vector locally_relevant_solution; + TrilinosWrappers::MPI::Vector system_rhs; + + ConditionalOStream pcout; + }; + + template + LaplaceProblem::LaplaceProblem() + : mpi_communicator(MPI_COMM_WORLD) + , triangulation(mpi_communicator, + typename Triangulation::MeshSmoothing( + Triangulation::smoothing_on_refinement | + Triangulation::smoothing_on_coarsening)) + , fe(2) + , dof_handler(triangulation) + , pcout(Utilities::MPI::this_mpi_process(mpi_communicator) == 0 ? + deallog.get_file_stream() : + std::cout, + (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)) + {} + + template + void + LaplaceProblem::setup_system() + { + dof_handler.distribute_dofs(fe); + + locally_owned_dofs = dof_handler.locally_owned_dofs(); + DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs); + + locally_relevant_solution.reinit(locally_owned_dofs, + locally_relevant_dofs, + mpi_communicator); + system_rhs.reinit(locally_owned_dofs, mpi_communicator); + + constraints.clear(); + constraints.reinit(locally_relevant_dofs); + DoFTools::make_hanging_node_constraints(dof_handler, constraints); + VectorTools::interpolate_boundary_values(dof_handler, + 0, + Functions::ZeroFunction(), + constraints); + constraints.close(); + + DynamicSparsityPattern dsp(locally_relevant_dofs); + + DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints, false); + SparsityTools::distribute_sparsity_pattern(dsp, + dof_handler.locally_owned_dofs(), + mpi_communicator, + locally_relevant_dofs); + + system_matrix.reinit(locally_owned_dofs, + locally_owned_dofs, + dsp, + mpi_communicator); + } + + template + void + LaplaceProblem::assemble_system() + { + const QGauss quadrature_formula(fe.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); + + for (const auto &cell : dof_handler.active_cell_iterators()) + if (cell->is_locally_owned()) + { + cell_matrix = 0.; + cell_rhs = 0.; + + fe_values.reinit(cell); + + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + { + const double rhs_value = + (fe_values.quadrature_point(q_point)[1] > + 0.5 + + 0.25 * std::sin(4.0 * numbers::PI * + fe_values.quadrature_point(q_point)[0]) ? + 1. : + -1.); + + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + for (unsigned int j = 0; j < dofs_per_cell; ++j) + cell_matrix(i, j) += fe_values.shape_grad(i, q_point) * + fe_values.shape_grad(j, q_point) * + fe_values.JxW(q_point); + + cell_rhs(i) += rhs_value * // + fe_values.shape_value(i, q_point) * // + 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); + } + + system_matrix.compress(VectorOperation::add); + system_rhs.compress(VectorOperation::add); + } + + template + void + LaplaceProblem::solve() + { + TrilinosWrappers::MPI::Vector completely_distributed_solution( + locally_owned_dofs, mpi_communicator); + + SolverControl solver_control(dof_handler.n_dofs(), 1e-12, false, false); + + TrilinosWrappers::SolverCG solver(solver_control); + + TrilinosWrappers::PreconditionAMG preconditioner; + + TrilinosWrappers::PreconditionAMG::AdditionalData data; + + preconditioner.initialize(system_matrix, data); + + auto K = linear_operator(system_matrix); + auto prec_K = + linear_operator(system_matrix, + preconditioner); + auto K_inv = inverse_operator(K, solver, prec_K); + + check_solver_within_range(completely_distributed_solution = + K_inv * system_rhs, + solver_control.last_step(), + 20, + 25); + + pcout << " Solved in " << solver_control.last_step() << " iterations." + << std::endl; + + constraints.distribute(completely_distributed_solution); + + locally_relevant_solution = completely_distributed_solution; + } + + template + void + LaplaceProblem::refine_grid() + { + triangulation.refine_global(1); + } + + template + void + LaplaceProblem::run() + { + const unsigned int n_cycles = 2; + for (unsigned int cycle = 0; cycle < n_cycles; ++cycle) + { + pcout << "Cycle " << cycle << ':' << std::endl; + + if (cycle == 0) + { + GridGenerator::hyper_cube(triangulation); + triangulation.refine_global(5); + } + else + refine_grid(); + + setup_system(); + + pcout << " Number of active cells: " + << triangulation.n_global_active_cells() << std::endl + << " "; + const auto n_locally_owned_active_cells_per_processor = + Utilities::MPI::all_gather( + triangulation.get_communicator(), + triangulation.n_locally_owned_active_cells()); + for (unsigned int i = 0; + i < Utilities::MPI::n_mpi_processes(mpi_communicator); + ++i) + pcout << n_locally_owned_active_cells_per_processor[i] << '+'; + pcout << std::endl; + + pcout << " Number of degrees of freedom: " << dof_handler.n_dofs() + << std::endl + << " "; + for (unsigned int i = 0; + i < Utilities::MPI::n_mpi_processes(mpi_communicator); + ++i) + pcout << Utilities::MPI::all_gather( + MPI_COMM_WORLD, dof_handler.n_locally_owned_dofs())[i] + << '+'; + pcout << std::endl; + + assemble_system(); + solve(); + + pcout << std::endl; + } + } +} // namespace Step40 + + +int +test_mpi() +{ + try + { + using namespace Step40; + + + { + LaplaceProblem<2> laplace_problem_2d; + laplace_problem_2d.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; +} + + + +int +main(int argc, char *argv[]) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + + if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0) + { + initlog(); + + deallog.push("mpi"); + test_mpi(); + deallog.pop(); + } + else + test_mpi(); +} diff --git a/tests/lac/step-40-linear_operator_02.with_trilinos=true.mpirun=3.with_p4est=true.output b/tests/lac/step-40-linear_operator_02.with_trilinos=true.mpirun=3.with_p4est=true.output new file mode 100644 index 0000000000..f1fc100130 --- /dev/null +++ b/tests/lac/step-40-linear_operator_02.with_trilinos=true.mpirun=3.with_p4est=true.output @@ -0,0 +1,17 @@ + +Cycle 0: + Number of active cells: 1024 + 340+344+340+ + Number of degrees of freedom: 4225 + 1453+1412+1360+ +DEAL:mpi::Solver stopped within 20 - 25 iterations + Solved in 21 iterations. + +Cycle 1: + Number of active cells: 4096 + 1364+1368+1364+ + Number of degrees of freedom: 16641 + 5645+5540+5456+ +DEAL:mpi::Solver stopped within 20 - 25 iterations + Solved in 22 iterations. + diff --git a/tests/lac/step-40-linear_operator_03.cc b/tests/lac/step-40-linear_operator_03.cc new file mode 100644 index 0000000000..d8438bb62c --- /dev/null +++ b/tests/lac/step-40-linear_operator_03.cc @@ -0,0 +1,377 @@ + + +// --------------------------------------------------------------------- +// +// Copyright (C) 2009 - 2021 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. +// +// --------------------------------------------------------------------- + +// A lightly adapted version of the step-40 tutorial program. +// Tests different ways of settings up a Trilinos inverse_operator(), +// specifically targeting the correct configuration of the +// TrilinosPayload::inverse_payload() when using different exemplar matrices. + + +#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 "../tests.h" + +namespace Step40 +{ + using namespace dealii; + + template + class LaplaceProblem + { + public: + LaplaceProblem(); + + void + run(); + + private: + void + setup_system(); + void + assemble_system(); + void + solve(); + void + refine_grid(); + + MPI_Comm mpi_communicator; + + parallel::distributed::Triangulation triangulation; + + FE_Q fe; + DoFHandler dof_handler; + + IndexSet locally_owned_dofs; + IndexSet locally_relevant_dofs; + + AffineConstraints constraints; + + TrilinosWrappers::SparseMatrix system_matrix; + TrilinosWrappers::MPI::Vector locally_relevant_solution; + TrilinosWrappers::MPI::Vector system_rhs; + + ConditionalOStream pcout; + }; + + template + LaplaceProblem::LaplaceProblem() + : mpi_communicator(MPI_COMM_WORLD) + , triangulation(mpi_communicator, + typename Triangulation::MeshSmoothing( + Triangulation::smoothing_on_refinement | + Triangulation::smoothing_on_coarsening)) + , fe(2) + , dof_handler(triangulation) + , pcout(Utilities::MPI::this_mpi_process(mpi_communicator) == 0 ? + deallog.get_file_stream() : + std::cout, + (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)) + {} + + template + void + LaplaceProblem::setup_system() + { + dof_handler.distribute_dofs(fe); + + locally_owned_dofs = dof_handler.locally_owned_dofs(); + DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs); + + locally_relevant_solution.reinit(locally_owned_dofs, + locally_relevant_dofs, + mpi_communicator); + system_rhs.reinit(locally_owned_dofs, mpi_communicator); + + constraints.clear(); + constraints.reinit(locally_relevant_dofs); + DoFTools::make_hanging_node_constraints(dof_handler, constraints); + VectorTools::interpolate_boundary_values(dof_handler, + 0, + Functions::ZeroFunction(), + constraints); + constraints.close(); + + DynamicSparsityPattern dsp(locally_relevant_dofs); + + DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints, false); + SparsityTools::distribute_sparsity_pattern(dsp, + dof_handler.locally_owned_dofs(), + mpi_communicator, + locally_relevant_dofs); + + system_matrix.reinit(locally_owned_dofs, + locally_owned_dofs, + dsp, + mpi_communicator); + } + + template + void + LaplaceProblem::assemble_system() + { + const QGauss quadrature_formula(fe.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); + + for (const auto &cell : dof_handler.active_cell_iterators()) + if (cell->is_locally_owned()) + { + cell_matrix = 0.; + cell_rhs = 0.; + + fe_values.reinit(cell); + + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + { + const double rhs_value = + (fe_values.quadrature_point(q_point)[1] > + 0.5 + + 0.25 * std::sin(4.0 * numbers::PI * + fe_values.quadrature_point(q_point)[0]) ? + 1. : + -1.); + + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + for (unsigned int j = 0; j < dofs_per_cell; ++j) + cell_matrix(i, j) += fe_values.shape_grad(i, q_point) * + fe_values.shape_grad(j, q_point) * + fe_values.JxW(q_point); + + cell_rhs(i) += rhs_value * // + fe_values.shape_value(i, q_point) * // + 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); + } + + system_matrix.compress(VectorOperation::add); + system_rhs.compress(VectorOperation::add); + } + + template + void + LaplaceProblem::solve() + { + TrilinosWrappers::MPI::Vector completely_distributed_solution( + locally_owned_dofs, mpi_communicator); + + SolverControl solver_control(dof_handler.n_dofs(), 1e-12, false, false); + + TrilinosWrappers::SolverCG solver(solver_control); + + TrilinosWrappers::PreconditionAMG preconditioner; + + TrilinosWrappers::PreconditionAMG::AdditionalData data; + + preconditioner.initialize(system_matrix, data); + + auto K = linear_operator(system_matrix); + auto prec_K = + TrilinosWrappers::linear_operator( + system_matrix, preconditioner); + auto K_inv = inverse_operator(K, solver, prec_K); + + check_solver_within_range(completely_distributed_solution = + K_inv * system_rhs, + solver_control.last_step(), + 20, + 25); + + pcout << " Solved in " << solver_control.last_step() << " iterations." + << std::endl; + + constraints.distribute(completely_distributed_solution); + + locally_relevant_solution = completely_distributed_solution; + } + + template + void + LaplaceProblem::refine_grid() + { + triangulation.refine_global(1); + } + + template + void + LaplaceProblem::run() + { + const unsigned int n_cycles = 2; + for (unsigned int cycle = 0; cycle < n_cycles; ++cycle) + { + pcout << "Cycle " << cycle << ':' << std::endl; + + if (cycle == 0) + { + GridGenerator::hyper_cube(triangulation); + triangulation.refine_global(5); + } + else + refine_grid(); + + setup_system(); + + pcout << " Number of active cells: " + << triangulation.n_global_active_cells() << std::endl + << " "; + const auto n_locally_owned_active_cells_per_processor = + Utilities::MPI::all_gather( + triangulation.get_communicator(), + triangulation.n_locally_owned_active_cells()); + for (unsigned int i = 0; + i < Utilities::MPI::n_mpi_processes(mpi_communicator); + ++i) + pcout << n_locally_owned_active_cells_per_processor[i] << '+'; + pcout << std::endl; + + pcout << " Number of degrees of freedom: " << dof_handler.n_dofs() + << std::endl + << " "; + for (unsigned int i = 0; + i < Utilities::MPI::n_mpi_processes(mpi_communicator); + ++i) + pcout << Utilities::MPI::all_gather( + MPI_COMM_WORLD, dof_handler.n_locally_owned_dofs())[i] + << '+'; + pcout << std::endl; + + assemble_system(); + solve(); + + pcout << std::endl; + } + } +} // namespace Step40 + + +int +test_mpi() +{ + try + { + using namespace Step40; + + + { + LaplaceProblem<2> laplace_problem_2d; + laplace_problem_2d.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; +} + + + +int +main(int argc, char *argv[]) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + + if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0) + { + initlog(); + + deallog.push("mpi"); + test_mpi(); + deallog.pop(); + } + else + test_mpi(); +} diff --git a/tests/lac/step-40-linear_operator_03.with_trilinos=true.mpirun=3.with_p4est=true.output b/tests/lac/step-40-linear_operator_03.with_trilinos=true.mpirun=3.with_p4est=true.output new file mode 100644 index 0000000000..f1fc100130 --- /dev/null +++ b/tests/lac/step-40-linear_operator_03.with_trilinos=true.mpirun=3.with_p4est=true.output @@ -0,0 +1,17 @@ + +Cycle 0: + Number of active cells: 1024 + 340+344+340+ + Number of degrees of freedom: 4225 + 1453+1412+1360+ +DEAL:mpi::Solver stopped within 20 - 25 iterations + Solved in 21 iterations. + +Cycle 1: + Number of active cells: 4096 + 1364+1368+1364+ + Number of degrees of freedom: 16641 + 5645+5540+5456+ +DEAL:mpi::Solver stopped within 20 - 25 iterations + Solved in 22 iterations. + diff --git a/tests/lac/step-40-linear_operator_04.cc b/tests/lac/step-40-linear_operator_04.cc new file mode 100644 index 0000000000..d7fe515bf4 --- /dev/null +++ b/tests/lac/step-40-linear_operator_04.cc @@ -0,0 +1,377 @@ + + +// --------------------------------------------------------------------- +// +// Copyright (C) 2009 - 2021 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. +// +// --------------------------------------------------------------------- + +// A lightly adapted version of the step-40 tutorial program. +// Tests different ways of settings up a Trilinos inverse_operator(), +// specifically targeting the correct configuration of the +// TrilinosPayload::inverse_payload() when using different exemplar matrices. + + +#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 "../tests.h" + +namespace Step40 +{ + using namespace dealii; + + template + class LaplaceProblem + { + public: + LaplaceProblem(); + + void + run(); + + private: + void + setup_system(); + void + assemble_system(); + void + solve(); + void + refine_grid(); + + MPI_Comm mpi_communicator; + + parallel::distributed::Triangulation triangulation; + + FE_Q fe; + DoFHandler dof_handler; + + IndexSet locally_owned_dofs; + IndexSet locally_relevant_dofs; + + AffineConstraints constraints; + + TrilinosWrappers::SparseMatrix system_matrix; + TrilinosWrappers::MPI::Vector locally_relevant_solution; + TrilinosWrappers::MPI::Vector system_rhs; + + ConditionalOStream pcout; + }; + + template + LaplaceProblem::LaplaceProblem() + : mpi_communicator(MPI_COMM_WORLD) + , triangulation(mpi_communicator, + typename Triangulation::MeshSmoothing( + Triangulation::smoothing_on_refinement | + Triangulation::smoothing_on_coarsening)) + , fe(2) + , dof_handler(triangulation) + , pcout(Utilities::MPI::this_mpi_process(mpi_communicator) == 0 ? + deallog.get_file_stream() : + std::cout, + (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)) + {} + + template + void + LaplaceProblem::setup_system() + { + dof_handler.distribute_dofs(fe); + + locally_owned_dofs = dof_handler.locally_owned_dofs(); + DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs); + + locally_relevant_solution.reinit(locally_owned_dofs, + locally_relevant_dofs, + mpi_communicator); + system_rhs.reinit(locally_owned_dofs, mpi_communicator); + + constraints.clear(); + constraints.reinit(locally_relevant_dofs); + DoFTools::make_hanging_node_constraints(dof_handler, constraints); + VectorTools::interpolate_boundary_values(dof_handler, + 0, + Functions::ZeroFunction(), + constraints); + constraints.close(); + + DynamicSparsityPattern dsp(locally_relevant_dofs); + + DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints, false); + SparsityTools::distribute_sparsity_pattern(dsp, + dof_handler.locally_owned_dofs(), + mpi_communicator, + locally_relevant_dofs); + + system_matrix.reinit(locally_owned_dofs, + locally_owned_dofs, + dsp, + mpi_communicator); + } + + template + void + LaplaceProblem::assemble_system() + { + const QGauss quadrature_formula(fe.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); + + for (const auto &cell : dof_handler.active_cell_iterators()) + if (cell->is_locally_owned()) + { + cell_matrix = 0.; + cell_rhs = 0.; + + fe_values.reinit(cell); + + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + { + const double rhs_value = + (fe_values.quadrature_point(q_point)[1] > + 0.5 + + 0.25 * std::sin(4.0 * numbers::PI * + fe_values.quadrature_point(q_point)[0]) ? + 1. : + -1.); + + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + for (unsigned int j = 0; j < dofs_per_cell; ++j) + cell_matrix(i, j) += fe_values.shape_grad(i, q_point) * + fe_values.shape_grad(j, q_point) * + fe_values.JxW(q_point); + + cell_rhs(i) += rhs_value * // + fe_values.shape_value(i, q_point) * // + 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); + } + + system_matrix.compress(VectorOperation::add); + system_rhs.compress(VectorOperation::add); + } + + template + void + LaplaceProblem::solve() + { + TrilinosWrappers::MPI::Vector completely_distributed_solution( + locally_owned_dofs, mpi_communicator); + + SolverControl solver_control(dof_handler.n_dofs(), 1e-12, false, false); + + TrilinosWrappers::SolverCG solver(solver_control); + + TrilinosWrappers::PreconditionAMG preconditioner; + + TrilinosWrappers::PreconditionAMG::AdditionalData data; + + preconditioner.initialize(system_matrix, data); + + auto K = linear_operator(system_matrix); + auto prec_K = + TrilinosWrappers::linear_operator( + K, preconditioner); + auto K_inv = inverse_operator(K, solver, prec_K); + + check_solver_within_range(completely_distributed_solution = + K_inv * system_rhs, + solver_control.last_step(), + 20, + 25); + + pcout << " Solved in " << solver_control.last_step() << " iterations." + << std::endl; + + constraints.distribute(completely_distributed_solution); + + locally_relevant_solution = completely_distributed_solution; + } + + template + void + LaplaceProblem::refine_grid() + { + triangulation.refine_global(1); + } + + template + void + LaplaceProblem::run() + { + const unsigned int n_cycles = 2; + for (unsigned int cycle = 0; cycle < n_cycles; ++cycle) + { + pcout << "Cycle " << cycle << ':' << std::endl; + + if (cycle == 0) + { + GridGenerator::hyper_cube(triangulation); + triangulation.refine_global(5); + } + else + refine_grid(); + + setup_system(); + + pcout << " Number of active cells: " + << triangulation.n_global_active_cells() << std::endl + << " "; + const auto n_locally_owned_active_cells_per_processor = + Utilities::MPI::all_gather( + triangulation.get_communicator(), + triangulation.n_locally_owned_active_cells()); + for (unsigned int i = 0; + i < Utilities::MPI::n_mpi_processes(mpi_communicator); + ++i) + pcout << n_locally_owned_active_cells_per_processor[i] << '+'; + pcout << std::endl; + + pcout << " Number of degrees of freedom: " << dof_handler.n_dofs() + << std::endl + << " "; + for (unsigned int i = 0; + i < Utilities::MPI::n_mpi_processes(mpi_communicator); + ++i) + pcout << Utilities::MPI::all_gather( + MPI_COMM_WORLD, dof_handler.n_locally_owned_dofs())[i] + << '+'; + pcout << std::endl; + + assemble_system(); + solve(); + + pcout << std::endl; + } + } +} // namespace Step40 + + +int +test_mpi() +{ + try + { + using namespace Step40; + + + { + LaplaceProblem<2> laplace_problem_2d; + laplace_problem_2d.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; +} + + + +int +main(int argc, char *argv[]) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + + if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0) + { + initlog(); + + deallog.push("mpi"); + test_mpi(); + deallog.pop(); + } + else + test_mpi(); +} diff --git a/tests/lac/step-40-linear_operator_04.with_trilinos=true.mpirun=3.with_p4est=true.output b/tests/lac/step-40-linear_operator_04.with_trilinos=true.mpirun=3.with_p4est=true.output new file mode 100644 index 0000000000..f1fc100130 --- /dev/null +++ b/tests/lac/step-40-linear_operator_04.with_trilinos=true.mpirun=3.with_p4est=true.output @@ -0,0 +1,17 @@ + +Cycle 0: + Number of active cells: 1024 + 340+344+340+ + Number of degrees of freedom: 4225 + 1453+1412+1360+ +DEAL:mpi::Solver stopped within 20 - 25 iterations + Solved in 21 iterations. + +Cycle 1: + Number of active cells: 4096 + 1364+1368+1364+ + Number of degrees of freedom: 16641 + 5645+5540+5456+ +DEAL:mpi::Solver stopped within 20 - 25 iterations + Solved in 22 iterations. + diff --git a/tests/lac/step-40-linear_operator_05.cc b/tests/lac/step-40-linear_operator_05.cc new file mode 100644 index 0000000000..d752d26dd3 --- /dev/null +++ b/tests/lac/step-40-linear_operator_05.cc @@ -0,0 +1,375 @@ + + +// --------------------------------------------------------------------- +// +// Copyright (C) 2009 - 2021 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. +// +// --------------------------------------------------------------------- + +// A lightly adapted version of the step-40 tutorial program. +// Tests different ways of settings up a Trilinos inverse_operator(), +// specifically targeting the correct configuration of the +// TrilinosPayload::inverse_payload() when using different exemplar matrices. + + +#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 "../tests.h" + +namespace Step40 +{ + using namespace dealii; + + template + class LaplaceProblem + { + public: + LaplaceProblem(); + + void + run(); + + private: + void + setup_system(); + void + assemble_system(); + void + solve(); + void + refine_grid(); + + MPI_Comm mpi_communicator; + + parallel::distributed::Triangulation triangulation; + + FE_Q fe; + DoFHandler dof_handler; + + IndexSet locally_owned_dofs; + IndexSet locally_relevant_dofs; + + AffineConstraints constraints; + + TrilinosWrappers::SparseMatrix system_matrix; + TrilinosWrappers::MPI::Vector locally_relevant_solution; + TrilinosWrappers::MPI::Vector system_rhs; + + ConditionalOStream pcout; + }; + + template + LaplaceProblem::LaplaceProblem() + : mpi_communicator(MPI_COMM_WORLD) + , triangulation(mpi_communicator, + typename Triangulation::MeshSmoothing( + Triangulation::smoothing_on_refinement | + Triangulation::smoothing_on_coarsening)) + , fe(2) + , dof_handler(triangulation) + , pcout(Utilities::MPI::this_mpi_process(mpi_communicator) == 0 ? + deallog.get_file_stream() : + std::cout, + (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)) + {} + + template + void + LaplaceProblem::setup_system() + { + dof_handler.distribute_dofs(fe); + + locally_owned_dofs = dof_handler.locally_owned_dofs(); + DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs); + + locally_relevant_solution.reinit(locally_owned_dofs, + locally_relevant_dofs, + mpi_communicator); + system_rhs.reinit(locally_owned_dofs, mpi_communicator); + + constraints.clear(); + constraints.reinit(locally_relevant_dofs); + DoFTools::make_hanging_node_constraints(dof_handler, constraints); + VectorTools::interpolate_boundary_values(dof_handler, + 0, + Functions::ZeroFunction(), + constraints); + constraints.close(); + + DynamicSparsityPattern dsp(locally_relevant_dofs); + + DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints, false); + SparsityTools::distribute_sparsity_pattern(dsp, + dof_handler.locally_owned_dofs(), + mpi_communicator, + locally_relevant_dofs); + + system_matrix.reinit(locally_owned_dofs, + locally_owned_dofs, + dsp, + mpi_communicator); + } + + template + void + LaplaceProblem::assemble_system() + { + const QGauss quadrature_formula(fe.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); + + for (const auto &cell : dof_handler.active_cell_iterators()) + if (cell->is_locally_owned()) + { + cell_matrix = 0.; + cell_rhs = 0.; + + fe_values.reinit(cell); + + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + { + const double rhs_value = + (fe_values.quadrature_point(q_point)[1] > + 0.5 + + 0.25 * std::sin(4.0 * numbers::PI * + fe_values.quadrature_point(q_point)[0]) ? + 1. : + -1.); + + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + for (unsigned int j = 0; j < dofs_per_cell; ++j) + cell_matrix(i, j) += fe_values.shape_grad(i, q_point) * + fe_values.shape_grad(j, q_point) * + fe_values.JxW(q_point); + + cell_rhs(i) += rhs_value * // + fe_values.shape_value(i, q_point) * // + 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); + } + + system_matrix.compress(VectorOperation::add); + system_rhs.compress(VectorOperation::add); + } + + template + void + LaplaceProblem::solve() + { + TrilinosWrappers::MPI::Vector completely_distributed_solution( + locally_owned_dofs, mpi_communicator); + + SolverControl solver_control(dof_handler.n_dofs(), 1e-12, false, false); + + TrilinosWrappers::SolverCG solver(solver_control); + + TrilinosWrappers::PreconditionAMG preconditioner; + + TrilinosWrappers::PreconditionAMG::AdditionalData data; + + preconditioner.initialize(system_matrix, data); + + auto K = linear_operator(system_matrix); + auto prec_K = TrilinosWrappers::linear_operator(K, preconditioner); + auto K_inv = inverse_operator(K, solver, prec_K); + + check_solver_within_range(completely_distributed_solution = + K_inv * system_rhs, + solver_control.last_step(), + 20, + 25); + + pcout << " Solved in " << solver_control.last_step() << " iterations." + << std::endl; + + constraints.distribute(completely_distributed_solution); + + locally_relevant_solution = completely_distributed_solution; + } + + template + void + LaplaceProblem::refine_grid() + { + triangulation.refine_global(1); + } + + template + void + LaplaceProblem::run() + { + const unsigned int n_cycles = 2; + for (unsigned int cycle = 0; cycle < n_cycles; ++cycle) + { + pcout << "Cycle " << cycle << ':' << std::endl; + + if (cycle == 0) + { + GridGenerator::hyper_cube(triangulation); + triangulation.refine_global(5); + } + else + refine_grid(); + + setup_system(); + + pcout << " Number of active cells: " + << triangulation.n_global_active_cells() << std::endl + << " "; + const auto n_locally_owned_active_cells_per_processor = + Utilities::MPI::all_gather( + triangulation.get_communicator(), + triangulation.n_locally_owned_active_cells()); + for (unsigned int i = 0; + i < Utilities::MPI::n_mpi_processes(mpi_communicator); + ++i) + pcout << n_locally_owned_active_cells_per_processor[i] << '+'; + pcout << std::endl; + + pcout << " Number of degrees of freedom: " << dof_handler.n_dofs() + << std::endl + << " "; + for (unsigned int i = 0; + i < Utilities::MPI::n_mpi_processes(mpi_communicator); + ++i) + pcout << Utilities::MPI::all_gather( + MPI_COMM_WORLD, dof_handler.n_locally_owned_dofs())[i] + << '+'; + pcout << std::endl; + + assemble_system(); + solve(); + + pcout << std::endl; + } + } +} // namespace Step40 + + +int +test_mpi() +{ + try + { + using namespace Step40; + + + { + LaplaceProblem<2> laplace_problem_2d; + laplace_problem_2d.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; +} + + + +int +main(int argc, char *argv[]) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + + if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0) + { + initlog(); + + deallog.push("mpi"); + test_mpi(); + deallog.pop(); + } + else + test_mpi(); +} diff --git a/tests/lac/step-40-linear_operator_05.with_trilinos=true.mpirun=3.with_p4est=true.output b/tests/lac/step-40-linear_operator_05.with_trilinos=true.mpirun=3.with_p4est=true.output new file mode 100644 index 0000000000..f1fc100130 --- /dev/null +++ b/tests/lac/step-40-linear_operator_05.with_trilinos=true.mpirun=3.with_p4est=true.output @@ -0,0 +1,17 @@ + +Cycle 0: + Number of active cells: 1024 + 340+344+340+ + Number of degrees of freedom: 4225 + 1453+1412+1360+ +DEAL:mpi::Solver stopped within 20 - 25 iterations + Solved in 21 iterations. + +Cycle 1: + Number of active cells: 4096 + 1364+1368+1364+ + Number of degrees of freedom: 16641 + 5645+5540+5456+ +DEAL:mpi::Solver stopped within 20 - 25 iterations + Solved in 22 iterations. + diff --git a/tests/lac/step-40-linear_operator_06.cc b/tests/lac/step-40-linear_operator_06.cc new file mode 100644 index 0000000000..91f3f253d8 --- /dev/null +++ b/tests/lac/step-40-linear_operator_06.cc @@ -0,0 +1,375 @@ + + +// --------------------------------------------------------------------- +// +// Copyright (C) 2009 - 2021 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. +// +// --------------------------------------------------------------------- + +// A lightly adapted version of the step-40 tutorial program. +// Tests different ways of settings up a Trilinos inverse_operator(), +// specifically targeting the correct configuration of the +// TrilinosPayload::inverse_payload() when using different exemplar matrices. + + +#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 "../tests.h" + +namespace Step40 +{ + using namespace dealii; + + template + class LaplaceProblem + { + public: + LaplaceProblem(); + + void + run(); + + private: + void + setup_system(); + void + assemble_system(); + void + solve(); + void + refine_grid(); + + MPI_Comm mpi_communicator; + + parallel::distributed::Triangulation triangulation; + + FE_Q fe; + DoFHandler dof_handler; + + IndexSet locally_owned_dofs; + IndexSet locally_relevant_dofs; + + AffineConstraints constraints; + + TrilinosWrappers::SparseMatrix system_matrix; + TrilinosWrappers::MPI::Vector locally_relevant_solution; + TrilinosWrappers::MPI::Vector system_rhs; + + ConditionalOStream pcout; + }; + + template + LaplaceProblem::LaplaceProblem() + : mpi_communicator(MPI_COMM_WORLD) + , triangulation(mpi_communicator, + typename Triangulation::MeshSmoothing( + Triangulation::smoothing_on_refinement | + Triangulation::smoothing_on_coarsening)) + , fe(2) + , dof_handler(triangulation) + , pcout(Utilities::MPI::this_mpi_process(mpi_communicator) == 0 ? + deallog.get_file_stream() : + std::cout, + (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)) + {} + + template + void + LaplaceProblem::setup_system() + { + dof_handler.distribute_dofs(fe); + + locally_owned_dofs = dof_handler.locally_owned_dofs(); + DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs); + + locally_relevant_solution.reinit(locally_owned_dofs, + locally_relevant_dofs, + mpi_communicator); + system_rhs.reinit(locally_owned_dofs, mpi_communicator); + + constraints.clear(); + constraints.reinit(locally_relevant_dofs); + DoFTools::make_hanging_node_constraints(dof_handler, constraints); + VectorTools::interpolate_boundary_values(dof_handler, + 0, + Functions::ZeroFunction(), + constraints); + constraints.close(); + + DynamicSparsityPattern dsp(locally_relevant_dofs); + + DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints, false); + SparsityTools::distribute_sparsity_pattern(dsp, + dof_handler.locally_owned_dofs(), + mpi_communicator, + locally_relevant_dofs); + + system_matrix.reinit(locally_owned_dofs, + locally_owned_dofs, + dsp, + mpi_communicator); + } + + template + void + LaplaceProblem::assemble_system() + { + const QGauss quadrature_formula(fe.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); + + for (const auto &cell : dof_handler.active_cell_iterators()) + if (cell->is_locally_owned()) + { + cell_matrix = 0.; + cell_rhs = 0.; + + fe_values.reinit(cell); + + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + { + const double rhs_value = + (fe_values.quadrature_point(q_point)[1] > + 0.5 + + 0.25 * std::sin(4.0 * numbers::PI * + fe_values.quadrature_point(q_point)[0]) ? + 1. : + -1.); + + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + for (unsigned int j = 0; j < dofs_per_cell; ++j) + cell_matrix(i, j) += fe_values.shape_grad(i, q_point) * + fe_values.shape_grad(j, q_point) * + fe_values.JxW(q_point); + + cell_rhs(i) += rhs_value * // + fe_values.shape_value(i, q_point) * // + 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); + } + + system_matrix.compress(VectorOperation::add); + system_rhs.compress(VectorOperation::add); + } + + template + void + LaplaceProblem::solve() + { + TrilinosWrappers::MPI::Vector completely_distributed_solution( + locally_owned_dofs, mpi_communicator); + + SolverControl solver_control(dof_handler.n_dofs(), 1e-12, false, false); + + TrilinosWrappers::SolverCG solver(solver_control); + + TrilinosWrappers::PreconditionAMG preconditioner; + + TrilinosWrappers::PreconditionAMG::AdditionalData data; + + preconditioner.initialize(system_matrix, data); + + auto K = linear_operator(system_matrix); + auto prec_K = linear_operator(K, preconditioner); + auto K_inv = inverse_operator(K, solver, prec_K); + + check_solver_within_range(completely_distributed_solution = + K_inv * system_rhs, + solver_control.last_step(), + 20, + 25); + + pcout << " Solved in " << solver_control.last_step() << " iterations." + << std::endl; + + constraints.distribute(completely_distributed_solution); + + locally_relevant_solution = completely_distributed_solution; + } + + template + void + LaplaceProblem::refine_grid() + { + triangulation.refine_global(1); + } + + template + void + LaplaceProblem::run() + { + const unsigned int n_cycles = 2; + for (unsigned int cycle = 0; cycle < n_cycles; ++cycle) + { + pcout << "Cycle " << cycle << ':' << std::endl; + + if (cycle == 0) + { + GridGenerator::hyper_cube(triangulation); + triangulation.refine_global(5); + } + else + refine_grid(); + + setup_system(); + + pcout << " Number of active cells: " + << triangulation.n_global_active_cells() << std::endl + << " "; + const auto n_locally_owned_active_cells_per_processor = + Utilities::MPI::all_gather( + triangulation.get_communicator(), + triangulation.n_locally_owned_active_cells()); + for (unsigned int i = 0; + i < Utilities::MPI::n_mpi_processes(mpi_communicator); + ++i) + pcout << n_locally_owned_active_cells_per_processor[i] << '+'; + pcout << std::endl; + + pcout << " Number of degrees of freedom: " << dof_handler.n_dofs() + << std::endl + << " "; + for (unsigned int i = 0; + i < Utilities::MPI::n_mpi_processes(mpi_communicator); + ++i) + pcout << Utilities::MPI::all_gather( + MPI_COMM_WORLD, dof_handler.n_locally_owned_dofs())[i] + << '+'; + pcout << std::endl; + + assemble_system(); + solve(); + + pcout << std::endl; + } + } +} // namespace Step40 + + +int +test_mpi() +{ + try + { + using namespace Step40; + + + { + LaplaceProblem<2> laplace_problem_2d; + laplace_problem_2d.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; +} + + + +int +main(int argc, char *argv[]) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + + if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0) + { + initlog(); + + deallog.push("mpi"); + test_mpi(); + deallog.pop(); + } + else + test_mpi(); +} diff --git a/tests/lac/step-40-linear_operator_06.with_trilinos=true.mpirun=3.with_p4est=true.output b/tests/lac/step-40-linear_operator_06.with_trilinos=true.mpirun=3.with_p4est=true.output new file mode 100644 index 0000000000..f1fc100130 --- /dev/null +++ b/tests/lac/step-40-linear_operator_06.with_trilinos=true.mpirun=3.with_p4est=true.output @@ -0,0 +1,17 @@ + +Cycle 0: + Number of active cells: 1024 + 340+344+340+ + Number of degrees of freedom: 4225 + 1453+1412+1360+ +DEAL:mpi::Solver stopped within 20 - 25 iterations + Solved in 21 iterations. + +Cycle 1: + Number of active cells: 4096 + 1364+1368+1364+ + Number of degrees of freedom: 16641 + 5645+5540+5456+ +DEAL:mpi::Solver stopped within 20 - 25 iterations + Solved in 22 iterations. +