From: kronbichler Date: Mon, 30 Jun 2014 10:23:34 +0000 (+0000) Subject: Commit bug in Trilinos precondition X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=a16e93b3183c660f7cfcc32ff4e1bb213ce61b10;p=dealii-svn.git Commit bug in Trilinos precondition git-svn-id: https://svn.dealii.org/trunk@33099 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index ba57ec5880..95df70da1b 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -148,6 +148,13 @@ inconvenience this causes.

Specific improvements

    +
  1. Fixed: TrilinosWrappers::PreconditionAMG did not read user-provided + constant modes (aka null space) when the null space dimension is one but not + just the trivial one vector. This is now fixed. +
    + (Martin Kronbichler, 2014/06/30) +
    2. +
  2. Simplified interfaces for FEEvaluation: Previously, the user had to select the appropriate kernel (FEEvaluation, FEEvaluationGeneral, FEEvaluationDGP, FEEvaluationGL) for the matrix-free evaluation diff --git a/deal.II/include/deal.II/lac/trilinos_precondition.h b/deal.II/include/deal.II/lac/trilinos_precondition.h index ba1d6ce23e..9ab30fb98f 100644 --- a/deal.II/include/deal.II/lac/trilinos_precondition.h +++ b/deal.II/include/deal.II/lac/trilinos_precondition.h @@ -1301,7 +1301,7 @@ namespace TrilinosWrappers const unsigned int n_cycles = 1, const bool w_cyle = false, const double aggregation_threshold = 1e-4, - const std::vector > &constant_modes = std::vector > (1), + const std::vector > &constant_modes = std::vector > (0), const unsigned int smoother_sweeps = 2, const unsigned int smoother_overlap = 0, const bool output_details = false, diff --git a/deal.II/source/lac/trilinos_precondition.cc b/deal.II/source/lac/trilinos_precondition.cc index a526b310d7..4c42b28e78 100644 --- a/deal.II/source/lac/trilinos_precondition.cc +++ b/deal.II/source/lac/trilinos_precondition.cc @@ -611,7 +611,7 @@ namespace TrilinosWrappers constant_modes_dimension); std::vector dummy (constant_modes_dimension); - if (constant_modes_dimension > 1) + if (constant_modes_dimension > 0) { const bool constant_modes_are_global = additional_data.constant_modes[0].size() == n_rows; diff --git a/tests/trilinos/precondition_amg_dgp.cc b/tests/trilinos/precondition_amg_dgp.cc new file mode 100644 index 0000000000..dacb279f96 --- /dev/null +++ b/tests/trilinos/precondition_amg_dgp.cc @@ -0,0 +1,254 @@ +// --------------------------------------------------------------------- +// $Id$ +// +// Copyright (C) 2013 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 at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + + +// solves a 2D Poisson equation for linear FE_DGP elements (SIP +// discretization) with AMG preconditioner + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + + +#include +#include + + +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_Q 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 (7); +} + + + +template +void Step4::setup_system () +{ + dof_handler.distribute_dofs (fe); + + CompressedSparsityPattern 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()); + + MappingQ1 mapping; + 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 +void Step4::solve () +{ + + // variant 1: solve with AMG + deallog.push(Utilities::int_to_string(dof_handler.n_dofs(),5)); + TrilinosWrappers::PreconditionAMG preconditioner; + TrilinosWrappers::PreconditionAMG::AdditionalData data; + DoFTools::extract_constant_modes(dof_handler, std::vector(1,true), + data.constant_modes); + 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) +{ + std::ofstream logfile("output"); + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv); + + 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.with_64bit_indices=off.output b/tests/trilinos/precondition_amg_dgp.with_64bit_indices=off.output new file mode 100644 index 0000000000..45b1d6cb98 --- /dev/null +++ b/tests/trilinos/precondition_amg_dgp.with_64bit_indices=off.output @@ -0,0 +1,5 @@ + +DEAL:16641:cg::Starting value 5039.73 +DEAL:16641:cg::Convergence step 15 value 0 +DEAL:66049:cg::Starting value 407234. +DEAL:66049:cg::Convergence step 17 value 0