From 66d5cc387ba9415c5e039e06fdd3f934a6d4aa50 Mon Sep 17 00:00:00 2001 From: kronbichler Date: Sat, 14 Dec 2013 16:10:22 +0000 Subject: [PATCH] Apply patch by Andrew Baker. git-svn-id: https://svn.dealii.org/trunk@32009 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/doc/news/changes.h | 6 + .../deal.II/lac/trilinos_precondition.h | 230 ++++++------- deal.II/source/lac/trilinos_precondition.cc | 52 ++- tests/trilinos/precondition_amg_smoother.cc | 302 ++++++++++++++++++ .../trilinos/precondition_amg_smoother.output | 9 + 5 files changed, 446 insertions(+), 153 deletions(-) create mode 100644 tests/trilinos/precondition_amg_smoother.cc create mode 100644 tests/trilinos/precondition_amg_smoother.output diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 1cdd7325ca..afa759f909 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -59,6 +59,12 @@ inconvenience this causes.

Specific improvements

    +
  1. New: It is now possible to select between different smoothers and coarse + solvers in the Trilinos AMG preconditioners by a string to the smoother's name. +
    + (Andrew Baker, 2013/12/14) +
    2. +
diff --git a/deal.II/include/deal.II/lac/trilinos_precondition.h b/deal.II/include/deal.II/lac/trilinos_precondition.h index 0f295a99e1..eb8747a34d 100644 --- a/deal.II/include/deal.II/lac/trilinos_precondition.h +++ b/deal.II/include/deal.II/lac/trilinos_precondition.h @@ -1293,10 +1293,8 @@ namespace TrilinosWrappers struct AdditionalData { /** - * Constructor. By default, we - * pretend to work on elliptic - * problems with linear finite - * elements on a scalar equation. + * Constructor. By default, we pretend to work on elliptic problems with + * linear finite elements on a scalar equation. */ AdditionalData (const bool elliptic = true, const bool higher_order_elements = false, @@ -1306,153 +1304,151 @@ namespace TrilinosWrappers const std::vector > &constant_modes = std::vector > (1), const unsigned int smoother_sweeps = 2, const unsigned int smoother_overlap = 0, - const bool output_details = false); + const bool output_details = false, + const char* smoother_type = "Chebyshev", + const char* coarse_type = "Amesos-KLU"); /** - * Determines whether the AMG - * preconditioner should be optimized - * for elliptic problems (ML option - * smoothed aggregation SA, using a - * Chebyshev smoother) or for - * non-elliptic problems (ML option - * non-symmetric smoothed aggregation - * NSSA, smoother is SSOR with + * Determines whether the AMG preconditioner should be optimized for + * elliptic problems (ML option smoothed aggregation SA, using a + * Chebyshev smoother) or for non-elliptic problems (ML option + * non-symmetric smoothed aggregation NSSA, smoother is SSOR with * underrelaxation). */ bool elliptic; /** - * Determines whether the matrix that - * the preconditioner is built upon - * is generated from linear or - * higher-order elements. + * Determines whether the matrix that the preconditioner is built upon + * is generated from linear or higher-order elements. */ bool higher_order_elements; /** - * Defines how many multigrid cycles - * should be performed by the + * Defines how many multigrid cycles should be performed by the * preconditioner. */ unsigned int n_cycles; /** - * Defines whether a w-cycle should be - * used instead of the standard setting - * of a v-cycle. + * Defines whether a w-cycle should be used instead of the standard + * setting of a v-cycle. */ bool w_cycle; /** - * This threshold tells the AMG setup - * how the coarsening should be - * performed. In the AMG used by ML, - * all points that strongly couple - * with the tentative coarse-level - * point form one aggregate. The term - * strong coupling is - * controlled by the variable - * aggregation_threshold, - * meaning that all elements that are - * not smaller than - * aggregation_threshold - * times the diagonal element do - * couple strongly. + * This threshold tells the AMG setup how the coarsening should be + * performed. In the AMG used by ML, all points that strongly couple + * with the tentative coarse-level point form one aggregate. The term + * strong coupling is controlled by the variable + * aggregation_threshold, meaning that all elements that are + * not smaller than aggregation_threshold times the diagonal + * element do couple strongly. */ double aggregation_threshold; /** - * Specifies the constant modes (near - * null space) of the matrix. This - * parameter tells AMG whether we - * work on a scalar equation (where - * the near null space only consists - * of ones) or on a vector-valued + * Specifies the constant modes (near null space) of the matrix. This + * parameter tells AMG whether we work on a scalar equation (where the + * near null space only consists of ones) or on a vector-valued * equation. */ std::vector > constant_modes; /** - * Determines how many sweeps of the - * smoother should be performed. When - * the flag elliptic is set - * to true, i.e., for - * elliptic or almost elliptic - * problems, the polynomial degree of - * the Chebyshev smoother is set to - * smoother_sweeps. The term - * sweeps refers to the number of - * matrix-vector products performed - * in the Chebyshev case. In the - * non-elliptic case, - * smoother_sweeps sets the - * number of SSOR relaxation sweeps - * for post-smoothing to be - * performed. + * Determines how many sweeps of the smoother should be performed. When + * the flag elliptic is set to true, i.e., for + * elliptic or almost elliptic problems, the polynomial degree of the + * Chebyshev smoother is set to smoother_sweeps. The term + * sweeps refers to the number of matrix-vector products performed in + * the Chebyshev case. In the non-elliptic case, + * smoother_sweeps sets the number of SSOR relaxation sweeps + * for post-smoothing to be performed. */ unsigned int smoother_sweeps; /** - * Determines the overlap in the - * SSOR/Chebyshev error smoother when - * run in parallel. + * Determines the overlap in the SSOR/Chebyshev error smoother when run + * in parallel. */ unsigned int smoother_overlap; /** - * If this flag is set to - * true, then internal - * information from the ML - * preconditioner is printed to - * screen. This can be useful when + * If this flag is set to true, then internal information from + * the ML preconditioner is printed to screen. This can be useful when * debugging the preconditioner. */ bool output_details; + + /** + * Determines which smoother to use for the AMG cycle. Possibilities + * for smoother_type are the following: + * "Aztec" + * "IFPACK" + * "Jacobi" + * "ML symmetric Gauss-Seidel" + * "symmetric Gauss-Seidel" + * "ML Gauss-Seidel" + * "Gauss-Seidel" + * "block Gauss-Seidel" + * "symmetric block Gauss-Seidel" + * "Chebyshev" + * "MLS" + * "Hiptmair" + * "Amesos-KLU" + * "Amesos-Superlu" + * "Amesos-UMFPACK" + * "Amesos-Superludist" + * "Amesos-MUMPS" + * "user-defined" + * "SuperLU" + * "IFPACK-Chebyshev" + * "self" + * "do-nothing" + * "IC" + * "ICT" + * "ILU" + * "ILUT" + * "Block Chebyshev" + * "IFPACK-Block Chebyshev" + */ + const char* smoother_type; + + /** + * Determines which solver to use on the coarsest level. The same + * settings as for the smoother type are possible. + */ + const char* coarse_type; }; + /** - * Let Trilinos compute a multilevel - * hierarchy for the solution of a - * linear system with the given - * matrix. The function uses the - * matrix format specified in - * TrilinosWrappers::SparseMatrix. + * Let Trilinos compute a multilevel hierarchy for the solution of a + * linear system with the given matrix. The function uses the matrix + * format specified in TrilinosWrappers::SparseMatrix. */ void initialize (const SparseMatrix &matrix, const AdditionalData &additional_data = AdditionalData()); /** - * Let Trilinos compute a multilevel - * hierarchy for the solution of a - * linear system with the given - * matrix. The function uses the - * matrix format specified in - * TrilinosWrappers::SparseMatrix. + * Let Trilinos compute a multilevel hierarchy for the solution of a + * linear system with the given matrix. The function uses the matrix + * format specified in TrilinosWrappers::SparseMatrix. * - * This function is similar to the one - * above, but allows the user to set - * all the options of the Trilinos ML - * preconditioner. In order to find out - * about all the options for ML, we - * refer to the ML - * user's guide. In particular, - * users need to follow the ML - * instructions in case a vector-valued - * problem ought to be solved. + * This function is similar to the one above, but allows the user to set + * all the options of the Trilinos ML preconditioner. In order to find out + * about all the options for ML, we refer to the ML user's + * guide. In particular, users need to follow the ML instructions in + * case a vector-valued problem ought to be solved. */ void initialize (const SparseMatrix &matrix, const Teuchos::ParameterList &ml_parameters); /** - * Let Trilinos compute a multilevel - * hierarchy for the solution of a - * linear system with the given - * matrix. This function takes a - * deal.ii matrix and copies the - * content into a Trilinos matrix, so - * the function can be considered - * rather inefficient. + * Let Trilinos compute a multilevel hierarchy for the solution of a + * linear system with the given matrix. This function takes a deal.ii + * matrix and copies the content into a Trilinos matrix, so the function + * can be considered rather inefficient. */ template void initialize (const ::dealii::SparseMatrix &deal_ii_sparse_matrix, @@ -1461,47 +1457,33 @@ namespace TrilinosWrappers const ::dealii::SparsityPattern *use_this_sparsity = 0); /** - * This function can be used for a - * faster recalculation of the - * preconditioner construction when - * the matrix entries underlying the - * preconditioner have changed, but - * the matrix sparsity pattern has - * remained the same. What this - * function does is taking the - * already generated coarsening - * structure, computing the AMG - * prolongation and restriction - * according to a smoothed - * aggregation strategy and then - * building the whole multilevel - * hiearchy. This function can be - * considerably faster than the - * initialize function, since the - * coarsening pattern is usually the - * most difficult thing to do when - * setting up the AMG ML - * preconditioner. + * This function can be used for a faster recalculation of the + * preconditioner construction when the matrix entries underlying the + * preconditioner have changed, but the matrix sparsity pattern has + * remained the same. What this function does is taking the already + * generated coarsening structure, computing the AMG prolongation and + * restriction according to a smoothed aggregation strategy and then + * building the whole multilevel hiearchy. This function can be + * considerably faster than the initialize function, since the coarsening + * pattern is usually the most difficult thing to do when setting up the + * AMG ML preconditioner. */ void reinit (); /** - * Destroys the preconditioner, leaving - * an object like just after having + * Destroys the preconditioner, leaving an object like just after having * called the constructor. */ void clear (); /** - * Prints an estimate of the memory - * consumption of this class. + * Prints an estimate of the memory consumption of this class. */ size_type memory_consumption () const; private: /** - * A copy of the deal.II matrix into - * Trilinos format. + * A copy of the deal.II matrix into Trilinos format. */ std_cxx1x::shared_ptr trilinos_matrix; }; diff --git a/deal.II/source/lac/trilinos_precondition.cc b/deal.II/source/lac/trilinos_precondition.cc index db67419c43..be35b10663 100644 --- a/deal.II/source/lac/trilinos_precondition.cc +++ b/deal.II/source/lac/trilinos_precondition.cc @@ -516,7 +516,9 @@ namespace TrilinosWrappers const std::vector > &constant_modes, const unsigned int smoother_sweeps, const unsigned int smoother_overlap, - const bool output_details) + const bool output_details, + const char* smoother_type, + const char* coarse_type) : elliptic (elliptic), higher_order_elements (higher_order_elements), @@ -526,7 +528,9 @@ namespace TrilinosWrappers constant_modes (constant_modes), smoother_sweeps (smoother_sweeps), smoother_overlap (smoother_overlap), - output_details (output_details) + output_details (output_details), + smoother_type (smoother_type), + coarse_type (coarse_type) {} @@ -543,26 +547,17 @@ namespace TrilinosWrappers if (additional_data.elliptic == true) { ML_Epetra::SetDefaults("SA",parameter_list); - parameter_list.set("smoother: type", "Chebyshev"); - - // uncoupled mode can give a lot of - // warnings or even fail when there - // are too many entries per row and - // aggreggation gets complicated, but - // MIS does not work if too few - // elements are located on one - // processor. work around these - // warnings by choosing the different - // strategies in different - // situations: for low order, always - // use the standard choice - // uncoupled. if higher order, right - // now we also just use Uncoupled, - // but we should be aware that maybe - // MIS might be needed + + // uncoupled mode can give a lot of warnings or even fail when there + // are too many entries per row and aggreggation gets complicated, but + // MIS does not work if too few elements are located on one + // processor. work around these warnings by choosing the different + // strategies in different situations: for low order, always use the + // standard choice uncoupled. if higher order, right now we also just + // use Uncoupled, but we should be aware that maybe MIS might be + // needed // - // TODO: Maybe there are any - // other/better options? + // TODO: Maybe there are any other/better options? if (additional_data.higher_order_elements) { //if (matrix.m()/matrix.matrix->Comm().NumProc() < 50000) @@ -578,6 +573,9 @@ namespace TrilinosWrappers parameter_list.set("aggregation: block scaling", true); } + parameter_list.set("smoother: type", additional_data.smoother_type); + parameter_list.set("coarse: type", additional_data.coarse_type); + parameter_list.set("smoother: sweeps", static_cast(additional_data.smoother_sweeps)); parameter_list.set("cycle applications", @@ -622,10 +620,8 @@ namespace TrilinosWrappers ExcDimensionMismatch(n_rows, global_length(distributed_constant_modes))); - // Reshape null space as a - // contiguous vector of - // doubles so that Trilinos - // can read from it. + // Reshape null space as a contiguous vector of doubles so that + // Trilinos can read from it. for (size_type d=0; d 0) parameter_list.set("null space: vectors", distributed_constant_modes.Values()); - // 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 + // 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 else parameter_list.set("null space: vectors", &dummy[0]); diff --git a/tests/trilinos/precondition_amg_smoother.cc b/tests/trilinos/precondition_amg_smoother.cc new file mode 100644 index 0000000000..11c1754858 --- /dev/null +++ b/tests/trilinos/precondition_amg_smoother.cc @@ -0,0 +1,302 @@ +// --------------------------------------------------------------------- +// $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 elements with AMG preconditioner +// and various smoothers + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + + +template +class Step4 +{ +public: + Step4 (); + void run (); + +private: + void make_grid (); + void setup_system(); + void assemble_system (); + void solve (); + + Triangulation triangulation; + FE_Q fe; + DoFHandler dof_handler; + + ConstraintMatrix constraints; + + TrilinosWrappers::SparseMatrix system_matrix; + + Vector solution; + Vector system_rhs; +}; + + +template +class RightHandSide : public Function +{ +public: + RightHandSide () : Function() {} + + virtual double value (const Point &p, + const unsigned int component = 0) const; +}; + + + +template +class BoundaryValues : public Function +{ +public: + BoundaryValues () : Function() {} + + virtual double value (const Point &p, + const unsigned int component = 0) const; +}; + + + + +template +double RightHandSide::value (const Point &p, + const unsigned int /*component*/) const +{ + double return_value = 0; + for (unsigned int i=0; i +double BoundaryValues::value (const Point &p, + const unsigned int /*component*/) const +{ + return p.square(); +} + + + +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); + + constraints.clear(); + std::map boundary_values; + VectorTools::interpolate_boundary_values (dof_handler, + 0, + BoundaryValues(), + constraints); + constraints.close(); + + CompressedSparsityPattern c_sparsity(dof_handler.n_dofs()); + DoFTools::make_sparsity_pattern (dof_handler, c_sparsity, constraints, false); + system_matrix.reinit (c_sparsity); + + solution.reinit (dof_handler.n_dofs()); + system_rhs.reinit (dof_handler.n_dofs()); +} + + +template +void Step4::assemble_system () +{ + QGauss quadrature_formula(fe.degree+1); + + const RightHandSide right_hand_side; + + 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); + + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + + for (; cell!=endc; ++cell) + { + fe_values.reinit (cell); + cell_matrix = 0; + cell_rhs = 0; + + for (unsigned int q_point=0; q_pointget_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); +} + + + +template +void Step4::solve () +{ + + // variant 1: solve with AMG + deallog.push(Utilities::int_to_string(dof_handler.n_dofs(),5)); + deallog.push("Chebyshev"); + TrilinosWrappers::PreconditionAMG preconditioner; + TrilinosWrappers::PreconditionAMG::AdditionalData data; + data.coarse_type = "Amesos-KLU"; + data.smoother_type = "Chebyshev"; + data.output_details = true; + data.aggregation_threshold = 1e-3; + data.smoother_sweeps = 3; + { + 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(); + + deallog.push("SGS"); + data.smoother_type = "symmetric Gauss-Seidel"; + data.smoother_sweeps = 2; + { + solution = 0; + SolverControl solver_control (1000, 1e-12); + SolverCG<> solver (solver_control); + preconditioner.initialize(system_matrix, data); + solver.solve (system_matrix, solution, system_rhs, + preconditioner); + } + deallog.pop(); + 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(); + assemble_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_smoother.output b/tests/trilinos/precondition_amg_smoother.output new file mode 100644 index 0000000000..ba3157c883 --- /dev/null +++ b/tests/trilinos/precondition_amg_smoother.output @@ -0,0 +1,9 @@ + +DEAL:04225:Chebyshev:cg::Starting value 21.8299 +DEAL:04225:Chebyshev:cg::Convergence step 10 value 0 +DEAL:04225:SGS:cg::Starting value 21.8299 +DEAL:04225:SGS:cg::Convergence step 10 value 0 +DEAL:16641:Chebyshev:cg::Starting value 30.8770 +DEAL:16641:Chebyshev:cg::Convergence step 8 value 0 +DEAL:16641:SGS:cg::Starting value 30.8770 +DEAL:16641:SGS:cg::Convergence step 7 value 0 -- 2.39.5