]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add MueLu preconditioner using MueLu::MLParameterListInterpreter. 214/head
authorBruno Turcksin <bruno.turcksin@gmail.com>
Mon, 6 Oct 2014 14:31:30 +0000 (09:31 -0500)
committerBruno Turcksin <bruno.turcksin@gmail.com>
Wed, 11 Mar 2015 21:35:36 +0000 (16:35 -0500)
14 files changed:
cmake/configure/configure_2_trilinos.cmake
doc/news/changes.h
include/deal.II/lac/trilinos_precondition.h
source/lac/trilinos_precondition.cc
tests/mpi/muelu_periodicity.cc [new file with mode: 0644]
tests/mpi/muelu_periodicity.with_trilinos=true.mpirun=3.output [new file with mode: 0644]
tests/mpi/muelu_periodicity.with_trilinos=true.mpirun=5.output [new file with mode: 0644]
tests/mpi/muelu_periodicity.with_trilinos=true.mpirun=7.output [new file with mode: 0644]
tests/trilinos/precondition_muelu_dgp.cc [new file with mode: 0644]
tests/trilinos/precondition_muelu_dgp.with_64bit_indices=off.output [new file with mode: 0644]
tests/trilinos/precondition_muelu_q_iso_q1.cc [new file with mode: 0644]
tests/trilinos/precondition_muelu_q_iso_q1.with_64bit_indices=off.output [new file with mode: 0644]
tests/trilinos/precondition_muelu_smoother.cc [new file with mode: 0644]
tests/trilinos/precondition_muelu_smoother.with_64bit_indices=off.output [new file with mode: 0644]

index 05fcc6872ca9360d6641cd2410dcced303eeeae0..1e251b6d7246a5180faedae15c4eaa21f7418417 100644 (file)
@@ -54,6 +54,20 @@ MACRO(FEATURE_TRILINOS_FIND_EXTERNAL var)
       ENDIF()
     ENDFOREACH()
 
+    IF((TRILINOS_VERSION_MAJOR EQUAL 11 AND
+        NOT (TRILINOS_VERSION_MINOR LESS 14)) 
+       OR
+       (NOT (TRILINOS_VERSION_MAJOR LESS 12)))
+        ITEM_MATCHES(_module_found MueLu ${Trilinos_PACKAGE_LIST})
+      IF(_module_found)
+        MESSAGE(STATUS "Found MueLu")
+      ELSE()
+        MESSAGE(STATUS "Module MueLu not found!")
+        SET(_modules_missing "${_modules_missing} MueLu")
+        SET(${var} FALSE)
+      ENDIF()
+    ENDIF()
+
     IF(NOT ${var})
       MESSAGE(STATUS "Could not find a sufficient Trilinos installation: "
         "Missing ${_modules_missing}"
index 3e672fb1baf6115e3cdc6af9b47852b7e73cf66f..909ba4aa0891d7012c2bd368c6b584a0c13a8733 100644 (file)
@@ -376,6 +376,13 @@ inconvenience this causes.
 <h3>Specific improvements</h3>
 
 <ol>
+  <li> New: Add MueLu preconditioner from Trilinos. This is a new algebraic 
+  multigrid package. The input parameters are almost the same as the ones 
+  from ML so that the two preconditioners can be easily swapped.
+  <br>
+  (Bruno Turcksin, 2015/03/11)
+  </li>
+
   <li> Fixed: Iterating over the elements of a TrilinosWrappers::SparseMatrix
   object previously led to errors if the matrix was in fact stored in
   parallel across multiple MPI processes. This is now fixed: rows not
index 079c2dd55114fb22041e729b7d13322719ead8d8..eeef76c23baf9f2bd7b3d05ccbd8da277ac38cdc 100644 (file)
@@ -1,6 +1,6 @@
 // ---------------------------------------------------------------------
 //
-// Copyright (C) 2008 - 2014 by the deal.II authors
+// Copyright (C) 2008 - 2015 by the deal.II authors
 //
 // This file is part of the deal.II library.
 //
 #  include <Epetra_RowMatrix.h>
 #  include <Epetra_Vector.h>
 
+#  if DEAL_II_TRILINOS_VERSION_GTE(11,14,0)
+#  include <MueLu.hpp>
+#  include <Teuchos_RCP.hpp>
+#  endif
+
 // forward declarations
 class Ifpack_Preconditioner;
 class Ifpack_Chebyshev;
@@ -1541,6 +1546,235 @@ namespace TrilinosWrappers
 
 
 
+#if DEAL_II_TRILINOS_VERSION_GTE(11,14,0)
+  /**
+   * This class implements an algebraic multigrid (AMG) preconditioner based
+   * on the Trilinos MueLu implementation, which is a black-box preconditioner
+   * that works well for many PDE-based linear problems. The interface of
+   * PreconditionerAMGMueLu is the same as the interface of PreconditionerAMG
+   * except for the higher_order_elements parameter which does not exist in
+   * PreconditionerAMGMueLu.
+   *
+   * This class requires Trilinos 11.14 or higher.
+   *
+   * This interface should not be considered as stable.
+   *
+   * @ingroup TrilinosWrappers
+   * @ingroup Preconditioners
+   * @author Bruno Turcksin, 2014
+   */
+  class PreconditionAMGMueLu : public PreconditionBase
+  {
+  public:
+
+
+    /**
+     * A data structure that is used to control details of how the algebraic
+     * multigrid is set up. The flags detailed in here are then passed to
+     * the Trilinos MueLu implementation. A structure of the current type are
+     * passed to the constructor of PreconditionAMGMueLu.
+     */
+    struct AdditionalData
+    {
+      /**
+       * Constructor. By default, we pretend to work on elliptic problems with
+       * linear finite elements on a scalar equation.
+       */
+      AdditionalData (const bool                             elliptic = true,
+                      const unsigned int                     n_cycles = 1,
+                      const bool                             w_cyle = false,
+                      const double                           aggregation_threshold = 1e-4,
+                      const std::vector<std::vector<bool> > &constant_modes = std::vector<std::vector<bool> > (0),
+                      const unsigned int                     smoother_sweeps = 2,
+                      const unsigned int                     smoother_overlap = 0,
+                      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 (MueLu option smoothed aggregation SA, using a
+       * Chebyshev smoother) or for non-elliptic problems (MueLu option
+       * non-symmetric smoothed aggregation NSSA, smoother is SSOR with
+       * underrelaxation).
+       */
+      bool elliptic;
+
+      /**
+       * 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.
+       */
+      bool w_cycle;
+
+      /**
+       * This threshold tells the AMG setup how the coarsening should be
+       * performed. In the AMG used by MueLu, 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
+       * 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.
+       */
+      unsigned int smoother_sweeps;
+
+      /**
+       * 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
+       * debugging the preconditioner.
+       */
+      bool output_details;
+
+      /**
+       * Determines which smoother to use for the AMG cycle. Possibilities
+       * for smoother_type are the following:
+       * <ul>
+       *   <li>  "Aztec" </li>
+       *   <li>  "IFPACK" </li>
+       *   <li>  "Jacobi" </li>
+       *   <li>  "ML symmetric Gauss-Seidel" </li>
+       *   <li>  "symmetric Gauss-Seidel" </li>
+       *   <li>  "ML Gauss-Seidel" </li>
+       *   <li>  "Gauss-Seidel" </li>
+       *   <li>  "block Gauss-Seidel" </li>
+       *   <li>  "symmetric block Gauss-Seidel" </li>
+       *   <li>  "Chebyshev" </li>
+       *   <li>  "MLS" </li>
+       *   <li>  "Hiptmair" </li>
+       *   <li>  "Amesos-KLU" </li>
+       *   <li>  "Amesos-Superlu" </li>
+       *   <li>  "Amesos-UMFPACK" </li>
+       *   <li>  "Amesos-Superludist" </li>
+       *   <li>  "Amesos-MUMPS" </li>
+       *   <li>  "user-defined" </li>
+       *   <li>  "SuperLU" </li>
+       *   <li>  "IFPACK-Chebyshev" </li>
+       *   <li>  "self" </li>
+       *   <li>  "do-nothing" </li>
+       *   <li>  "IC" </li>
+       *   <li>  "ICT" </li>
+       *   <li>  "ILU" </li>
+       *   <li>  "ILUT" </li>
+       *   <li>  "Block Chebyshev" </li>
+       *   <li>  "IFPACK-Block Chebyshev" </li>
+       * </ul>
+       */
+      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;
+    };
+
+    /**
+     * Destructor.
+     */
+    ~PreconditionAMGMueLu();
+
+    /**
+     * 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. As opposed to the other initialize
+     * function above, this function uses an object of type Epetra_CrsMatrixCrs.
+     */
+    void initialize (const Epetra_CrsMatrix &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.
+     *
+     * This function is similar to the one above, but allows the user to set
+     * most of 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>. Not all ML options have a corresponding MueLu option.
+     */
+    void initialize (const SparseMatrix   &matrix,
+                     Teuchos::ParameterList &muelu_parameters);
+
+    /**
+     * Let Trilinos compute a multilevel hierarchy for the solution of a
+     * linear system with the given matrix. As opposed to the other initialize
+     * function above, this function uses an object of type Epetra_CrsMatrix.
+     */
+    void initialize (const Epetra_CrsMatrix &matrix,
+                     Teuchos::ParameterList &muelu_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.
+     */
+    template <typename number>
+    void initialize (const ::dealii::SparseMatrix<number> &deal_ii_sparse_matrix,
+                     const AdditionalData                 &additional_data = AdditionalData(),
+                     const double                          drop_tolerance = 1e-13,
+                     const ::dealii::SparsityPattern      *use_this_sparsity = 0);
+
+    /**
+     * 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.
+     */
+    size_type memory_consumption () const;
+
+  private:
+    /**
+     * A copy of the deal.II matrix into Trilinos format.
+     */
+    std_cxx11::shared_ptr<SparseMatrix> trilinos_matrix;
+  };
+#endif
+
+
+
   /**
    * A wrapper class for an identity preconditioner for Trilinos matrices.
    *
index f3e223df81a24fe16331aa031c8aedb141bdaee6..8012fda8ff6aee6466c83b5ea769b1b906a9811c 100644 (file)
 #  include <ml_include.h>
 #  include <ml_MultiLevelPreconditioner.h>
 
+#if DEAL_II_TRILINOS_VERSION_GTE(11,14,0)
+#include <MueLu_EpetraOperator.hpp>
+#include <MueLu_MLParameterListInterpreter.hpp>
+#endif
 
 DEAL_II_NAMESPACE_OPEN
 
@@ -955,8 +959,258 @@ namespace TrilinosWrappers
 
 
 
+#if DEAL_II_TRILINOS_VERSION_GTE(11,14,0)
+  /* -------------------------- PreconditionAMGMueLu --------------------- */
 
-  /* -------------------------- PreconditionAMG -------------------------- */
+  PreconditionAMGMueLu::AdditionalData::
+  AdditionalData (const bool                             elliptic,
+                  const unsigned int                     n_cycles,
+                  const bool                             w_cycle,
+                  const double                           aggregation_threshold,
+                  const std::vector<std::vector<bool> > &constant_modes,
+                  const unsigned int                     smoother_sweeps,
+                  const unsigned int                     smoother_overlap,
+                  const bool                             output_details,
+                  const char                            *smoother_type,
+                  const char                            *coarse_type)
+    :
+    elliptic (elliptic),
+    n_cycles (n_cycles),
+    w_cycle (w_cycle),
+    aggregation_threshold (aggregation_threshold),
+    constant_modes (constant_modes),
+    smoother_sweeps (smoother_sweeps),
+    smoother_overlap (smoother_overlap),
+    output_details (output_details),
+    smoother_type (smoother_type),
+    coarse_type (coarse_type)
+  {}
+
+
+  PreconditionAMGMueLu::~PreconditionAMGMueLu()
+  {
+    preconditioner.reset();
+    trilinos_matrix.reset();
+  }
+
+
+
+  void
+  PreconditionAMGMueLu::initialize (const SparseMatrix   &matrix,
+                                    const AdditionalData &additional_data)
+  {
+    initialize(matrix.trilinos_matrix(), additional_data);
+  }
+
+
+
+  void
+  PreconditionAMGMueLu::initialize (const Epetra_CrsMatrix &matrix,
+                                    const AdditionalData   &additional_data)
+  {
+    // Build the AMG preconditioner.
+    Teuchos::ParameterList parameter_list;
+
+    if (additional_data.elliptic == true)
+      ML_Epetra::SetDefaults("SA",parameter_list);
+    else
+      {
+        ML_Epetra::SetDefaults("NSSA",parameter_list);
+        parameter_list.set("aggregation: block scaling", true);
+      }
+    // MIS does not exist anymore, only choice are uncoupled and coupled. When using
+    // uncoupled, aggregates cannot span multiple processes. When using coupled
+    // aggregates can span multiple processes.
+    parameter_list.set("aggregation: type", "Uncoupled");
+
+    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",
+                       static_cast<int>(additional_data.n_cycles));
+    if (additional_data.w_cycle == true)
+      parameter_list.set("prec type", "MGW");
+    else
+      parameter_list.set("prec type", "MGV");
+
+    parameter_list.set("smoother: Chebyshev alpha",10.);
+    parameter_list.set("smoother: ifpack overlap",
+                       static_cast<int>(additional_data.smoother_overlap));
+    parameter_list.set("aggregation: threshold",
+                       additional_data.aggregation_threshold);
+    parameter_list.set("coarse: max size", 2000);
+
+    if (additional_data.output_details)
+      parameter_list.set("ML output", 10);
+    else
+      parameter_list.set("ML output", 0);
+
+    const Epetra_Map &domain_map = matrix.OperatorDomainMap();
+
+    const size_type constant_modes_dimension =
+      additional_data.constant_modes.size();
+    Epetra_MultiVector distributed_constant_modes (domain_map,
+                                                   constant_modes_dimension > 0 ?
+                                                   constant_modes_dimension : 1);
+    std::vector<double> dummy (constant_modes_dimension);
+
+    if (constant_modes_dimension > 0)
+      {
+        const size_type n_rows = n_global_rows(matrix);
+        const bool constant_modes_are_global =
+          additional_data.constant_modes[0].size() == n_rows;
+        const size_type n_relevant_rows =
+          constant_modes_are_global ? n_rows : additional_data.constant_modes[0].size();
+        const size_type my_size = domain_map.NumMyElements();
+        if (constant_modes_are_global == false)
+          Assert (n_relevant_rows == my_size,
+                  ExcDimensionMismatch(n_relevant_rows, my_size));
+        Assert (n_rows ==
+                static_cast<size_type>(global_length(distributed_constant_modes)),
+                ExcDimensionMismatch(n_rows,
+                                     global_length(distributed_constant_modes)));
+
+        // 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)
+            {
+              TrilinosWrappers::types::int_type global_row_id =
+                constant_modes_are_global ? gid(domain_map,row) : row;
+              distributed_constant_modes[d][row] =
+                additional_data.constant_modes[d][global_row_id];
+            }
+
+        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());
+        // 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]);
+      }
+
+    initialize (matrix, parameter_list);
+  }
+
+
+
+  void
+  PreconditionAMGMueLu::initialize (const SparseMatrix           &matrix,
+                                    Teuchos::ParameterList &muelu_parameters)
+  {
+    initialize(matrix.trilinos_matrix(), muelu_parameters);
+  }
+
+
+
+  void
+  PreconditionAMGMueLu::initialize (const Epetra_CrsMatrix       &matrix,
+                                    Teuchos::ParameterList &muelu_parameters)
+  {
+    // We cannot use MueLu::CreateEpetraOperator directly because, we cannot
+    // transfer ownership of MueLu::EpetraOperator from Teuchos::RCP to
+    // std::shared_ptr.
+
+    // For now, just use serial node, i.e. no multithreaing or GPU.
+    typedef KokkosClassic::DefaultNode::DefaultNodeType node;
+    preconditioner.reset ();
+
+    // Cast matrix into a MueLu::Matrix. The constness needs to be cast away.
+    // MueLu uses Teuchos::RCP which are Trilinos version of std::shared_ptr.
+    Teuchos::RCP<Epetra_CrsMatrix> rcp_matrix = Teuchos::rcpFromRef(
+                                                  *(const_cast<Epetra_CrsMatrix *>(&matrix)));
+    Teuchos::RCP<Xpetra::CrsMatrix<double,int,int,node> > muelu_crs_matrix =
+      Teuchos::rcp(new Xpetra::EpetraCrsMatrix (rcp_matrix));
+    Teuchos::RCP<Xpetra::Matrix<double,int,int,node> > muelu_matrix =
+      Teuchos::rcp(new Xpetra::CrsMatrixWrap<double,int,int,node> (muelu_crs_matrix));
+
+    // Create the multigrid hierarchy using ML parameters.
+    Teuchos::RCP<MueLu::HierarchyManager<double,int,int,node> > hierarchy_factory;
+    hierarchy_factory = Teuchos::rcp(
+                          new MueLu::MLParameterListInterpreter<double,int,int,node> (muelu_parameters));
+    Teuchos::RCP<MueLu::Hierarchy<double,int,int,node> > hierarchy = hierarchy_factory->CreateHierarchy();
+    hierarchy->GetLevel(0)->Set("A",muelu_matrix);
+    hierarchy_factory->SetupHierarchy(*hierarchy);
+
+    // MueLu::EpetraOperator is just a wrapper around a "standard"
+    // Epetra_Operator.
+    preconditioner.reset(new MueLu::EpetraOperator(hierarchy));
+  }
+
+
+
+  template <typename number>
+  void
+  PreconditionAMGMueLu::
+  initialize (const ::dealii::SparseMatrix<number> &deal_ii_sparse_matrix,
+              const AdditionalData                 &additional_data,
+              const double                          drop_tolerance,
+              const ::dealii::SparsityPattern      *use_this_sparsity)
+  {
+    preconditioner.reset();
+    const size_type n_rows = deal_ii_sparse_matrix.m();
+
+    // Init Epetra Matrix using an
+    // equidistributed map; avoid
+    // storing the nonzero
+    // elements.
+    vector_distributor.reset (new Epetra_Map(static_cast<TrilinosWrappers::types::int_type>(n_rows),
+                                             0, communicator));
+
+    if (trilinos_matrix.get() == 0)
+      trilinos_matrix.reset (new SparseMatrix());
+
+    trilinos_matrix->reinit (*vector_distributor, *vector_distributor,
+                             deal_ii_sparse_matrix, drop_tolerance, true,
+                             use_this_sparsity);
+
+    initialize (*trilinos_matrix, additional_data);
+  }
+
+
+
+  void PreconditionAMGMueLu::clear ()
+  {
+    PreconditionBase::clear();
+    trilinos_matrix.reset();
+  }
+
+
+
+  PreconditionAMGMueLu::size_type
+  PreconditionAMGMueLu::memory_consumption() const
+  {
+    unsigned int memory = sizeof(this);
+
+    // todo: find a way to read out ML's data
+    // sizes
+    if (trilinos_matrix.get() != 0)
+      memory += trilinos_matrix->memory_consumption();
+    return memory;
+  }
+
+
+
+
+  // explicit instantiations
+  template void PreconditionAMGMueLu::initialize (const ::dealii::SparseMatrix<double> &,
+                                                  const AdditionalData &, const double,
+                                                  const ::dealii::SparsityPattern *);
+  template void PreconditionAMGMueLu::initialize (const ::dealii::SparseMatrix<float> &,
+                                                  const AdditionalData &, const double,
+                                                  const ::dealii::SparsityPattern *);
+#endif
+
+
+
+  /* -------------------------- PreconditionIdentity --------------------- */
 
   void
   PreconditionIdentity::vmult(VectorBase       &dst,
diff --git a/tests/mpi/muelu_periodicity.cc b/tests/mpi/muelu_periodicity.cc
new file mode 100644 (file)
index 0000000..18b0ef6
--- /dev/null
@@ -0,0 +1,819 @@
+/* ---------------------------------------------------------------------
+ *
+ * Copyright (C) 2008 - 2014 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.
+ *
+ * ---------------------------------------------------------------------
+ */
+
+//
+// Check MueLue in parallel. This test is the same as periodicity_02 except that
+// the preconditioner is MueLu.
+//
+
+#define PERIODIC
+
+#include "../tests.h"
+#include <deal.II/base/conditional_ostream.h>
+
+#include <deal.II/distributed/grid_refinement.h>
+
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/constraint_matrix.h>
+
+#include <deal.II/lac/trilinos_solver.h>
+#include <deal.II/lac/trilinos_precondition.h>
+#include <deal.II/lac/trilinos_block_sparse_matrix.h>
+#include <deal.II/lac/trilinos_parallel_block_vector.h>
+#include <deal.II/lac/block_sparsity_pattern.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/grid/grid_tools.h>
+
+#include <deal.II/dofs/dof_renumbering.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/error_estimator.h>
+
+namespace Step22
+{
+  using namespace dealii;
+
+  template <int dim>
+  class StokesProblem
+  {
+  public:
+    StokesProblem (const unsigned int degree);
+    void run ();
+
+  private:
+    void setup_dofs ();
+    void assemble_system ();
+    void solve ();
+    void get_point_value (const Point<dim> point, const int proc,
+                          Vector<double> &value) const;
+    void check_periodicity(const unsigned int cycle) const;
+    void output_results (const unsigned int refinement_cycle) const;
+    void refine_mesh ();
+
+    const unsigned int   degree;
+
+    MPI_Comm                                    mpi_communicator;
+
+    HyperShellBoundary<dim>                     boundary;
+    parallel::distributed::Triangulation<dim>   triangulation;
+    FESystem<dim>                               fe;
+    DoFHandler<dim>                             dof_handler;
+
+    ConstraintMatrix                            constraints;
+    std::vector<IndexSet>                       owned_partitioning;
+    std::vector<IndexSet>                       relevant_partitioning;
+
+    TrilinosWrappers::BlockSparsityPattern      sparsity_pattern;
+    TrilinosWrappers::BlockSparseMatrix         system_matrix;
+
+    TrilinosWrappers::MPI::BlockVector          solution;
+    TrilinosWrappers::MPI::BlockVector          system_rhs;
+
+    ConditionalOStream                          pcout;
+  };
+
+
+
+  template <int dim>
+  class BoundaryValues : public Function<dim>
+  {
+  public:
+    BoundaryValues () : Function<dim>(dim+1) {}
+
+    virtual double value (const Point<dim>   &p,
+                          const unsigned int  component = 0) const;
+
+    virtual void vector_value (const Point<dim> &p,
+                               Vector<double>   &value) const;
+  };
+
+
+  template <int dim>
+  double
+  BoundaryValues<dim>::value (const Point<dim>  &p,
+                              const unsigned int component) const
+  {
+    Assert (component < this->n_components,
+            ExcIndexRange (component, 0, this->n_components));
+
+//    if (component == 0)
+//      return (p[1]-.5)*(1.-p[1]);
+//    return 0;
+    return (1-2*(component==0))*p[(component+1)%2]/(p[0]*p[0]+p[1]*p[1]);
+  }
+
+
+  template <int dim>
+  void
+  BoundaryValues<dim>::vector_value (const Point<dim> &p,
+                                     Vector<double>   &values) const
+  {
+    for (unsigned int c=0; c<this->n_components; ++c)
+      values(c) = BoundaryValues<dim>::value (p, c);
+  }
+
+
+
+  template <int dim>
+  class RightHandSide : public Function<dim>
+  {
+  public:
+    RightHandSide () : Function<dim>(dim+1) {}
+
+    virtual double value (const Point<dim>   &p,
+                          const unsigned int  component = 0) const;
+
+    virtual void vector_value (const Point<dim> &p,
+                               Vector<double>   &value) const;
+
+  };
+
+
+  template <int dim>
+  double
+  RightHandSide<dim>::value (const Point<dim>  &/*p*/,
+                             const unsigned int /*component*/) const
+  {
+    return 0;
+  }
+
+
+  template <int dim>
+  void
+  RightHandSide<dim>::vector_value (const Point<dim> &p,
+                                    Vector<double>   &values) const
+  {
+    for (unsigned int c=0; c<this->n_components; ++c)
+      values(c) = RightHandSide<dim>::value (p, c);
+  }
+
+
+
+
+
+  template <class Matrix, class Preconditioner>
+  class InverseMatrix : public Preconditioner
+  {
+  public:
+    InverseMatrix (const Matrix         &m,
+                   const Preconditioner &preconditioner,
+                   const IndexSet       &locally_owned,
+                   const MPI_Comm       &mpi_communicator);
+
+    void vmult (TrilinosWrappers::MPI::Vector       &dst,
+                const TrilinosWrappers::MPI::Vector &src) const;
+
+  private:
+    const SmartPointer<const Matrix> matrix;
+    const SmartPointer<const Preconditioner> preconditioner;
+
+    const MPI_Comm *mpi_communicator;
+    mutable TrilinosWrappers::MPI::Vector tmp;
+  };
+
+
+  template <class Matrix, class Preconditioner>
+  InverseMatrix<Matrix,Preconditioner>::InverseMatrix
+  (const Matrix         &m,
+   const Preconditioner &preconditioner,
+   const IndexSet       &locally_owned,
+   const MPI_Comm       &mpi_communicator)
+    :
+    matrix (&m),
+    preconditioner (&preconditioner),
+    mpi_communicator (&mpi_communicator),
+    tmp(locally_owned, mpi_communicator)
+  {}
+
+
+
+  template <class Matrix, class Preconditioner>
+  void InverseMatrix<Matrix,Preconditioner>::vmult
+  (TrilinosWrappers::MPI::Vector       &dst,
+   const TrilinosWrappers::MPI::Vector &src) const
+  {
+    SolverControl solver_control (src.size(), 1e-6*src.l2_norm(), false, false);
+    TrilinosWrappers::SolverCG    cg (solver_control,
+                                      TrilinosWrappers::SolverCG::AdditionalData());
+
+    tmp = 0.;
+    cg.solve (*matrix, tmp, src, *preconditioner);
+    dst = tmp;
+  }
+
+
+
+  template <class Preconditioner>
+  class SchurComplement : public TrilinosWrappers::SparseMatrix
+  {
+  public:
+    SchurComplement ( const TrilinosWrappers::BlockSparseMatrix &system_matrix,
+                      const InverseMatrix<TrilinosWrappers::SparseMatrix,
+                      Preconditioner> &A_inverse,
+                      const IndexSet &owned_pres,
+                      const IndexSet &relevant_pres,
+                      const MPI_Comm &mpi_communicator);
+
+    void vmult (TrilinosWrappers::MPI::Vector       &dst,
+                const TrilinosWrappers::MPI::Vector &src) const;
+
+  private:
+    const SmartPointer<const TrilinosWrappers::BlockSparseMatrix> system_matrix;
+    const SmartPointer<const InverseMatrix<TrilinosWrappers::SparseMatrix,
+          Preconditioner> > A_inverse;
+    mutable TrilinosWrappers::MPI::Vector tmp1, tmp2;
+  };
+
+
+
+  template <class Preconditioner>
+  SchurComplement<Preconditioner>::
+  SchurComplement (const TrilinosWrappers::BlockSparseMatrix &system_matrix,
+                   const InverseMatrix<TrilinosWrappers::SparseMatrix,
+                   Preconditioner> &A_inverse,
+                   const IndexSet &owned_vel,
+                   const IndexSet &relevant_vel,
+                   const MPI_Comm &mpi_communicator)
+    :
+    system_matrix (&system_matrix),
+    A_inverse (&A_inverse),
+    tmp1 (owned_vel, mpi_communicator),
+    tmp2 (tmp1)
+  {}
+
+
+  template <class Preconditioner>
+  void SchurComplement<Preconditioner>::vmult
+  (TrilinosWrappers::MPI::Vector       &dst,
+   const TrilinosWrappers::MPI::Vector &src) const
+  {
+    system_matrix->block(0,1).vmult (tmp1, src);
+    A_inverse->vmult (tmp2, tmp1);
+    system_matrix->block(1,0).vmult (dst, tmp2);
+  }
+
+
+
+
+  template <int dim>
+  StokesProblem<dim>::StokesProblem (const unsigned int degree)
+    :
+    degree (degree),
+    mpi_communicator (MPI_COMM_WORLD),
+    triangulation (mpi_communicator/*,
+                   Triangulation<dim>::maximum_smoothing*/),
+    fe (FE_Q<dim>(degree+1), dim,
+        FE_Q<dim>(degree)  , 1),
+    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 <int dim>
+  void StokesProblem<dim>::setup_dofs ()
+  {
+    dof_handler.distribute_dofs (fe);
+
+    std::vector<unsigned int> block_component (dim+1,0);
+    block_component[dim] = 1;
+    DoFRenumbering::component_wise (dof_handler, block_component);
+
+    std::vector<types::global_dof_index> dofs_per_block (2);
+    DoFTools::count_dofs_per_block (dof_handler, dofs_per_block, block_component);
+    const unsigned int n_u = dofs_per_block[0],
+                       n_p = dofs_per_block[1];
+
+    {
+      owned_partitioning.clear();
+      IndexSet locally_owned_dofs = dof_handler.locally_owned_dofs();
+      owned_partitioning.push_back(locally_owned_dofs.get_view(0, n_u));
+      owned_partitioning.push_back(locally_owned_dofs.get_view(n_u, n_u+n_p));
+
+      relevant_partitioning.clear();
+      IndexSet locally_relevant_dofs;
+      DoFTools::extract_locally_relevant_dofs (dof_handler, locally_relevant_dofs);
+      relevant_partitioning.push_back(locally_relevant_dofs.get_view(0, n_u));
+      relevant_partitioning.push_back(locally_relevant_dofs.get_view(n_u, n_u+n_p));
+    }
+
+    {
+      constraints.clear ();
+
+      FEValuesExtractors::Vector velocities(0);
+      FEValuesExtractors::Scalar pressure(dim);
+
+      DoFTools::make_hanging_node_constraints (dof_handler,
+                                               constraints);
+#ifdef PERIODIC
+      VectorTools::interpolate_boundary_values (dof_handler,
+                                                3,
+                                                BoundaryValues<dim>(),
+                                                constraints,
+                                                fe.component_mask(velocities));
+#endif
+      VectorTools::interpolate_boundary_values (dof_handler,
+                                                0,
+                                                BoundaryValues<dim>(),//ZeroFunction<dim>(dim+1),
+                                                constraints,
+                                                fe.component_mask(velocities));
+      VectorTools::interpolate_boundary_values (dof_handler,
+                                                1,
+                                                BoundaryValues<dim>(),//ZeroFunction<dim>(dim+1),
+                                                constraints,
+                                                fe.component_mask(velocities));
+
+#ifdef PERIODIC
+      std::vector<GridTools::PeriodicFacePair<typename DoFHandler<dim>::cell_iterator> >
+      periodicity_vector;
+
+      FullMatrix<double> matrix(dim);
+      matrix[0][1]=1.;
+      matrix[1][0]=-1.;
+
+      std::vector<unsigned int> first_vector_components;
+      first_vector_components.push_back(0);
+
+      GridTools::collect_periodic_faces
+      (dof_handler, 2, 3, 1, periodicity_vector, Tensor<1,dim>(),
+       matrix, first_vector_components);
+
+      DoFTools::make_periodicity_constraints<DoFHandler<dim> >
+      (periodicity_vector, constraints, fe.component_mask(velocities));
+#endif
+    }
+
+    constraints.close ();
+
+    {
+      TrilinosWrappers::BlockSparsityPattern bsp
+      (owned_partitioning, owned_partitioning,
+       relevant_partitioning, mpi_communicator);
+
+      DoFTools::make_sparsity_pattern
+      (dof_handler, bsp, constraints, false,
+       Utilities::MPI::this_mpi_process(mpi_communicator));
+
+      bsp.compress();
+
+      system_matrix.reinit (bsp);
+    }
+
+    system_rhs.reinit (owned_partitioning,
+                       mpi_communicator);
+    solution.reinit (owned_partitioning, relevant_partitioning,
+                     mpi_communicator);
+  }
+
+
+
+  template <int dim>
+  void StokesProblem<dim>::assemble_system ()
+  {
+    system_matrix=0.;
+    system_rhs=0.;
+
+    QGauss<dim>   quadrature_formula(degree+2);
+
+    FEValues<dim> fe_values (fe, quadrature_formula,
+                             update_values    |
+                             update_quadrature_points  |
+                             update_JxW_values |
+                             update_gradients);
+
+    const unsigned int   dofs_per_cell   = fe.dofs_per_cell;
+
+    const unsigned int   n_q_points      = quadrature_formula.size();
+
+    FullMatrix<double>   local_matrix (dofs_per_cell, dofs_per_cell);
+    Vector<double>       local_rhs (dofs_per_cell);
+
+    std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
+
+    const RightHandSide<dim>          right_hand_side;
+    std::vector<Vector<double> >      rhs_values (n_q_points,
+                                                  Vector<double>(dim+1));
+
+    const FEValuesExtractors::Vector velocities (0);
+    const FEValuesExtractors::Scalar pressure (dim);
+
+    std::vector<SymmetricTensor<2,dim> > symgrad_phi_u (dofs_per_cell);
+    std::vector<double>                  div_phi_u   (dofs_per_cell);
+    std::vector<double>                  phi_p       (dofs_per_cell);
+
+    typename DoFHandler<dim>::active_cell_iterator
+    cell = dof_handler.begin_active(),
+    endc = dof_handler.end();
+    for (; cell!=endc; ++cell)
+      if (cell->is_locally_owned())
+        {
+          fe_values.reinit (cell);
+          local_matrix = 0;
+          local_rhs = 0;
+
+          right_hand_side.vector_value_list(fe_values.get_quadrature_points(),
+                                            rhs_values);
+
+          for (unsigned int q=0; q<n_q_points; ++q)
+            {
+              for (unsigned int k=0; k<dofs_per_cell; ++k)
+                {
+                  symgrad_phi_u[k] = fe_values[velocities].symmetric_gradient (k, q);
+                  div_phi_u[k]     = fe_values[velocities].divergence (k, q);
+                  phi_p[k]         = fe_values[pressure].value (k, q);
+                }
+
+              for (unsigned int i=0; i<dofs_per_cell; ++i)
+                {
+                  for (unsigned int j=0; j<=i; ++j)
+                    {
+                      local_matrix(i,j) += (symgrad_phi_u[i] * symgrad_phi_u[j]
+                                            - div_phi_u[i] * phi_p[j]
+                                            - phi_p[i] * div_phi_u[j]
+                                            + phi_p[i] * phi_p[j])
+                                           * fe_values.JxW(q);
+
+                    }
+
+                  const unsigned int component_i =
+                    fe.system_to_component_index(i).first;
+                  local_rhs(i) += fe_values.shape_value(i,q) *
+                                  rhs_values[q](component_i) *
+                                  fe_values.JxW(q);
+                }
+            }
+
+
+          for (unsigned int i=0; i<dofs_per_cell; ++i)
+            for (unsigned int j=i+1; j<dofs_per_cell; ++j)
+              local_matrix(i,j) = local_matrix(j,i);
+
+          cell->get_dof_indices (local_dof_indices);
+          constraints.distribute_local_to_global (local_matrix, local_rhs,
+                                                  local_dof_indices,
+                                                  system_matrix, system_rhs);
+        }
+
+    system_matrix.compress (VectorOperation::add);
+    system_rhs.compress (VectorOperation::add);
+
+    pcout << "   Computing preconditioner..." << std::endl << std::flush;
+  }
+
+
+
+
+  template <int dim>
+  void StokesProblem<dim>::solve ()
+  {
+    TrilinosWrappers::PreconditionJacobi A_preconditioner;
+    A_preconditioner.initialize(system_matrix.block(0,0));
+
+    const InverseMatrix<TrilinosWrappers::SparseMatrix,
+          TrilinosWrappers::PreconditionJacobi>
+          A_inverse (system_matrix.block(0,0),
+                     A_preconditioner,
+                     owned_partitioning[0],
+                     mpi_communicator);
+
+    TrilinosWrappers::MPI::BlockVector tmp (owned_partitioning,
+                                            mpi_communicator);
+
+    {
+      TrilinosWrappers::MPI::Vector schur_rhs (owned_partitioning[1],
+                                               mpi_communicator);
+      A_inverse.vmult (tmp.block(0), system_rhs.block(0));
+      system_matrix.block(1,0).vmult (schur_rhs, tmp.block(0));
+      schur_rhs -= system_rhs.block(1);
+
+      SchurComplement<TrilinosWrappers::PreconditionJacobi>
+      schur_complement (system_matrix, A_inverse,
+                        owned_partitioning[0], relevant_partitioning[0],
+                        mpi_communicator);
+
+      SolverControl solver_control (solution.block(1).size(),
+                                    1e-6*schur_rhs.l2_norm(), false, false);
+      SolverCG<TrilinosWrappers::MPI::Vector> cg(solver_control);
+
+      TrilinosWrappers::PreconditionAMGMueLu preconditioner;
+      preconditioner.initialize (system_matrix.block(1,1));
+
+      InverseMatrix<TrilinosWrappers::SparseMatrix,
+                    TrilinosWrappers::PreconditionAMGMueLu>
+                    m_inverse (system_matrix.block(1,1), preconditioner,
+                               owned_partitioning[1], mpi_communicator);
+
+      cg.solve (schur_complement,
+                tmp.block(1),
+                schur_rhs,
+                preconditioner);
+
+      constraints.distribute (tmp);
+      solution.block(1)=tmp.block(1);
+    }
+
+    {
+      system_matrix.block(0,1).vmult (tmp.block(0), tmp.block(1));
+      tmp.block(0) *= -1;
+      tmp.block(0) += system_rhs.block(0);
+
+      A_inverse.vmult (tmp.block(0), tmp.block(0));
+
+      constraints.distribute (tmp);
+      solution.block(0)=tmp.block(0);
+    }
+  }
+
+  template <int dim>
+  void StokesProblem<dim>::get_point_value
+  (const Point<dim> point, const int proc, Vector<double> &value) const
+  {
+    typename DoFHandler<dim>::active_cell_iterator cell
+      = GridTools::find_active_cell_around_point (dof_handler, point);
+
+    if (cell->is_locally_owned())
+      VectorTools::point_value (dof_handler, solution,
+                                point, value);
+
+    std::vector<double> tmp (value.size());
+    for (unsigned int i=0; i<value.size(); ++i)
+      tmp[i]=value[i];
+
+    std::vector<double> tmp2 (value.size());
+    MPI_Reduce(&(tmp[0]), &(tmp2[0]), value.size(), MPI_DOUBLE,
+               MPI_SUM, proc, mpi_communicator);
+
+    for (unsigned int i=0; i<value.size(); ++i)
+      value[i]=tmp2[i];
+  }
+
+  template <int dim>
+  void StokesProblem<dim>::check_periodicity (const unsigned int cycle) const
+  {}
+
+  template <>
+  void StokesProblem<2>::check_periodicity (const unsigned int cycle) const
+  {
+    unsigned int n_points = 4;
+    for (unsigned int i = 0; i<cycle; i++)
+      n_points*=2;
+
+    for (unsigned int i=1; i< n_points; i++)
+      {
+        Vector<double> value1(3);
+        Vector<double> value2(3);
+
+        Point<2> point1;
+        point1(0)=0;
+        point1(1)=.5*(1.+1.*i/n_points);
+        Point<2> point2;
+        point2(0)=.5*(1.+1.*i/n_points);
+        point2(1)=0.;
+
+        get_point_value (point1, 0, value1);
+        get_point_value (point2, 0, value2);
+
+        if (Utilities::MPI::this_mpi_process(mpi_communicator)==0)
+          {
+            pcout << point1 << "\t" << value1[0] << "\t" << value1[1] << std::endl;
+            if (std::abs(value2[0]-value1[1])>1e-8)
+              {
+                std::cout<<point1<< "\t" << value1[1] << std::endl;
+                std::cout<<point2<< "\t" << value2[0] << std::endl;
+                Assert(false, ExcInternalError());
+              }
+            if (std::abs(value2[1]+value1[0])>1e-8)
+              {
+                std::cout<<point1<< "\t" << value1[0] << std::endl;
+                std::cout<<point2<< "\t" << value2[1] << std::endl;
+                Assert(false, ExcInternalError());
+              }
+          }
+      }
+  }
+
+
+  template <int dim>
+  void
+  StokesProblem<dim>::output_results (const unsigned int refinement_cycle)  const
+  {
+    std::vector<std::string> solution_names (dim, "velocity");
+    solution_names.push_back ("pressure");
+
+    std::vector<DataComponentInterpretation::DataComponentInterpretation>
+    data_component_interpretation
+    (dim, DataComponentInterpretation::component_is_part_of_vector);
+    data_component_interpretation
+    .push_back (DataComponentInterpretation::component_is_scalar);
+
+    DataOut<dim> data_out;
+    data_out.attach_dof_handler (dof_handler);
+    data_out.add_data_vector (solution, solution_names,
+                              DataOut<dim>::type_dof_data,
+                              data_component_interpretation);
+    Vector<float> subdomain (triangulation.n_active_cells());
+    for (unsigned int i=0; i<subdomain.size(); ++i)
+      subdomain(i) = triangulation.locally_owned_subdomain();
+    data_out.add_data_vector (subdomain, "subdomain");
+    data_out.build_patches (degree+1);
+
+    std::ostringstream filename;
+    filename << "solution-"
+             << Utilities::int_to_string (refinement_cycle, 2)
+             << "."
+             << Utilities::int_to_string (triangulation.locally_owned_subdomain(),2)
+             << ".vtu";
+
+    std::ofstream output (filename.str().c_str());
+    data_out.write_vtu (output);
+
+    if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
+      {
+        std::vector<std::string> filenames;
+        for (unsigned int i=0; i<Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD); ++i)
+          filenames.push_back (std::string("solution-") +
+                               Utilities::int_to_string (refinement_cycle, 2) +
+                               "." +
+                               Utilities::int_to_string(i, 2) +
+                               ".vtu");
+        const std::string
+        pvtu_master_filename = ("solution-" +
+                                Utilities::int_to_string (refinement_cycle, 2) +
+                                ".pvtu");
+        std::ofstream pvtu_master (pvtu_master_filename.c_str());
+        data_out.write_pvtu_record (pvtu_master, filenames);
+      }
+  }
+
+
+
+  template <int dim>
+  void
+  StokesProblem<dim>::refine_mesh ()
+  {
+
+    Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
+
+    FEValuesExtractors::Scalar pressure(dim);
+    KellyErrorEstimator<dim>::estimate (dof_handler,
+                                        QGauss<dim-1>(degree+1),
+                                        typename FunctionMap<dim>::type(),
+                                        solution,
+                                        estimated_error_per_cell,
+                                        fe.component_mask(pressure));
+
+    parallel::distributed::GridRefinement::
+    refine_and_coarsen_fixed_number (triangulation,
+                                     estimated_error_per_cell,
+                                     0.3, 0.0);
+    triangulation.execute_coarsening_and_refinement ();
+  }
+
+
+
+  template <int dim>
+  void StokesProblem<dim>::run ()
+  {
+    Point<dim> center;
+    const double inner_radius = .5;
+    const double outer_radius = 1.;
+
+    GridGenerator::quarter_hyper_shell (triangulation,
+                                        center,
+                                        inner_radius,
+                                        outer_radius,
+                                        0,
+                                        true);
+
+#ifdef PERIODIC
+    std::vector<GridTools::PeriodicFacePair<typename parallel::distributed::Triangulation<dim>::cell_iterator> >
+    periodicity_vector;
+
+    FullMatrix<double> matrix(dim);
+    matrix[0][1]=1.;
+    matrix[1][0]=-1.;
+
+    std::vector<unsigned int> first_vector_components;
+    first_vector_components.push_back(0);
+
+    GridTools::collect_periodic_faces
+    (triangulation, 2, 3, 1, periodicity_vector, Tensor<1,dim>(),
+     matrix, first_vector_components);
+
+    triangulation.add_periodicity(periodicity_vector);
+#endif
+
+
+    triangulation.set_boundary(0, boundary);
+    triangulation.set_boundary(1, boundary);
+
+    triangulation.refine_global (4-dim);
+
+    for (unsigned int refinement_cycle = 0; refinement_cycle<3;
+         ++refinement_cycle)
+      {
+        pcout << "Refinement cycle " << refinement_cycle << std::endl;
+
+        if (refinement_cycle > 0)
+          refine_mesh ();
+
+        setup_dofs ();
+
+        pcout << "   Assembling..." << std::endl << std::flush;
+        assemble_system ();
+
+        pcout << "   Solving..." << std::endl << std::flush;
+        solve ();
+
+        check_periodicity(refinement_cycle);
+//         output_results (refinement_cycle);
+
+        pcout << std::endl;
+      }
+  }
+}
+
+
+
+int main (int argc, char *argv[])
+{
+  try
+    {
+      using namespace dealii;
+      using namespace Step22;
+
+      Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+      deallog.depth_console (0);
+
+      if (Utilities::MPI::this_mpi_process (MPI_COMM_WORLD)==0)
+        {
+          std::ofstream logfile("output");
+          deallog.attach(logfile, false);
+          deallog.depth_console(0);
+          deallog.threshold_double(1.e-10);
+          {
+            StokesProblem<2> flow_problem(1);
+            flow_problem.run ();
+          }
+        }
+      else
+        {
+          StokesProblem<2> flow_problem(1);
+          flow_problem.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;
+}
diff --git a/tests/mpi/muelu_periodicity.with_trilinos=true.mpirun=3.output b/tests/mpi/muelu_periodicity.with_trilinos=true.mpirun=3.output
new file mode 100644 (file)
index 0000000..aa60359
--- /dev/null
@@ -0,0 +1,40 @@
+Refinement cycle 0
+   Assembling...
+   Computing preconditioner...
+   Solving...
+0.00000 0.625000       -1.60000        0.00000
+0.00000 0.750000       -1.33333        0.00000
+0.00000 0.875000       -1.14286        0.00000
+
+Refinement cycle 1
+   Assembling...
+   Computing preconditioner...
+   Solving...
+0.00000 0.562500       -1.77778        0.00000
+0.00000 0.625000       -1.60000        0.00000
+0.00000 0.687500       -1.45455        0.00000
+0.00000 0.750000       -1.33333        5.49010e-13
+0.00000 0.812500       -1.23077        2.80537e-13
+0.00000 0.875000       -1.14286        0.00000
+0.00000 0.937500       -1.06667        0.00000
+
+Refinement cycle 2
+   Assembling...
+   Computing preconditioner...
+   Solving...
+0.00000 0.531250       -1.88235        0.00000
+0.00000 0.562500       -1.77778        0.00000
+0.00000 0.593750       -1.68421        0.00000
+0.00000 0.625000       -1.60000        0.00000
+0.00000 0.656250       -1.52381        0.00000
+0.00000 0.687500       -1.45455        0.00000
+0.00000 0.718750       -1.39130        1.95566e-17
+0.00000 0.750000       -1.33333        5.49063e-13
+0.00000 0.781250       -1.28000        0.00000
+0.00000 0.812500       -1.23077        2.80586e-13
+0.00000 0.843750       -1.18519        1.69980e-17
+0.00000 0.875000       -1.14286        0.00000
+0.00000 0.906250       -1.10357        0.00000
+0.00000 0.937500       -1.06667        0.00000
+0.00000 0.968750       -1.03214        0.00000
+
diff --git a/tests/mpi/muelu_periodicity.with_trilinos=true.mpirun=5.output b/tests/mpi/muelu_periodicity.with_trilinos=true.mpirun=5.output
new file mode 100644 (file)
index 0000000..aa60359
--- /dev/null
@@ -0,0 +1,40 @@
+Refinement cycle 0
+   Assembling...
+   Computing preconditioner...
+   Solving...
+0.00000 0.625000       -1.60000        0.00000
+0.00000 0.750000       -1.33333        0.00000
+0.00000 0.875000       -1.14286        0.00000
+
+Refinement cycle 1
+   Assembling...
+   Computing preconditioner...
+   Solving...
+0.00000 0.562500       -1.77778        0.00000
+0.00000 0.625000       -1.60000        0.00000
+0.00000 0.687500       -1.45455        0.00000
+0.00000 0.750000       -1.33333        5.49010e-13
+0.00000 0.812500       -1.23077        2.80537e-13
+0.00000 0.875000       -1.14286        0.00000
+0.00000 0.937500       -1.06667        0.00000
+
+Refinement cycle 2
+   Assembling...
+   Computing preconditioner...
+   Solving...
+0.00000 0.531250       -1.88235        0.00000
+0.00000 0.562500       -1.77778        0.00000
+0.00000 0.593750       -1.68421        0.00000
+0.00000 0.625000       -1.60000        0.00000
+0.00000 0.656250       -1.52381        0.00000
+0.00000 0.687500       -1.45455        0.00000
+0.00000 0.718750       -1.39130        1.95566e-17
+0.00000 0.750000       -1.33333        5.49063e-13
+0.00000 0.781250       -1.28000        0.00000
+0.00000 0.812500       -1.23077        2.80586e-13
+0.00000 0.843750       -1.18519        1.69980e-17
+0.00000 0.875000       -1.14286        0.00000
+0.00000 0.906250       -1.10357        0.00000
+0.00000 0.937500       -1.06667        0.00000
+0.00000 0.968750       -1.03214        0.00000
+
diff --git a/tests/mpi/muelu_periodicity.with_trilinos=true.mpirun=7.output b/tests/mpi/muelu_periodicity.with_trilinos=true.mpirun=7.output
new file mode 100644 (file)
index 0000000..aa60359
--- /dev/null
@@ -0,0 +1,40 @@
+Refinement cycle 0
+   Assembling...
+   Computing preconditioner...
+   Solving...
+0.00000 0.625000       -1.60000        0.00000
+0.00000 0.750000       -1.33333        0.00000
+0.00000 0.875000       -1.14286        0.00000
+
+Refinement cycle 1
+   Assembling...
+   Computing preconditioner...
+   Solving...
+0.00000 0.562500       -1.77778        0.00000
+0.00000 0.625000       -1.60000        0.00000
+0.00000 0.687500       -1.45455        0.00000
+0.00000 0.750000       -1.33333        5.49010e-13
+0.00000 0.812500       -1.23077        2.80537e-13
+0.00000 0.875000       -1.14286        0.00000
+0.00000 0.937500       -1.06667        0.00000
+
+Refinement cycle 2
+   Assembling...
+   Computing preconditioner...
+   Solving...
+0.00000 0.531250       -1.88235        0.00000
+0.00000 0.562500       -1.77778        0.00000
+0.00000 0.593750       -1.68421        0.00000
+0.00000 0.625000       -1.60000        0.00000
+0.00000 0.656250       -1.52381        0.00000
+0.00000 0.687500       -1.45455        0.00000
+0.00000 0.718750       -1.39130        1.95566e-17
+0.00000 0.750000       -1.33333        5.49063e-13
+0.00000 0.781250       -1.28000        0.00000
+0.00000 0.812500       -1.23077        2.80586e-13
+0.00000 0.843750       -1.18519        1.69980e-17
+0.00000 0.875000       -1.14286        0.00000
+0.00000 0.906250       -1.10357        0.00000
+0.00000 0.937500       -1.06667        0.00000
+0.00000 0.968750       -1.03214        0.00000
+
diff --git a/tests/trilinos/precondition_muelu_dgp.cc b/tests/trilinos/precondition_muelu_dgp.cc
new file mode 100644 (file)
index 0000000..1435129
--- /dev/null
@@ -0,0 +1,252 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2015 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 MueLu preconditioner
+
+#include "../tests.h"
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <deal.II/fe/fe_dgp.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 <deal.II/meshworker/dof_info.h>
+#include <deal.II/meshworker/integration_info.h>
+#include <deal.II/meshworker/assembler.h>
+#include <deal.II/meshworker/loop.h>
+#include <deal.II/integrators/laplace.h>
+
+
+#include <fstream>
+#include <iomanip>
+
+
+template <int dim>
+class MatrixIntegrator : public MeshWorker::LocalIntegrator<dim>
+{
+public:
+  void cell(MeshWorker::DoFInfo<dim> &dinfo,
+            typename MeshWorker::IntegrationInfo<dim> &info) const;
+  void boundary(MeshWorker::DoFInfo<dim> &dinfo,
+                typename MeshWorker::IntegrationInfo<dim> &info) const;
+  void face(MeshWorker::DoFInfo<dim> &dinfo1,
+            MeshWorker::DoFInfo<dim> &dinfo2,
+            typename MeshWorker::IntegrationInfo<dim> &info1,
+            typename MeshWorker::IntegrationInfo<dim> &info2) const;
+};
+
+template <int dim>
+void MatrixIntegrator<dim>
+::cell(MeshWorker::DoFInfo<dim> &dinfo,
+       typename MeshWorker::IntegrationInfo<dim> &info) const
+{
+  LocalIntegrators::Laplace::cell_matrix(dinfo.matrix(0,false).matrix,
+                                         info.fe_values());
+}
+
+template <int dim>
+void MatrixIntegrator<dim>
+::boundary(MeshWorker::DoFInfo<dim> &dinfo,
+           typename MeshWorker::IntegrationInfo<dim> &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 <int dim>
+void MatrixIntegrator<dim>
+::face(MeshWorker::DoFInfo<dim> &dinfo1,
+       MeshWorker::DoFInfo<dim> &dinfo2,
+       typename MeshWorker::IntegrationInfo<dim> &info1,
+       typename MeshWorker::IntegrationInfo<dim> &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 <int dim>
+class Step4
+{
+public:
+  Step4 ();
+  void run ();
+
+private:
+  void make_grid ();
+  void setup_system();
+  void solve ();
+
+  Triangulation<dim>   triangulation;
+  FE_DGP<dim>          fe;
+  DoFHandler<dim>      dof_handler;
+
+  TrilinosWrappers::SparseMatrix system_matrix;
+
+  Vector<double>       solution;
+  Vector<double>       system_rhs;
+};
+
+
+
+
+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);
+
+  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<dim> mapping;
+  MeshWorker::IntegrationInfoBox<dim> info_box;
+  UpdateFlags update_flags = update_values | update_gradients;
+  info_box.add_update_flags_all(update_flags);
+  info_box.initialize(fe, mapping);
+
+  MeshWorker::DoFInfo<dim> dof_info(dof_handler);
+  MeshWorker::Assembler::MatrixSimple<TrilinosWrappers::SparseMatrix> assembler;
+  assembler.initialize(system_matrix);
+  MatrixIntegrator<dim> integrator;
+  MeshWorker::integration_loop<dim, dim>(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 <int dim>
+void Step4<dim>::solve ()
+{
+
+  deallog.push(Utilities::int_to_string(dof_handler.n_dofs(),5));
+  TrilinosWrappers::PreconditionAMGMueLu preconditioner;
+  TrilinosWrappers::PreconditionAMGMueLu::AdditionalData data;
+  DoFTools::extract_constant_modes(dof_handler, std::vector<bool>(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 <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();
+      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_muelu_dgp.with_64bit_indices=off.output b/tests/trilinos/precondition_muelu_dgp.with_64bit_indices=off.output
new file mode 100644 (file)
index 0000000..f7ee5d0
--- /dev/null
@@ -0,0 +1,5 @@
+
+DEAL:12288:cg::Starting value 1970.22
+DEAL:12288:cg::Convergence step 19 value 0
+DEAL:49152:cg::Starting value 179304.
+DEAL:49152:cg::Convergence step 31 value 0
diff --git a/tests/trilinos/precondition_muelu_q_iso_q1.cc b/tests/trilinos/precondition_muelu_q_iso_q1.cc
new file mode 100644 (file)
index 0000000..c4e1366
--- /dev/null
@@ -0,0 +1,354 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2015 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 FE_Q elements using FE_Q_iso_Q1 elements
+// for the MueLu preconditioner. The problem is taken from step-4
+
+#include "../tests.h"
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_q_iso_q1.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 assemble_preconditioner ();
+  void solve ();
+
+  Triangulation<dim>   triangulation;
+  FE_Q<dim>            fe;
+  DoFHandler<dim>      dof_handler;
+
+  FE_Q_iso_Q1<dim>     fe_precondition;
+  DoFHandler<dim>      dof_handler_precondition;
+
+  ConstraintMatrix     constraints;
+
+  TrilinosWrappers::SparseMatrix system_matrix;
+  TrilinosWrappers::SparseMatrix preconditioner_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 (3),
+  dof_handler (triangulation),
+  fe_precondition (3),
+  dof_handler_precondition(triangulation)
+{}
+
+
+template <int dim>
+void Step4<dim>::make_grid ()
+{
+  GridGenerator::hyper_cube (triangulation, -1, 1);
+  triangulation.refine_global (4);
+}
+
+
+
+template <int dim>
+void Step4<dim>::setup_system ()
+{
+  dof_handler.distribute_dofs (fe);
+  dof_handler_precondition.distribute_dofs (fe_precondition);
+
+  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);
+  preconditioner_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>::assemble_preconditioner ()
+{
+  QIterated<dim>  quadrature_formula(QGauss<1>(2), fe.degree);
+
+  const RightHandSide<dim> right_hand_side;
+
+  FEValues<dim> fe_values (fe_precondition, quadrature_formula,
+                           update_values   | update_gradients |
+                           update_quadrature_points | update_JxW_values);
+
+  const unsigned int   dofs_per_cell = fe_precondition.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_precondition.begin_active(),
+  endc = dof_handler_precondition.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->get_dof_indices (local_dof_indices);
+      constraints.distribute_local_to_global(cell_matrix,
+                                             local_dof_indices,
+                                             preconditioner_matrix);
+    }
+  preconditioner_matrix.compress(VectorOperation::add);
+}
+
+
+template <int dim>
+void Step4<dim>::solve ()
+{
+
+  // variant 1: solve with MueLu
+  deallog.push(Utilities::int_to_string(dof_handler.n_dofs(),5));
+  deallog.push("MueLu_Q");
+  {
+    solution = 0;
+    SolverControl           solver_control (1000, 1e-12);
+    SolverCG<>              solver (solver_control);
+    TrilinosWrappers::PreconditionAMGMueLu preconditioner;
+    preconditioner.initialize(system_matrix);
+    solver.solve (system_matrix, solution, system_rhs,
+                  preconditioner);
+  }
+  deallog.pop();
+
+  deallog.push("MueLu_Q_iso_Q1");
+  {
+    solution = 0;
+    SolverControl           solver_control (1000, 1e-12);
+    SolverCG<>              solver (solver_control);
+    TrilinosWrappers::PreconditionAMGMueLu preconditioner;
+    preconditioner.initialize(preconditioner_matrix);
+    solver.solve (system_matrix, solution, system_rhs,
+                  preconditioner);
+  }
+  deallog.pop();
+  deallog << std::endl;
+  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();
+      assemble_preconditioner();
+      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_muelu_q_iso_q1.with_64bit_indices=off.output b/tests/trilinos/precondition_muelu_q_iso_q1.with_64bit_indices=off.output
new file mode 100644 (file)
index 0000000..c592763
--- /dev/null
@@ -0,0 +1,11 @@
+
+DEAL:02401:MueLu_Q:cg::Starting value 31.4643
+DEAL:02401:MueLu_Q:cg::Convergence step 45 value 0
+DEAL:02401:MueLu_Q_iso_Q1:cg::Starting value 31.4643
+DEAL:02401:MueLu_Q_iso_Q1:cg::Convergence step 34 value 0
+DEAL:02401::
+DEAL:09409:MueLu_Q:cg::Starting value 44.5271
+DEAL:09409:MueLu_Q:cg::Convergence step 47 value 0
+DEAL:09409:MueLu_Q_iso_Q1:cg::Starting value 44.5271
+DEAL:09409:MueLu_Q_iso_Q1:cg::Convergence step 34 value 0
+DEAL:09409::
diff --git a/tests/trilinos/precondition_muelu_smoother.cc b/tests/trilinos/precondition_muelu_smoother.cc
new file mode 100644 (file)
index 0000000..0d5ab7b
--- /dev/null
@@ -0,0 +1,300 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2015 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 MueLu 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 MueLu
+  deallog.push(Utilities::int_to_string(dof_handler.n_dofs(),5));
+  deallog.push("Chebyshev");
+  TrilinosWrappers::PreconditionAMGMueLu preconditioner;
+  TrilinosWrappers::PreconditionAMGMueLu::AdditionalData data;
+  data.coarse_type = "Amesos-KLU";
+  data.smoother_type = "Chebyshev";
+  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;
+    };
+}
diff --git a/tests/trilinos/precondition_muelu_smoother.with_64bit_indices=off.output b/tests/trilinos/precondition_muelu_smoother.with_64bit_indices=off.output
new file mode 100644 (file)
index 0000000..2195715
--- /dev/null
@@ -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 8 value 0

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.