From: Bruno Turcksin Date: Mon, 6 Oct 2014 14:31:30 +0000 (-0500) Subject: Add MueLu preconditioner using MueLu::MLParameterListInterpreter. X-Git-Tag: v8.3.0-rc1~379^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F214%2Fhead;p=dealii.git Add MueLu preconditioner using MueLu::MLParameterListInterpreter. --- diff --git a/cmake/configure/configure_2_trilinos.cmake b/cmake/configure/configure_2_trilinos.cmake index 05fcc6872c..1e251b6d72 100644 --- a/cmake/configure/configure_2_trilinos.cmake +++ b/cmake/configure/configure_2_trilinos.cmake @@ -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}" diff --git a/doc/news/changes.h b/doc/news/changes.h index 3e672fb1ba..909ba4aa08 100644 --- a/doc/news/changes.h +++ b/doc/news/changes.h @@ -376,6 +376,13 @@ inconvenience this causes.

Specific improvements

    +
  1. 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. +
    + (Bruno Turcksin, 2015/03/11) +
    2. +
  2. 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 diff --git a/include/deal.II/lac/trilinos_precondition.h b/include/deal.II/lac/trilinos_precondition.h index 079c2dd551..eeef76c23b 100644 --- a/include/deal.II/lac/trilinos_precondition.h +++ b/include/deal.II/lac/trilinos_precondition.h @@ -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. // @@ -38,6 +38,11 @@ # include # include +# if DEAL_II_TRILINOS_VERSION_GTE(11,14,0) +# include +# include +# 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 > &constant_modes = std::vector > (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 + * strong coupling is controlled by the variable + * aggregation_threshold, meaning that all elements that are + * not smaller than aggregation_threshold times the diagonal + * element do couple strongly. + */ + double aggregation_threshold; + + /** + * Specifies the constant modes (near null space) of the matrix. This + * parameter tells AMG whether we work on a scalar equation (where the + * near null space only consists of ones) or on a vector-valued + * equation. + */ + std::vector > constant_modes; + + /** + * Determines how many sweeps of the smoother should be performed. When + * the flag elliptic is set to true, i.e., for + * elliptic or almost elliptic problems, the polynomial degree of the + * Chebyshev smoother is set to smoother_sweeps. The term + * sweeps refers to the number of matrix-vector products performed in + * the Chebyshev case. In the non-elliptic case, + * smoother_sweeps sets the number of SSOR relaxation sweeps + * for post-smoothing to be performed. + */ + 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 true, then internal information from + * the ML preconditioner is printed to screen. This can be useful when + * debugging the preconditioner. + */ + bool output_details; + + /** + * Determines which smoother to use for the AMG cycle. Possibilities + * for smoother_type are the following: + *
      + *
    • "Aztec"
      • + *
    • "IFPACK"
      • + *
    • "Jacobi"
      • + *
    • "ML symmetric Gauss-Seidel"
      • + *
    • "symmetric Gauss-Seidel"
      • + *
    • "ML Gauss-Seidel"
      • + *
    • "Gauss-Seidel"
      • + *
    • "block Gauss-Seidel"
      • + *
    • "symmetric block Gauss-Seidel"
      • + *
    • "Chebyshev"
      • + *
    • "MLS"
      • + *
    • "Hiptmair"
      • + *
    • "Amesos-KLU"
      • + *
    • "Amesos-Superlu"
      • + *
    • "Amesos-UMFPACK"
      • + *
    • "Amesos-Superludist"
      • + *
    • "Amesos-MUMPS"
      • + *
    • "user-defined"
      • + *
    • "SuperLU"
      • + *
    • "IFPACK-Chebyshev"
      • + *
    • "self"
      • + *
    • "do-nothing"
      • + *
    • "IC"
      • + *
    • "ICT"
      • + *
    • "ILU"
      • + *
    • "ILUT"
      • + *
    • "Block Chebyshev"
      • + *
    • "IFPACK-Block Chebyshev"
      • + *
    + */ + const char *smoother_type; + + /** + * Determines which solver to use on the coarsest level. The same + * settings as for the smoother type are possible. + */ + const char *coarse_type; + }; + + /** + * 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 ML user's + * guide. 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 + void initialize (const ::dealii::SparseMatrix &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 trilinos_matrix; + }; +#endif + + + /** * A wrapper class for an identity preconditioner for Trilinos matrices. * diff --git a/source/lac/trilinos_precondition.cc b/source/lac/trilinos_precondition.cc index f3e223df81..8012fda8ff 100644 --- a/source/lac/trilinos_precondition.cc +++ b/source/lac/trilinos_precondition.cc @@ -28,6 +28,10 @@ # include # include +#if DEAL_II_TRILINOS_VERSION_GTE(11,14,0) +#include +#include +#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 > &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(additional_data.smoother_sweeps)); + parameter_list.set("cycle applications", + static_cast(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(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 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(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 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 rcp_matrix = Teuchos::rcpFromRef( + *(const_cast(&matrix))); + Teuchos::RCP > muelu_crs_matrix = + Teuchos::rcp(new Xpetra::EpetraCrsMatrix (rcp_matrix)); + Teuchos::RCP > muelu_matrix = + Teuchos::rcp(new Xpetra::CrsMatrixWrap (muelu_crs_matrix)); + + // Create the multigrid hierarchy using ML parameters. + Teuchos::RCP > hierarchy_factory; + hierarchy_factory = Teuchos::rcp( + new MueLu::MLParameterListInterpreter (muelu_parameters)); + Teuchos::RCP > 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 + void + PreconditionAMGMueLu:: + initialize (const ::dealii::SparseMatrix &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(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 &, + const AdditionalData &, const double, + const ::dealii::SparsityPattern *); + template void PreconditionAMGMueLu::initialize (const ::dealii::SparseMatrix &, + 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 index 0000000000..18b0ef6ddb --- /dev/null +++ b/tests/mpi/muelu_periodicity.cc @@ -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 + +#include + +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include + +#include +#include + +#include +#include + +#include +#include +#include + +namespace Step22 +{ + using namespace dealii; + + template + 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 point, const int proc, + Vector &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 boundary; + parallel::distributed::Triangulation triangulation; + FESystem fe; + DoFHandler dof_handler; + + ConstraintMatrix constraints; + std::vector owned_partitioning; + std::vector relevant_partitioning; + + TrilinosWrappers::BlockSparsityPattern sparsity_pattern; + TrilinosWrappers::BlockSparseMatrix system_matrix; + + TrilinosWrappers::MPI::BlockVector solution; + TrilinosWrappers::MPI::BlockVector system_rhs; + + ConditionalOStream pcout; + }; + + + + template + class BoundaryValues : public Function + { + public: + BoundaryValues () : Function(dim+1) {} + + virtual double value (const Point &p, + const unsigned int component = 0) const; + + virtual void vector_value (const Point &p, + Vector &value) const; + }; + + + template + double + BoundaryValues::value (const Point &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 + void + BoundaryValues::vector_value (const Point &p, + Vector &values) const + { + for (unsigned int c=0; cn_components; ++c) + values(c) = BoundaryValues::value (p, c); + } + + + + template + class RightHandSide : public Function + { + public: + RightHandSide () : Function(dim+1) {} + + virtual double value (const Point &p, + const unsigned int component = 0) const; + + virtual void vector_value (const Point &p, + Vector &value) const; + + }; + + + template + double + RightHandSide::value (const Point &/*p*/, + const unsigned int /*component*/) const + { + return 0; + } + + + template + void + RightHandSide::vector_value (const Point &p, + Vector &values) const + { + for (unsigned int c=0; cn_components; ++c) + values(c) = RightHandSide::value (p, c); + } + + + + + + template + 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 matrix; + const SmartPointer preconditioner; + + const MPI_Comm *mpi_communicator; + mutable TrilinosWrappers::MPI::Vector tmp; + }; + + + template + InverseMatrix::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 + void InverseMatrix::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 SchurComplement : public TrilinosWrappers::SparseMatrix + { + public: + SchurComplement ( const TrilinosWrappers::BlockSparseMatrix &system_matrix, + const InverseMatrix &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 system_matrix; + const SmartPointer > A_inverse; + mutable TrilinosWrappers::MPI::Vector tmp1, tmp2; + }; + + + + template + SchurComplement:: + SchurComplement (const TrilinosWrappers::BlockSparseMatrix &system_matrix, + const InverseMatrix &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 + void SchurComplement::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 + StokesProblem::StokesProblem (const unsigned int degree) + : + degree (degree), + mpi_communicator (MPI_COMM_WORLD), + triangulation (mpi_communicator/*, + Triangulation::maximum_smoothing*/), + fe (FE_Q(degree+1), dim, + FE_Q(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 + void StokesProblem::setup_dofs () + { + dof_handler.distribute_dofs (fe); + + std::vector block_component (dim+1,0); + block_component[dim] = 1; + DoFRenumbering::component_wise (dof_handler, block_component); + + std::vector 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(), + constraints, + fe.component_mask(velocities)); +#endif + VectorTools::interpolate_boundary_values (dof_handler, + 0, + BoundaryValues(),//ZeroFunction(dim+1), + constraints, + fe.component_mask(velocities)); + VectorTools::interpolate_boundary_values (dof_handler, + 1, + BoundaryValues(),//ZeroFunction(dim+1), + constraints, + fe.component_mask(velocities)); + +#ifdef PERIODIC + std::vector::cell_iterator> > + periodicity_vector; + + FullMatrix matrix(dim); + matrix[0][1]=1.; + matrix[1][0]=-1.; + + std::vector 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 > + (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 + void StokesProblem::assemble_system () + { + system_matrix=0.; + system_rhs=0.; + + QGauss quadrature_formula(degree+2); + + FEValues 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 local_matrix (dofs_per_cell, dofs_per_cell); + Vector local_rhs (dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); + + const RightHandSide right_hand_side; + std::vector > rhs_values (n_q_points, + Vector(dim+1)); + + const FEValuesExtractors::Vector velocities (0); + const FEValuesExtractors::Scalar pressure (dim); + + std::vector > symgrad_phi_u (dofs_per_cell); + std::vector div_phi_u (dofs_per_cell); + std::vector phi_p (dofs_per_cell); + + typename DoFHandler::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; qget_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 + void StokesProblem::solve () + { + TrilinosWrappers::PreconditionJacobi A_preconditioner; + A_preconditioner.initialize(system_matrix.block(0,0)); + + const InverseMatrix + 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 + 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 cg(solver_control); + + TrilinosWrappers::PreconditionAMGMueLu preconditioner; + preconditioner.initialize (system_matrix.block(1,1)); + + InverseMatrix + 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 + void StokesProblem::get_point_value + (const Point point, const int proc, Vector &value) const + { + typename DoFHandler::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 tmp (value.size()); + for (unsigned int i=0; i tmp2 (value.size()); + MPI_Reduce(&(tmp[0]), &(tmp2[0]), value.size(), MPI_DOUBLE, + MPI_SUM, proc, mpi_communicator); + + for (unsigned int i=0; i + void StokesProblem::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 value1(3); + Vector 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<1e-8) + { + std::cout< + void + StokesProblem::output_results (const unsigned int refinement_cycle) const + { + std::vector solution_names (dim, "velocity"); + solution_names.push_back ("pressure"); + + std::vector + data_component_interpretation + (dim, DataComponentInterpretation::component_is_part_of_vector); + data_component_interpretation + .push_back (DataComponentInterpretation::component_is_scalar); + + DataOut data_out; + data_out.attach_dof_handler (dof_handler); + data_out.add_data_vector (solution, solution_names, + DataOut::type_dof_data, + data_component_interpretation); + Vector subdomain (triangulation.n_active_cells()); + for (unsigned int i=0; i filenames; + for (unsigned int i=0; i + void + StokesProblem::refine_mesh () + { + + Vector estimated_error_per_cell (triangulation.n_active_cells()); + + FEValuesExtractors::Scalar pressure(dim); + KellyErrorEstimator::estimate (dof_handler, + QGauss(degree+1), + typename FunctionMap::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 + void StokesProblem::run () + { + Point 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::cell_iterator> > + periodicity_vector; + + FullMatrix matrix(dim); + matrix[0][1]=1.; + matrix[1][0]=-1.; + + std::vector 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 index 0000000000..aa60359f7a --- /dev/null +++ b/tests/mpi/muelu_periodicity.with_trilinos=true.mpirun=3.output @@ -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 index 0000000000..aa60359f7a --- /dev/null +++ b/tests/mpi/muelu_periodicity.with_trilinos=true.mpirun=5.output @@ -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 index 0000000000..aa60359f7a --- /dev/null +++ b/tests/mpi/muelu_periodicity.with_trilinos=true.mpirun=7.output @@ -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 index 0000000000..1435129757 --- /dev/null +++ b/tests/trilinos/precondition_muelu_dgp.cc @@ -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 +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + + +#include +#include + + +template +class MatrixIntegrator : public MeshWorker::LocalIntegrator +{ +public: + void cell(MeshWorker::DoFInfo &dinfo, + typename MeshWorker::IntegrationInfo &info) const; + void boundary(MeshWorker::DoFInfo &dinfo, + typename MeshWorker::IntegrationInfo &info) const; + void face(MeshWorker::DoFInfo &dinfo1, + MeshWorker::DoFInfo &dinfo2, + typename MeshWorker::IntegrationInfo &info1, + typename MeshWorker::IntegrationInfo &info2) const; +}; + +template +void MatrixIntegrator +::cell(MeshWorker::DoFInfo &dinfo, + typename MeshWorker::IntegrationInfo &info) const +{ + LocalIntegrators::Laplace::cell_matrix(dinfo.matrix(0,false).matrix, + info.fe_values()); +} + +template +void MatrixIntegrator +::boundary(MeshWorker::DoFInfo &dinfo, + typename MeshWorker::IntegrationInfo &info) const +{ + const unsigned int deg = info.fe_values(0).get_fe().degree; + LocalIntegrators::Laplace + ::nitsche_matrix(dinfo.matrix(0,false).matrix, info.fe_values(0), + LocalIntegrators::Laplace:: + compute_penalty(dinfo, dinfo, deg, deg)); +} + +template +void MatrixIntegrator +::face(MeshWorker::DoFInfo &dinfo1, + MeshWorker::DoFInfo &dinfo2, + typename MeshWorker::IntegrationInfo &info1, + typename MeshWorker::IntegrationInfo &info2) const +{ + const unsigned int deg = info1.fe_values(0).get_fe().degree; + LocalIntegrators::Laplace + ::ip_matrix(dinfo1.matrix(0,false).matrix, dinfo1.matrix(0,true).matrix, + dinfo2.matrix(0,true).matrix, dinfo2.matrix(0,false).matrix, + info1.fe_values(0), info2.fe_values(0), + LocalIntegrators::Laplace::compute_penalty(dinfo1, dinfo2, deg, deg)); +} + + +template +class Step4 +{ +public: + Step4 (); + void run (); + +private: + void make_grid (); + void setup_system(); + void solve (); + + Triangulation triangulation; + FE_DGP fe; + DoFHandler dof_handler; + + TrilinosWrappers::SparseMatrix system_matrix; + + Vector solution; + Vector system_rhs; +}; + + + + +template +Step4::Step4 () + : + fe (1), + dof_handler (triangulation) +{} + + +template +void Step4::make_grid () +{ + GridGenerator::hyper_cube (triangulation, -1, 1); + triangulation.refine_global (6); +} + + + +template +void Step4::setup_system () +{ + dof_handler.distribute_dofs (fe); + + CompressedSparsityPattern c_sparsity(dof_handler.n_dofs()); + DoFTools::make_flux_sparsity_pattern (dof_handler, c_sparsity); + system_matrix.reinit (c_sparsity); + + solution.reinit (dof_handler.n_dofs()); + system_rhs.reinit (dof_handler.n_dofs()); + + MappingQ1 mapping; + MeshWorker::IntegrationInfoBox info_box; + UpdateFlags update_flags = update_values | update_gradients; + info_box.add_update_flags_all(update_flags); + info_box.initialize(fe, mapping); + + MeshWorker::DoFInfo dof_info(dof_handler); + MeshWorker::Assembler::MatrixSimple assembler; + assembler.initialize(system_matrix); + MatrixIntegrator integrator; + MeshWorker::integration_loop(dof_handler.begin_active(), + dof_handler.end(), + dof_info, info_box, + integrator, assembler); + + system_matrix.compress(VectorOperation::add); + + for (unsigned int i=0; i +void Step4::solve () +{ + + 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(1,true), + data.constant_modes); + data.smoother_sweeps = 2; + { + solution = 0; + SolverControl solver_control (1000, 1e-10); + SolverCG<> solver (solver_control); + preconditioner.initialize(system_matrix, data); + solver.solve (system_matrix, solution, system_rhs, + preconditioner); + } + deallog.pop(); +} + + + +template +void Step4::run() +{ + for (unsigned int cycle = 0; cycle < 2; ++cycle) + { + if (cycle == 0) + make_grid(); + else + triangulation.refine_global(1); + + setup_system(); + solve(); + } +} + + +int main (int argc, char **argv) +{ + std::ofstream logfile("output"); + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv); + + try + { + Step4<2> test; + test.run(); + } + catch (std::exception &exc) + { + deallog << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + deallog << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/trilinos/precondition_muelu_dgp.with_64bit_indices=off.output b/tests/trilinos/precondition_muelu_dgp.with_64bit_indices=off.output new file mode 100644 index 0000000000..f7ee5d0889 --- /dev/null +++ b/tests/trilinos/precondition_muelu_dgp.with_64bit_indices=off.output @@ -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 index 0000000000..c4e1366043 --- /dev/null +++ b/tests/trilinos/precondition_muelu_q_iso_q1.cc @@ -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 +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + + +template +class Step4 +{ +public: + Step4 (); + void run (); + +private: + void make_grid (); + void setup_system(); + void assemble_system (); + void assemble_preconditioner (); + void solve (); + + Triangulation triangulation; + FE_Q fe; + DoFHandler dof_handler; + + FE_Q_iso_Q1 fe_precondition; + DoFHandler dof_handler_precondition; + + ConstraintMatrix constraints; + + TrilinosWrappers::SparseMatrix system_matrix; + TrilinosWrappers::SparseMatrix preconditioner_matrix; + + Vector solution; + Vector system_rhs; +}; + + +template +class RightHandSide : public Function +{ +public: + RightHandSide () : Function() {} + + virtual double value (const Point &p, + const unsigned int component = 0) const; +}; + + + +template +class BoundaryValues : public Function +{ +public: + BoundaryValues () : Function() {} + + virtual double value (const Point &p, + const unsigned int component = 0) const; +}; + + + + +template +double RightHandSide::value (const Point &p, + const unsigned int /*component*/) const +{ + double return_value = 0; + for (unsigned int i=0; i +double BoundaryValues::value (const Point &p, + const unsigned int /*component*/) const +{ + return p.square(); +} + + + +template +Step4::Step4 () + : + fe (3), + dof_handler (triangulation), + fe_precondition (3), + dof_handler_precondition(triangulation) +{} + + +template +void Step4::make_grid () +{ + GridGenerator::hyper_cube (triangulation, -1, 1); + triangulation.refine_global (4); +} + + + +template +void Step4::setup_system () +{ + dof_handler.distribute_dofs (fe); + dof_handler_precondition.distribute_dofs (fe_precondition); + + constraints.clear(); + std::map boundary_values; + VectorTools::interpolate_boundary_values (dof_handler, + 0, + BoundaryValues(), + constraints); + constraints.close(); + + CompressedSparsityPattern c_sparsity(dof_handler.n_dofs()); + DoFTools::make_sparsity_pattern (dof_handler, c_sparsity, constraints, false); + system_matrix.reinit (c_sparsity); + preconditioner_matrix.reinit(c_sparsity); + + solution.reinit (dof_handler.n_dofs()); + system_rhs.reinit (dof_handler.n_dofs()); +} + + +template +void Step4::assemble_system () +{ + QGauss quadrature_formula(fe.degree+1); + + const RightHandSide right_hand_side; + + FEValues fe_values (fe, quadrature_formula, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); + + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); + + FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); + Vector cell_rhs (dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); + + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + + for (; cell!=endc; ++cell) + { + fe_values.reinit (cell); + cell_matrix = 0; + cell_rhs = 0; + + for (unsigned int q_point=0; q_pointget_dof_indices (local_dof_indices); + constraints.distribute_local_to_global(cell_matrix, cell_rhs, + local_dof_indices, + system_matrix, system_rhs); + } + system_matrix.compress(VectorOperation::add); +} + + + + +template +void Step4::assemble_preconditioner () +{ + QIterated quadrature_formula(QGauss<1>(2), fe.degree); + + const RightHandSide right_hand_side; + + FEValues 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 cell_matrix (dofs_per_cell, dofs_per_cell); + Vector cell_rhs (dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); + + typename DoFHandler::active_cell_iterator + cell = dof_handler_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_pointget_dof_indices (local_dof_indices); + constraints.distribute_local_to_global(cell_matrix, + local_dof_indices, + preconditioner_matrix); + } + preconditioner_matrix.compress(VectorOperation::add); +} + + +template +void Step4::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 +void Step4::run() +{ + for (unsigned int cycle = 0; cycle < 2; ++cycle) + { + if (cycle == 0) + make_grid(); + else + triangulation.refine_global(1); + + setup_system(); + assemble_system(); + 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 index 0000000000..c592763f04 --- /dev/null +++ b/tests/trilinos/precondition_muelu_q_iso_q1.with_64bit_indices=off.output @@ -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 index 0000000000..0d5ab7bed6 --- /dev/null +++ b/tests/trilinos/precondition_muelu_smoother.cc @@ -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 +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + + +template +class Step4 +{ +public: + Step4 (); + void run (); + +private: + void make_grid (); + void setup_system(); + void assemble_system (); + void solve (); + + Triangulation triangulation; + FE_Q fe; + DoFHandler dof_handler; + + ConstraintMatrix constraints; + + TrilinosWrappers::SparseMatrix system_matrix; + + Vector solution; + Vector system_rhs; +}; + + +template +class RightHandSide : public Function +{ +public: + RightHandSide () : Function() {} + + virtual double value (const Point &p, + const unsigned int component = 0) const; +}; + + + +template +class BoundaryValues : public Function +{ +public: + BoundaryValues () : Function() {} + + virtual double value (const Point &p, + const unsigned int component = 0) const; +}; + + + + +template +double RightHandSide::value (const Point &p, + const unsigned int /*component*/) const +{ + double return_value = 0; + for (unsigned int i=0; i +double BoundaryValues::value (const Point &p, + const unsigned int /*component*/) const +{ + return p.square(); +} + + + +template +Step4::Step4 () + : + fe (1), + dof_handler (triangulation) +{} + + +template +void Step4::make_grid () +{ + GridGenerator::hyper_cube (triangulation, -1, 1); + triangulation.refine_global (6); +} + + + +template +void Step4::setup_system () +{ + dof_handler.distribute_dofs (fe); + + constraints.clear(); + std::map boundary_values; + VectorTools::interpolate_boundary_values (dof_handler, + 0, + BoundaryValues(), + constraints); + constraints.close(); + + CompressedSparsityPattern c_sparsity(dof_handler.n_dofs()); + DoFTools::make_sparsity_pattern (dof_handler, c_sparsity, constraints, false); + system_matrix.reinit (c_sparsity); + + solution.reinit (dof_handler.n_dofs()); + system_rhs.reinit (dof_handler.n_dofs()); +} + + +template +void Step4::assemble_system () +{ + QGauss quadrature_formula(fe.degree+1); + + const RightHandSide right_hand_side; + + FEValues fe_values (fe, quadrature_formula, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); + + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); + + FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); + Vector cell_rhs (dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); + + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + + for (; cell!=endc; ++cell) + { + fe_values.reinit (cell); + cell_matrix = 0; + cell_rhs = 0; + + for (unsigned int q_point=0; q_pointget_dof_indices (local_dof_indices); + constraints.distribute_local_to_global(cell_matrix, cell_rhs, + local_dof_indices, + system_matrix, system_rhs); + } + system_matrix.compress(VectorOperation::add); +} + + + +template +void Step4::solve () +{ + + // variant 1: solve with 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 +void Step4::run() +{ + for (unsigned int cycle = 0; cycle < 2; ++cycle) + { + if (cycle == 0) + make_grid(); + else + triangulation.refine_global(1); + + setup_system(); + assemble_system(); + solve(); + } +} + + +int main (int argc, char **argv) +{ + std::ofstream logfile("output"); + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv); + + try + { + Step4<2> test; + test.run(); + } + catch (std::exception &exc) + { + deallog << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + deallog << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/trilinos/precondition_muelu_smoother.with_64bit_indices=off.output b/tests/trilinos/precondition_muelu_smoother.with_64bit_indices=off.output new file mode 100644 index 0000000000..2195715169 --- /dev/null +++ b/tests/trilinos/precondition_muelu_smoother.with_64bit_indices=off.output @@ -0,0 +1,9 @@ + +DEAL:04225:Chebyshev:cg::Starting value 21.8299 +DEAL:04225:Chebyshev:cg::Convergence step 10 value 0 +DEAL:04225:SGS:cg::Starting value 21.8299 +DEAL:04225:SGS:cg::Convergence step 10 value 0 +DEAL:16641:Chebyshev:cg::Starting value 30.8770 +DEAL:16641:Chebyshev:cg::Convergence step 8 value 0 +DEAL:16641:SGS:cg::Starting value 30.8770 +DEAL:16641:SGS:cg::Convergence step 8 value 0