<h3>Specific improvements</h3>
<ol>
+ <li> 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.
+ <br>
+ (Andrew Baker, 2013/12/14)
+ </li>
+
</ol>
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,
const std::vector<std::vector<bool> > &constant_modes = std::vector<std::vector<bool> > (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
- * <em>strong coupling</em> is
- * controlled by the variable
- * <tt>aggregation_threshold</tt>,
- * meaning that all elements that are
- * not smaller than
- * <tt>aggregation_threshold</tt>
- * 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
+ * <em>strong coupling</em> is controlled by the variable
+ * <tt>aggregation_threshold</tt>, meaning that all elements that are
+ * not smaller than <tt>aggregation_threshold</tt> 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<std::vector<bool> > constant_modes;
/**
- * Determines how many sweeps of the
- * smoother should be performed. When
- * the flag <tt>elliptic</tt> is set
- * to <tt>true</tt>, i.e., for
- * elliptic or almost elliptic
- * problems, the polynomial degree of
- * the Chebyshev smoother is set to
- * <tt>smoother_sweeps</tt>. The term
- * sweeps refers to the number of
- * matrix-vector products performed
- * in the Chebyshev case. In the
- * non-elliptic case,
- * <tt>smoother_sweeps</tt> 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 <tt>elliptic</tt> is set to <tt>true</tt>, i.e., for
+ * elliptic or almost elliptic problems, the polynomial degree of the
+ * Chebyshev smoother is set to <tt>smoother_sweeps</tt>. The term
+ * sweeps refers to the number of matrix-vector products performed in
+ * the Chebyshev case. In the non-elliptic case,
+ * <tt>smoother_sweeps</tt> 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
- * <tt>true</tt>, then internal
- * information from the ML
- * preconditioner is printed to
- * screen. This can be useful when
+ * If this flag is set to <tt>true</tt>, 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 <a
- * href=http://trilinos.sandia.gov/packages/ml/mlguide5.pdf>ML
- * user's guide</a>. 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 <a
+ * href=http://trilinos.sandia.gov/packages/ml/mlguide5.pdf>ML user's
+ * guide</a>. 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 <typename number>
void initialize (const ::dealii::SparseMatrix<number> &deal_ii_sparse_matrix,
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<SparseMatrix> trilinos_matrix;
};
const std::vector<std::vector<bool> > &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),
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)
{}
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)
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<int>(additional_data.smoother_sweeps));
parameter_list.set("cycle applications",
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<constant_modes_dimension; ++d)
for (size_type row=0; row<my_size; ++row)
{
if (my_size > 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]);
--- /dev/null
+// ---------------------------------------------------------------------
+// $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 <deal.II/lac/trilinos_sparse_matrix.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/lac/compressed_sparsity_pattern.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/lac/trilinos_precondition.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/base/function.h>
+#include <deal.II/grid/tria.h>
+
+#include <fstream>
+#include <iomanip>
+
+
+template <int dim>
+class Step4
+{
+public:
+ Step4 ();
+ void run ();
+
+private:
+ void make_grid ();
+ void setup_system();
+ void assemble_system ();
+ void solve ();
+
+ Triangulation<dim> triangulation;
+ FE_Q<dim> fe;
+ DoFHandler<dim> dof_handler;
+
+ ConstraintMatrix constraints;
+
+ TrilinosWrappers::SparseMatrix system_matrix;
+
+ Vector<double> solution;
+ Vector<double> system_rhs;
+};
+
+
+template <int dim>
+class RightHandSide : public Function<dim>
+{
+public:
+ RightHandSide () : Function<dim>() {}
+
+ virtual double value (const Point<dim> &p,
+ const unsigned int component = 0) const;
+};
+
+
+
+template <int dim>
+class BoundaryValues : public Function<dim>
+{
+public:
+ BoundaryValues () : Function<dim>() {}
+
+ virtual double value (const Point<dim> &p,
+ const unsigned int component = 0) const;
+};
+
+
+
+
+template <int dim>
+double RightHandSide<dim>::value (const Point<dim> &p,
+ const unsigned int /*component*/) const
+{
+ double return_value = 0;
+ for (unsigned int i=0; i<dim; ++i)
+ return_value += 4*std::pow(p(i), 4);
+
+ return return_value;
+}
+
+
+
+template <int dim>
+double BoundaryValues<dim>::value (const Point<dim> &p,
+ const unsigned int /*component*/) const
+{
+ return p.square();
+}
+
+
+
+template <int dim>
+Step4<dim>::Step4 ()
+ :
+ fe (1),
+ dof_handler (triangulation)
+{}
+
+
+template <int dim>
+void Step4<dim>::make_grid ()
+{
+ GridGenerator::hyper_cube (triangulation, -1, 1);
+ triangulation.refine_global (6);
+}
+
+
+
+template <int dim>
+void Step4<dim>::setup_system ()
+{
+ dof_handler.distribute_dofs (fe);
+
+ constraints.clear();
+ std::map<unsigned int,double> boundary_values;
+ VectorTools::interpolate_boundary_values (dof_handler,
+ 0,
+ BoundaryValues<dim>(),
+ 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 <int dim>
+void Step4<dim>::assemble_system ()
+{
+ QGauss<dim> quadrature_formula(fe.degree+1);
+
+ const RightHandSide<dim> right_hand_side;
+
+ FEValues<dim> 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<double> cell_matrix (dofs_per_cell, dofs_per_cell);
+ Vector<double> cell_rhs (dofs_per_cell);
+
+ std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
+
+ typename DoFHandler<dim>::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_point<n_q_points; ++q_point)
+ 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) += (fe_values.shape_value (i, q_point) *
+ right_hand_side.value (fe_values.quadrature_point (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);
+}
+
+
+
+template <int dim>
+void Step4<dim>::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 <int dim>
+void Step4<dim>::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;
+ };
+}
--- /dev/null
+
+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
projects / dealii-svn.git / commitdiff