From 91b02361986948ffd0bc82daebd17ef7fe8576d7 Mon Sep 17 00:00:00 2001 From: Jean-Paul Pelteret Date: Fri, 15 Mar 2019 09:40:03 +0100 Subject: [PATCH] Remove dummy vector when initializing AMG constant modes --- source/lac/trilinos_precondition_ml.cc | 14 +- tests/trilinos/precondition_amg_dgp_03.cc | 273 ++++++++++++++++++ ...n_amg_dgp_03.with_64bit_indices=off.output | 5 + 3 files changed, 280 insertions(+), 12 deletions(-) create mode 100644 tests/trilinos/precondition_amg_dgp_03.cc create mode 100644 tests/trilinos/precondition_amg_dgp_03.with_64bit_indices=off.output diff --git a/source/lac/trilinos_precondition_ml.cc b/source/lac/trilinos_precondition_ml.cc index 848315efaa..1e65a430e9 100644 --- a/source/lac/trilinos_precondition_ml.cc +++ b/source/lac/trilinos_precondition_ml.cc @@ -178,18 +178,8 @@ namespace TrilinosWrappers parameter_list.set("null space: type", "pre-computed"); parameter_list.set("null space: dimension", distributed_constant_modes.NumVectors()); - if (my_size > 0) - parameter_list.set("null space: vectors", - distributed_constant_modes.Values()); - else - { - // We need to set a valid pointer to data even if there is no data - // on the current processor. Therefore, pass a dummy in that case - static std::vector dummy; - if (dummy.size() != constant_modes_dimension) - dummy.resize(constant_modes_dimension); - parameter_list.set("null space: vectors", dummy.data()); - } + parameter_list.set("null space: vectors", + distributed_constant_modes.Values()); } } diff --git a/tests/trilinos/precondition_amg_dgp_03.cc b/tests/trilinos/precondition_amg_dgp_03.cc new file mode 100644 index 0000000000..1fcf871403 --- /dev/null +++ b/tests/trilinos/precondition_amg_dgp_03.cc @@ -0,0 +1,273 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2013 - 2019 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. +// +// --------------------------------------------------------------------- + + + +// solves a 2D Poisson equation for linear FE_DGP elements (SIP +// discretization) with AMG preconditioner. +// +// This variant of the test trilinos/precondition_amg_dgp_01 does not initialize +// the constant modes vector passes to the AMG preconditioner. + +#include + +#include +#include + +#include + +#include +#include + +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include + +#include "../tests.h" + + + +template +class MatrixIntegrator : public MeshWorker::LocalIntegrator +{ +public: + void + cell(MeshWorker::DoFInfo & dinfo, + typename MeshWorker::IntegrationInfo &info) const; + void + boundary(MeshWorker::DoFInfo & dinfo, + typename MeshWorker::IntegrationInfo &info) const; + void + face(MeshWorker::DoFInfo & dinfo1, + MeshWorker::DoFInfo & dinfo2, + typename MeshWorker::IntegrationInfo &info1, + typename MeshWorker::IntegrationInfo &info2) const; +}; + +template +void +MatrixIntegrator::cell( + MeshWorker::DoFInfo & dinfo, + typename MeshWorker::IntegrationInfo &info) const +{ + LocalIntegrators::Laplace::cell_matrix(dinfo.matrix(0, false).matrix, + info.fe_values()); +} + +template +void +MatrixIntegrator::boundary( + MeshWorker::DoFInfo & dinfo, + typename MeshWorker::IntegrationInfo &info) const +{ + const unsigned int deg = info.fe_values(0).get_fe().degree; + LocalIntegrators::Laplace ::nitsche_matrix( + dinfo.matrix(0, false).matrix, + info.fe_values(0), + LocalIntegrators::Laplace::compute_penalty(dinfo, dinfo, deg, deg)); +} + +template +void +MatrixIntegrator::face( + MeshWorker::DoFInfo & dinfo1, + MeshWorker::DoFInfo & dinfo2, + typename MeshWorker::IntegrationInfo &info1, + typename MeshWorker::IntegrationInfo &info2) const +{ + const unsigned int deg = info1.fe_values(0).get_fe().degree; + LocalIntegrators::Laplace ::ip_matrix( + dinfo1.matrix(0, false).matrix, + dinfo1.matrix(0, true).matrix, + dinfo2.matrix(0, true).matrix, + dinfo2.matrix(0, false).matrix, + info1.fe_values(0), + info2.fe_values(0), + LocalIntegrators::Laplace::compute_penalty(dinfo1, dinfo2, deg, deg)); +} + + +template +class Step4 +{ +public: + Step4(); + void + run(); + +private: + void + make_grid(); + void + setup_system(); + void + solve(); + + Triangulation triangulation; + FE_DGP fe; + DoFHandler dof_handler; + + TrilinosWrappers::SparseMatrix system_matrix; + + Vector solution; + Vector system_rhs; +}; + + + +template +Step4::Step4() + : fe(1) + , dof_handler(triangulation) +{} + + +template +void +Step4::make_grid() +{ + GridGenerator::hyper_cube(triangulation, -1, 1); + triangulation.refine_global(6); +} + + + +template +void +Step4::setup_system() +{ + dof_handler.distribute_dofs(fe); + + DynamicSparsityPattern c_sparsity(dof_handler.n_dofs()); + DoFTools::make_flux_sparsity_pattern(dof_handler, c_sparsity); + system_matrix.reinit(c_sparsity); + + solution.reinit(dof_handler.n_dofs()); + system_rhs.reinit(dof_handler.n_dofs()); + + MappingQGeneric mapping(1); + MeshWorker::IntegrationInfoBox info_box; + UpdateFlags update_flags = update_values | update_gradients; + info_box.add_update_flags_all(update_flags); + info_box.initialize(fe, mapping); + + MeshWorker::DoFInfo dof_info(dof_handler); + MeshWorker::Assembler::MatrixSimple assembler; + assembler.initialize(system_matrix); + MatrixIntegrator integrator; + MeshWorker::integration_loop(dof_handler.begin_active(), + dof_handler.end(), + dof_info, + info_box, + integrator, + assembler); + + system_matrix.compress(VectorOperation::add); + + for (unsigned int i = 0; i < system_rhs.size(); ++i) + system_rhs(i) = 0.01 * i - 0.000001 * i * i; +} + + + +template +void +Step4::solve() +{ + deallog.push(Utilities::int_to_string(dof_handler.n_dofs(), 5)); + TrilinosWrappers::PreconditionAMG preconditioner; + TrilinosWrappers::PreconditionAMG::AdditionalData data; + data.smoother_sweeps = 2; + { + solution = 0; + SolverControl solver_control(1000, 1e-10); + SolverCG<> solver(solver_control); + preconditioner.initialize(system_matrix, data); + solver.solve(system_matrix, solution, system_rhs, preconditioner); + } + deallog.pop(); +} + + + +template +void +Step4::run() +{ + for (unsigned int cycle = 0; cycle < 2; ++cycle) + { + if (cycle == 0) + make_grid(); + else + triangulation.refine_global(1); + + setup_system(); + solve(); + } +} + + +int +main(int argc, char **argv) +{ + initlog(); + + Utilities::MPI::MPI_InitFinalize mpi_initialization( + argc, argv, testing_max_num_threads()); + + try + { + Step4<2> test; + test.run(); + } + catch (std::exception &exc) + { + deallog << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + deallog << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/trilinos/precondition_amg_dgp_03.with_64bit_indices=off.output b/tests/trilinos/precondition_amg_dgp_03.with_64bit_indices=off.output new file mode 100644 index 0000000000..5399b9277c --- /dev/null +++ b/tests/trilinos/precondition_amg_dgp_03.with_64bit_indices=off.output @@ -0,0 +1,5 @@ + +DEAL:12288:cg::Starting value 1970.22 +DEAL:12288:cg::Convergence step 121 value 0 +DEAL:49152:cg::Starting value 179304. +DEAL:49152:cg::Convergence step 271 value 0 -- 2.39.5