From: Jan Philipp Thiele Date: Mon, 12 Feb 2024 19:53:17 +0000 (+0100) Subject: Add initial support for Amesos2 X-Git-Tag: relicensing~28^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F16602%2Fhead;p=dealii.git Add initial support for Amesos2 --- diff --git a/cmake/configure/configure_20_trilinos.cmake b/cmake/configure/configure_20_trilinos.cmake index f8ba7e57d7..d9f3c280f7 100644 --- a/cmake/configure/configure_20_trilinos.cmake +++ b/cmake/configure/configure_20_trilinos.cmake @@ -23,7 +23,7 @@ set(FEATURE_TRILINOS_DEPENDS MPI) # A list of optional Trilinos modules we use: # set(_deal_ii_trilinos_optional_modules - Belos EpetraExt Kokkos MueLu NOX ROL Sacado SEACAS Tpetra Zoltan + Amesos2 Belos EpetraExt Kokkos MueLu NOX ROL Sacado SEACAS Tpetra Zoltan ) # diff --git a/doc/doxygen/options.dox.in b/doc/doxygen/options.dox.in index 0e8058e970..82d336b9ee 100644 --- a/doc/doxygen/options.dox.in +++ b/doc/doxygen/options.dox.in @@ -241,6 +241,7 @@ PREDEFINED = DOXYGEN=1 \ DEAL_II_SYMENGINE_WITH_LLVM=1 \ DEAL_II_WITH_TBB=1 \ DEAL_II_WITH_TRILINOS=1 \ + DEAL_II_TRILINOS_WITH_AMESOS2=1 \ DEAL_II_TRILINOS_WITH_EPETRAEXT=1 \ DEAL_II_TRILINOS_WITH_MUELU=1 \ DEAL_II_TRILINOS_WITH_NOX=1 \ diff --git a/include/deal.II/base/config.h.in b/include/deal.II/base/config.h.in index c9943b3efe..25367eab2a 100644 --- a/include/deal.II/base/config.h.in +++ b/include/deal.II/base/config.h.in @@ -298,6 +298,7 @@ /* cmake/configure/configure_2_trilinos.cmake */ #cmakedefine DEAL_II_TRILINOS_CXX_SUPPORTS_SACADO_COMPLEX_RAD +#cmakedefine DEAL_II_TRILINOS_WITH_AMESOS2 #cmakedefine DEAL_II_TRILINOS_WITH_BELOS #cmakedefine DEAL_II_TRILINOS_WITH_EPETRAEXT #cmakedefine DEAL_II_TRILINOS_WITH_MUELU diff --git a/include/deal.II/lac/trilinos_tpetra_solver_direct.h b/include/deal.II/lac/trilinos_tpetra_solver_direct.h new file mode 100644 index 0000000000..b422f43956 --- /dev/null +++ b/include/deal.II/lac/trilinos_tpetra_solver_direct.h @@ -0,0 +1,368 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2008 - 2024 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +#ifndef dealii_trilinos_tpetra_solver_direct_h +#define dealii_trilinos_tpetra_solver_direct_h + +#include + +#include "deal.II/grid/grid_generator.h" + +#include + +#ifdef DEAL_II_TRILINOS_WITH_TPETRA +# ifdef DEAL_II_TRILINOS_WITH_AMESOS2 + +# include "deal.II/base/types.h" + +# include +# include +# include + +# include +# include +# include +# include + + +DEAL_II_NAMESPACE_OPEN + +namespace LinearAlgebra +{ + namespace TpetraWrappers + { + // forward declarations +# ifndef DOXYGEN + template + class SparseMatrix; +# endif + + /// The chosen Solver is not supported or configured with Amesos2. + DeclException1( + ExcTrilinosAmesos2SolverUnsupported, + std::string, + << "You tried to select the solver type <" << arg1 << ">\n" + << "but this solver is not supported by Trilinos/Amesos2\n" + << "due to one of the following reasons:\n" + << "* This solver does not exist\n" + << "* This solver is not (yet) supported by Trilinos/Amesos2\n" + << "* Trilinos/Amesos2 was not configured for its use."); + + + /** + * The base class for all direct solvers based on the Amesos2 package + * of Trilinos. + * + * @ingroup TpetraWrappers + */ + template + class SolverDirectBase + { + public: + /** + * Declare the type for container size. + */ + using size_type = dealii::types::signed_global_dof_index; + + /** + * Typedef for the NodeType + */ +# if DEAL_II_TRILINOS_VERSION_GTE(14, 2, 0) + using NodeType = Tpetra::KokkosCompat::KokkosDeviceWrapperNode< + typename MemorySpace::kokkos_space::execution_space, + typename MemorySpace::kokkos_space>; +# else + using NodeType = Kokkos::Compat::KokkosDeviceWrapperNode< + typename MemorySpace::kokkos_space::execution_space, + typename MemorySpace::kokkos_space>; +# endif + + /** + * Typedef for the LinearOperator type + */ + using LinearOperator = Tpetra::Operator; + + /** + * Typedef for the MultiVector type + */ + using MultiVector = Tpetra::MultiVector; + + /** + * Destructor. + */ + virtual ~SolverDirectBase() = default; + + /** + * Initializes the direct solver for the matrix A and creates a + * factorization for it with the package chosen from the additional + * data structure. Note that there is no need for a preconditioner + * here and solve() is not called. + */ + void + initialize(const SparseMatrix &A); + + /** + * Solve the linear system Ax=b based on the + * package set in initialize(). Note the matrix is not refactorized during + * this call. + */ + void + solve(Vector &x, + const Vector &b); + + /** + * Solve the linear system Ax=b. Creates a factorization of the + * matrix with the package chosen from the additional data structure and + * performs the solve. Note that there is no need for a preconditioner + * here. + */ + void + solve(const SparseMatrix &A, + Vector &x, + const Vector &b); + + + /** + * Access to object that controls convergence. + */ + SolverControl & + control() const; + + /** + * Exception + */ + DeclException1(ExcTrilinosError, + int, + << "An error with error number " << arg1 + << " occurred while calling a Trilinos function"); + + protected: + /** + * Constructor. Takes the solver control object and name + * and creates the solver. + */ + SolverDirectBase(SolverControl &cn, + const std::string &solver_type, + const bool output_solver_details = false); + /** + * Actually performs the operations for solving the linear system, + * including the factorization and forward and backward substitution. + */ + void + do_solve(); + + /** + * Reference to the object that controls convergence of the iterative + * solver. In fact, for these Trilinos wrappers, Trilinos does so itself, + * but we copy the data from this object before starting the solution + * process, and copy the data back into it afterwards. + */ + SolverControl &solver_control; + + /** + * A structure that contains the Trilinos solver object. + */ + Teuchos::RCP< + Amesos2::Solver::MatrixType, + MultiVector>> + solver; + + /* + * The set solver type to be handed to the solver factory of Amesos2. + */ + std::string solver_type; + + /** + * Enables/disables the output of solver details (residual in each + * iterations etc.). + */ + bool output_solver_details; + + + /** + * An optional Teuchos::ParameterList for fine tuning the solver. + * Please refer to the Amesos2 manual to see which parameters + * to set for each individual solver. + */ + Teuchos::ParameterList parameter_list; + }; // Base + + + + /** + * A general purpose class to support any solver that the Amesos2 package + * of Trilinos provides. + * + * Notes for users switching from TrilinosWrappers to TpetraWrappers: + * The general interface of this class is kept identical to the Amesos(1) + * wrapper with the notable exception of being templated. + * + * A further addition is the option to fine tune each solver with a + * Teuchos::ParameterList, to change default parameters. + * + * @ingroup TpetraWrappers + */ + template + class SolverDirect : public SolverDirectBase + { + public: + /** + * A structure whose member variables describe details of the algorithm + * to be employed by the surrounding class. An object of this type can be + * used to initialize an object of the surrounding class. + */ + struct AdditionalData + { + AdditionalData(const std::string &solver_name, + const bool output_solver_details = false); + + + /** + * Set the solver type (for third party solver support of Trilinos + * Amesos package). Current possibilities are: + *
    + *
  • "Basker"
    • + *
  • "Cholmod"
    • + *
  • "cuSOLVER"
    • + *
  • "KLU2"
    • + *
  • "LAPACK"
    • + *
  • "MUMPS"
    • + *
  • "PardisoMKL"
    • + *
  • "ShyLUBasker"
    • + *
  • "SuperLU"
    • + *
  • "SuperLU_DIST"
    • + *
  • "SuperLU_MT"
    • + *
  • "Tacho"
    • + *
  • "UMFPACK"
    • + *
+ * Note that the availability of these solvers in deal.II depends on + * which solvers were set when configuring Trilinos. + * Additionally, Amesos2 may add support for solvers not listed here + * that can nevertheless be used if configured. + */ + std::string solver_name; + + /** + * Enables/disables the output of solver details (residual in each + * iterations etc.). + */ + bool output_solver_details; + }; + + /** + * Constructor. Takes the solver control object and creates the solver. + */ + SolverDirect(SolverControl &cn, + const AdditionalData &additional_data = AdditionalData()); + + /** + * Set a parameter list to fine tune the solver. + * The valid parameters depend on the solver you use + * and the possible parameters can be found in the + * documentation of Amesos2. + */ + void + set_pararameter_list(Teuchos::ParameterList ¶meter_list); + }; // SolverDirect + + + /** + * A wrapper class for the solver KLU2 that works in serial and parallel. + * This solver is part of Amesos2 and enabled by default. + * + * The AdditionalData structure allows to pass options specific to this + * solver and the default will result in the same solver as constructing + * SolverDirect with solver_name KLU2. + * + * @ingroup TpetraWrappers + */ + template + class SolverDirectKLU2 : public SolverDirectBase + { + public: + /** + * A structure whose member variables describe details of the algorithm + * to be employed by the surrounding class. An object of this type can be + * used to initialize an object of the surrounding class. + */ + struct AdditionalData + { + AdditionalData(const std::string &transpose_mode = "NOTRANS", + const bool symmetric_mode = false, + const bool equilibrate_matrix = true, + const std::string &column_permutation = "COLAMD", + const std::string &iterative_refinement = "NO", + const bool output_solver_details = false); + /** + * Decide which system to solve + * "NOTRANS": Ax=b (default) + * "TRANS": A^Tx=b + * "CONJ": A^*x=b + */ + std::string transpose_mode; + + /** + * Is A symmetric or not? + */ + bool symmetric_mode; + + /** + * Equilibrate matrix before solving? (default: true) + */ + bool equilibrate_matrix; + + /** + * Choose an ordering strategy + * "COLAMD": approximate minimum degree column ordering (default). + * "NATURAL": natural ordering. + * "MMD_ATA": minimum degree ordering on the structure of A^TA. + * "MMD_AT_PLUS_A": minimum degree ordering on the structure of A^T+A. + */ + std::string column_permutation; + + /** + * Perform iterative refinement? + * "NO": no (default) + * "SINGLE": yes, use single precision for the residual + * "DOUBLE": yes, use double precision for the residual + * "EXTRA": ?? + */ + std::string iterative_refinement; + + + /** + * Enables/disables the output of solver details (residual in each + * iterations etc.). + */ + bool output_solver_details; + }; + + /** + * Constructor. Takes the solver control object and creates the solver. + */ + SolverDirectKLU2( + SolverControl &cn, + const AdditionalData &additional_data = AdditionalData()); + }; // KLU2 + + } // namespace TpetraWrappers +} // namespace LinearAlgebra + +DEAL_II_NAMESPACE_CLOSE + +# endif // DEAL_II_TRILINOS_WITH_AMESOS2 +#endif // DEAL_II_TRILINOS_WITH_TPETRA + +#endif diff --git a/include/deal.II/lac/trilinos_tpetra_solver_direct.templates.h b/include/deal.II/lac/trilinos_tpetra_solver_direct.templates.h new file mode 100644 index 0000000000..238f7ef359 --- /dev/null +++ b/include/deal.II/lac/trilinos_tpetra_solver_direct.templates.h @@ -0,0 +1,245 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2008 - 2024 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +#ifndef dealii_trilinos_tpetra_solver_templates_h +#define dealii_trilinos_tpetra_solver_templates_h + +#include + +#include "deal.II/lac/solver_control.h" + +#include + + +#ifdef DEAL_II_TRILINOS_WITH_TPETRA +# ifdef DEAL_II_TRILINOS_WITH_AMESOS2 + +# include "deal.II/base/template_constraints.h" +# include + +# include + + +DEAL_II_NAMESPACE_OPEN + +namespace LinearAlgebra +{ + namespace TpetraWrappers + { + /* ---------------------- SolverDirectBase ------------------------ */ + + template + SolverDirectBase::SolverDirectBase( + SolverControl &cn, + const std::string &solver_type, + const bool output_solver_details) + : solver_control(cn) + , solver_type(solver_type) + , output_solver_details(output_solver_details) + { + AssertThrow(Amesos2::query(solver_type), + ExcTrilinosAmesos2SolverUnsupported(solver_type)); + } + + + + template + SolverControl & + SolverDirectBase::control() const + { + return solver_control; + } + + + + template + void + SolverDirectBase::initialize( + const SparseMatrix &A) + { + // First set whether we want to print the solver information to screen or + // not. + ConditionalOStream verbose_cout(std::cout, output_solver_details); + + // Next allocate the Amesos2 solver with the concrete solver, if possible. + solver = + Amesos2::create::MatrixType, + MultiVector>(solver_type, A.trilinos_rcp()); + + solver->setParameters(Teuchos::rcpFromRef(parameter_list)); + // Now do the actual factorization, which is a two step procedure. + // The symbolic factorization determines the structure of the inverse, + // while the numeric factorization does to actual computation of L and U + verbose_cout << "Starting symbolic factorization" << std::endl; + solver->symbolicFactorization(); + + verbose_cout << "Starting numeric factorization" << std::endl; + solver->numericFactorization(); + } + + + + template + void + SolverDirectBase::solve( + Vector &x, + const Vector &b) + { + // Assign the empty solution vector + solver->setX(x.trilinos_rcp()); + + // Assign the RHS vector + solver->setB(b.trilinos_rcp()); + // First set whether we want to print the solver information to screen + // or not. + ConditionalOStream verbose_cout(std::cout, output_solver_details); + + verbose_cout << "Starting solve" << std::endl; + solver->solve(); + + // Finally, force the SolverControl object to report convergence + solver_control.check(0, 0); + if (solver_control.last_check() != SolverControl::success) + AssertThrow(false, + SolverControl::NoConvergence(solver_control.last_step(), + solver_control.last_value())); + } + + + + template + void + SolverDirectBase::do_solve() + { + // First set whether we want to print the solver information to screen or + // not. + ConditionalOStream verbose_cout(std::cout, output_solver_details); + + solver->setParameters(Teuchos::rcpFromRef(parameter_list)); + // Now do the actual factorization, which is a two step procedure. + // The symbolic factorization determines the structure of the inverse, + // while the numeric factorization does to actual computation of L and U + verbose_cout << "Starting symbolic factorization" << std::endl; + solver->symbolicFactorization(); + + verbose_cout << "Starting numeric factorization" << std::endl; + solver->numericFactorization(); + + verbose_cout << "Starting solve" << std::endl; + solver->solve(); + + // Finally, force the SolverControl object to report convergence + solver_control.check(0, 0); + + if (solver_control.last_check() != SolverControl::success) + AssertThrow(false, + SolverControl::NoConvergence(solver_control.last_step(), + solver_control.last_value())); + } + + + + template + void + SolverDirectBase::solve( + const SparseMatrix &A, + Vector &x, + const Vector &b) + { + solver = + Amesos2::create::MatrixType, + MultiVector>(solver_type, + A.trilinos_rcp(), + x.trilinos_rcp(), + b.trilinos_rcp()); + do_solve(); + } + + + + template + SolverDirect::AdditionalData::AdditionalData( + const std::string &solver_name, + const bool output_solver_details) + : solver_name(solver_name) + , output_solver_details(output_solver_details) + {} + + + + template + SolverDirect::SolverDirect(SolverControl &cn, + const AdditionalData &ad) + : SolverDirectBase(cn, + ad.solver_name, + ad.output_solver_details) + {} + + + + template + void + SolverDirect::set_pararameter_list( + Teuchos::ParameterList ¶meter_list) + { + this->parameter_list.setParameters(parameter_list); + } + + + + template + SolverDirectKLU2::AdditionalData::AdditionalData( + const std::string &transpose_mode, + const bool symmetric_mode, + const bool equilibrate_matrix, + const std::string &column_permutation, + const std::string &iterative_refinement, + const bool output_solver_details) + : transpose_mode(transpose_mode) + , symmetric_mode(symmetric_mode) + , equilibrate_matrix(equilibrate_matrix) + , column_permutation(column_permutation) + , iterative_refinement(iterative_refinement) + , output_solver_details(output_solver_details) + {} + + + + template + SolverDirectKLU2::SolverDirectKLU2( + SolverControl &cn, + const AdditionalData &ad) + : SolverDirectBase(cn, + "KLU2", + ad.output_solver_details) + { + this->parameter_list = Teuchos::ParameterList("Amesos2"); + Teuchos::ParameterList klu2_params = this->parameter_list.sublist("KLU2"); + klu2_params.set("Trans", ad.transpose_mode); + klu2_params.set("Equil", ad.equilibrate_matrix); + klu2_params.set("IterRefine", ad.iterative_refinement); + klu2_params.set("SymmetricMode", ad.symmetric_mode); + klu2_params.set("ColPerm", ad.column_permutation); + } + + } // namespace TpetraWrappers +} // namespace LinearAlgebra + +DEAL_II_NAMESPACE_CLOSE + +# endif // DEAL_II_TRILINOS_WITH_AMESOS2 +#endif // DEAL_II_TRILINOS_WITH_TPETRA + +#endif diff --git a/source/lac/CMakeLists.txt b/source/lac/CMakeLists.txt index 207197951e..f2548f3924 100644 --- a/source/lac/CMakeLists.txt +++ b/source/lac/CMakeLists.txt @@ -144,6 +144,7 @@ if(DEAL_II_WITH_TRILINOS) trilinos_sparse_matrix.cc trilinos_sparsity_pattern.cc trilinos_tpetra_communication_pattern.cc + trilinos_tpetra_solver_direct.cc trilinos_tpetra_sparse_matrix.cc trilinos_tpetra_sparsity_pattern.cc trilinos_tpetra_vector.cc diff --git a/source/lac/trilinos_tpetra_solver_direct.cc b/source/lac/trilinos_tpetra_solver_direct.cc new file mode 100644 index 0000000000..512baf815b --- /dev/null +++ b/source/lac/trilinos_tpetra_solver_direct.cc @@ -0,0 +1,61 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2008 - 2024 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +#include "deal.II/lac/trilinos_tpetra_solver_direct.templates.h" + +#ifdef DEAL_II_TRILINOS_WITH_TPETRA +# ifdef DEAL_II_TRILINOS_WITH_AMESOS2 + +DEAL_II_NAMESPACE_OPEN + +# ifndef DOXYGEN +// explicit instantiations +namespace LinearAlgebra +{ + namespace TpetraWrappers + { +# ifdef HAVE_TPETRA_INST_FLOAT + template class SolverDirectBase; + template class SolverDirect; + template class SolverDirectKLU2; +# endif + +# ifdef HAVE_TPETRA_INST_DOUBLE + template class SolverDirectBase; + template class SolverDirect; + template class SolverDirectKLU2; +# endif +# ifdef DEAL_II_WITH_COMPLEX_VALUES +# ifdef HAVE_TPETRA_INST_COMPLEX_FLOAT + template class SolverDirectBase>; + template class SolverDirect>; + template class SolverDirectKLU2>; +# endif +# ifdef HAVE_TPETRA_INST_COMPLEX_DOUBLE + template class SolverDirectBase>; + template class SolverDirect>; + template class SolverDirectKLU2>; +# endif +# endif + +# endif + + } // namespace TpetraWrappers + +} // namespace LinearAlgebra +DEAL_II_NAMESPACE_CLOSE + +# endif // DEAL_II_TRILINOS_WITH_AMESOS2 +#endif // DEAL_II_TRILINOS_WITH_TPETRA diff --git a/tests/trilinos_tpetra/direct_solver_01.cc b/tests/trilinos_tpetra/direct_solver_01.cc new file mode 100644 index 0000000000..8889251123 --- /dev/null +++ b/tests/trilinos_tpetra/direct_solver_01.cc @@ -0,0 +1,463 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2013 - 2024 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +// tests Trilinos direct solvers on a 2D Poisson equation for linear elements + +#include "deal.II/base/logstream.h" +#include +#include +#include +#include + +#include + +#include +#include + +#include +#include + +#include +#include + +#include "deal.II/lac/exceptions.h" +#include +#include +#include +#include +#include +// #include +#include +#include + +#include + +#include "../tests.h" + + + +template +class Step4 +{ +public: + Step4(); + void + run(); + +private: + void + make_grid(); + void + setup_system(); + void + assemble_system(); + void + solve(); + + parallel::distributed::Triangulation triangulation; + FE_Q fe; + DoFHandler dof_handler; + + AffineConstraints constraints; + SparsityPattern sparsity_pattern; + + LinearAlgebra::TpetraWrappers::SparseMatrix system_matrix; + + LinearAlgebra::TpetraWrappers::Vector solution; + LinearAlgebra::TpetraWrappers::Vector system_rhs; + LinearAlgebra::TpetraWrappers::Vector system_rhs_two; +}; + + +template +class RightHandSide : public Function +{ +public: + RightHandSide() + : Function() + {} + + virtual double + value(const Point &p, const unsigned int component = 0) const; +}; + + +template +class RightHandSideTwo : public Function +{ +public: + RightHandSideTwo() + : 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 < dim; ++i) + return_value += 2 * std::pow(p[i], 2); + + return return_value; +} + + +template +double +RightHandSideTwo::value(const Point &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 +double +BoundaryValues::value(const Point &p, + const unsigned int /*component*/) const +{ + return p.square(); +} + + + +template +Step4::Step4() + : triangulation(MPI_COMM_WORLD, + typename Triangulation::MeshSmoothing( + Triangulation::smoothing_on_refinement | + Triangulation::smoothing_on_coarsening)) + , 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(); + + const IndexSet &locally_owned_dofs = dof_handler.locally_owned_dofs(); + const IndexSet locally_relevant_dofs = + DoFTools::extract_locally_relevant_dofs(dof_handler); + + + DynamicSparsityPattern dsp(dof_handler.n_dofs()); + DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints, false); + SparsityTools::distribute_sparsity_pattern(dsp, + locally_owned_dofs, + MPI_COMM_WORLD, + locally_relevant_dofs); + + system_matrix.reinit(locally_owned_dofs, + locally_owned_dofs, + dsp, + MPI_COMM_WORLD); + + solution.reinit(locally_relevant_dofs, MPI_COMM_WORLD); + + system_rhs.reinit(locally_owned_dofs, + locally_relevant_dofs, + MPI_COMM_WORLD, + true); + + system_rhs_two.reinit(locally_owned_dofs, + locally_relevant_dofs, + MPI_COMM_WORLD, + true); +} + + +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); + Vector cell_rhs_two(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) + { + if (cell->is_locally_owned()) + { + fe_values.reinit(cell); + cell_matrix = 0; + cell_rhs = 0; + cell_rhs_two = 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_rhs_two(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); + + constraints.distribute_local_to_global(cell_rhs_two, + local_dof_indices, + system_rhs_two); + } + } + system_matrix.compress(VectorOperation::add); + system_rhs.compress(VectorOperation::add); + system_rhs_two.compress(VectorOperation::add); +} + + + +template +void +Step4::solve() +{ + // TODO: implement precondition SSOR Wrappers + // // Compute 'reference' solution with CG solver and SSOR preconditioner + // TrilinosWrappers::PreconditionSSOR preconditioner; + // preconditioner.initialize(system_matrix); + LinearAlgebra::TpetraWrappers::Vector temp_solution(system_rhs); + temp_solution = 0; + SolverControl solver_control(1000, 1e-12); + SolverCG> solver( + solver_control); + + solver.solve(system_matrix, + temp_solution, + system_rhs, + PreconditionIdentity()); + + constraints.distribute(temp_solution); + solution = temp_solution; + + LinearAlgebra::TpetraWrappers::Vector output(temp_solution); + + // do CG solve for new rhs + temp_solution = 0; + solution = 0; + solver.solve(system_matrix, + temp_solution, + system_rhs_two, + PreconditionIdentity()); + + constraints.distribute(temp_solution); + solution = temp_solution; + + LinearAlgebra::TpetraWrappers::Vector output_two(temp_solution); + + // factorize matrix for direct solver + temp_solution = 0; + solution = 0; + + { + deallog.push("DirectKLU2"); + LinearAlgebra::TpetraWrappers::SolverDirect::AdditionalData data( + "KLU2", true); + LinearAlgebra::TpetraWrappers::SolverDirect direct_solver( + solver_control, data); + + Teuchos::ParameterList amesos2_params("Amesos2"); + Teuchos::ParameterList klu2_params = amesos2_params.sublist("KLU2"); + klu2_params.set("ColPerm", "NATURAL"); + + direct_solver.set_pararameter_list(amesos2_params); + + direct_solver.initialize(system_matrix); + + // do solve 1 + direct_solver.solve(temp_solution, system_rhs); + + constraints.distribute(temp_solution); + solution = temp_solution; + + // calculate l2 errors + output.add(-1.0, temp_solution); + + const double local_error = output.l2_norm(); + const double global_error = + std::sqrt(Utilities::MPI::sum(local_error * local_error, MPI_COMM_WORLD)); + + deallog << "Norm of error in direct solve 1: " << global_error << std::endl; + + // do solve 2 without refactorizing + temp_solution = 0; + solution = 0; + direct_solver.solve(temp_solution, system_rhs_two); + + constraints.distribute(temp_solution); + solution = temp_solution; + + // calculate l2 errors + output_two.add(-1.0, temp_solution); + + const double local_error_two = output_two.l2_norm(); + const double global_error_two = std::sqrt( + Utilities::MPI::sum(local_error_two * local_error_two, MPI_COMM_WORLD)); + + deallog << "Norm of error in direct solve 2: " << global_error_two + << std::endl; + + deallog.pop(); + } + { + // Ensure that the exception for an unavailable solver is thrown. + // MEASLES is a non-existent solver so it throws the exception + // independent of the Amesos2 configuration. + deallog.push("DirectMEASLES"); + LinearAlgebra::TpetraWrappers::SolverDirect::AdditionalData data( + "MEASLES", true); + try + { + LinearAlgebra::TpetraWrappers::SolverDirect direct_solver( + solver_control, data); + } + catch (dealii::ExceptionBase &exc) + { + deallog << "Error: " << std::endl; + exc.print_info(deallog.get_file_stream()); + } + } +} + + + +template +void +Step4::run() +{ + make_grid(); + setup_system(); + assemble_system(); + solve(); +} + + +int +main(int argc, char **argv) +{ + initlog(); + + Utilities::MPI::MPI_InitFinalize mpi_initialization( + argc, argv, testing_max_num_threads()); + + try + { + Step4<2> test; + test.run(); + } + catch (const 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_tpetra/direct_solver_01.with_p4est=true.with_trilinos_with_amesos2=true.mpirun=1.output b/tests/trilinos_tpetra/direct_solver_01.with_p4est=true.with_trilinos_with_amesos2=true.mpirun=1.output new file mode 100644 index 0000000000..1facc1d50a --- /dev/null +++ b/tests/trilinos_tpetra/direct_solver_01.with_p4est=true.with_trilinos_with_amesos2=true.mpirun=1.output @@ -0,0 +1,18 @@ + +DEAL:cg::Starting value 21.8027 +DEAL:cg::Convergence step 109 value 0 +DEAL:cg::Starting value 0.0940512 +DEAL:cg::Convergence step 98 value 0 +DEAL:DirectKLU2::Starting value 0 +DEAL:DirectKLU2::Convergence step 0 value 0 +DEAL:DirectKLU2::Norm of error in direct solve 1: 0 +DEAL:DirectKLU2::Starting value 0 +DEAL:DirectKLU2::Convergence step 0 value 0 +DEAL:DirectKLU2::Norm of error in direct solve 2: 0 +DEAL:DirectMEASLES::Error: + You tried to select the solver type +but this solver is not supported by Trilinos/Amesos2 +due to one of the following reasons: +* This solver does not exist +* This solver is not (yet) supported by Trilinos/Amesos2 +* Trilinos/Amesos2 was not configured for its use. \ No newline at end of file diff --git a/tests/trilinos_tpetra/direct_solver_02.cc b/tests/trilinos_tpetra/direct_solver_02.cc new file mode 100644 index 0000000000..6b39431f8a --- /dev/null +++ b/tests/trilinos_tpetra/direct_solver_02.cc @@ -0,0 +1,433 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2013 - 2024 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +// tests Trilinos direct solvers on a 2D Poisson equation for linear elements + +#include +#include +#include +#include + +#include + +#include +#include + +#include +#include + +#include +#include + +#include +#include +#include +#include +#include +// #include +#include +#include + +#include + +#include "../tests.h" + + + +template +class Step4 +{ +public: + Step4(); + void + run(); + +private: + void + make_grid(); + void + setup_system(); + void + assemble_system(); + void + solve(); + + parallel::distributed::Triangulation triangulation; + FE_Q fe; + DoFHandler dof_handler; + + AffineConstraints constraints; + SparsityPattern sparsity_pattern; + + LinearAlgebra::TpetraWrappers::SparseMatrix system_matrix; + + LinearAlgebra::TpetraWrappers::Vector solution; + LinearAlgebra::TpetraWrappers::Vector system_rhs; + LinearAlgebra::TpetraWrappers::Vector system_rhs_two; +}; + + +template +class RightHandSide : public Function +{ +public: + RightHandSide() + : Function() + {} + + virtual double + value(const Point &p, const unsigned int component = 0) const; +}; + + +template +class RightHandSideTwo : public Function +{ +public: + RightHandSideTwo() + : 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 < dim; ++i) + return_value += 2 * std::pow(p[i], 2); + + return return_value; +} + + +template +double +RightHandSideTwo::value(const Point &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 +double +BoundaryValues::value(const Point &p, + const unsigned int /*component*/) const +{ + return p.square(); +} + + + +template +Step4::Step4() + : triangulation(MPI_COMM_WORLD, + typename Triangulation::MeshSmoothing( + Triangulation::smoothing_on_refinement | + Triangulation::smoothing_on_coarsening)) + , 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(); + + const IndexSet &locally_owned_dofs = dof_handler.locally_owned_dofs(); + const IndexSet locally_relevant_dofs = + DoFTools::extract_locally_relevant_dofs(dof_handler); + + + DynamicSparsityPattern dsp(dof_handler.n_dofs()); + DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints, false); + SparsityTools::distribute_sparsity_pattern(dsp, + locally_owned_dofs, + MPI_COMM_WORLD, + locally_relevant_dofs); + + system_matrix.reinit(locally_owned_dofs, + locally_owned_dofs, + dsp, + MPI_COMM_WORLD); + + solution.reinit(locally_relevant_dofs, MPI_COMM_WORLD); + + system_rhs.reinit(locally_owned_dofs, + locally_relevant_dofs, + MPI_COMM_WORLD, + true); + + system_rhs_two.reinit(locally_owned_dofs, + locally_relevant_dofs, + MPI_COMM_WORLD, + true); +} + + +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); + Vector cell_rhs_two(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) + { + if (cell->is_locally_owned()) + { + fe_values.reinit(cell); + cell_matrix = 0; + cell_rhs = 0; + cell_rhs_two = 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_rhs_two(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); + + constraints.distribute_local_to_global(cell_rhs_two, + local_dof_indices, + system_rhs_two); + } + } + system_matrix.compress(VectorOperation::add); + system_rhs.compress(VectorOperation::add); + system_rhs_two.compress(VectorOperation::add); +} + + + +template +void +Step4::solve() +{ + // TODO: implement precondition SSOR Wrappers + // // Compute 'reference' solution with CG solver and SSOR preconditioner + // TrilinosWrappers::PreconditionSSOR preconditioner; + // preconditioner.initialize(system_matrix); + LinearAlgebra::TpetraWrappers::Vector temp_solution(system_rhs); + temp_solution = 0; + SolverControl solver_control(1000, 1e-12); + SolverCG> solver( + solver_control); + + solver.solve(system_matrix, + temp_solution, + system_rhs, + PreconditionIdentity()); + + constraints.distribute(temp_solution); + solution = temp_solution; + + LinearAlgebra::TpetraWrappers::Vector output(temp_solution); + + // do CG solve for new rhs + temp_solution = 0; + solution = 0; + solver.solve(system_matrix, + temp_solution, + system_rhs_two, + PreconditionIdentity()); + + constraints.distribute(temp_solution); + solution = temp_solution; + + LinearAlgebra::TpetraWrappers::Vector output_two(temp_solution); + + // factorize matrix for direct solver + temp_solution = 0; + solution = 0; + + deallog.push("DirectKLU2"); + LinearAlgebra::TpetraWrappers::SolverDirectKLU2::AdditionalData data; + LinearAlgebra::TpetraWrappers::SolverDirectKLU2 direct_solver( + solver_control, data); + direct_solver.initialize(system_matrix); + + // do solve 1 + direct_solver.solve(temp_solution, system_rhs); + + constraints.distribute(temp_solution); + solution = temp_solution; + + // calculate l2 errors + output.add(-1.0, temp_solution); + + const double local_error = output.l2_norm(); + const double global_error = + std::sqrt(Utilities::MPI::sum(local_error * local_error, MPI_COMM_WORLD)); + + deallog << "Norm of error in direct solve 1: " << global_error << std::endl; + + // do solve 2 without refactorizing + temp_solution = 0; + solution = 0; + direct_solver.solve(temp_solution, system_rhs_two); + + constraints.distribute(temp_solution); + solution = temp_solution; + + // calculate l2 errors + output_two.add(-1.0, temp_solution); + + const double local_error_two = output_two.l2_norm(); + const double global_error_two = std::sqrt( + Utilities::MPI::sum(local_error_two * local_error_two, MPI_COMM_WORLD)); + + deallog << "Norm of error in direct solve 2: " << global_error_two + << std::endl; + + deallog.pop(); +} + + + +template +void +Step4::run() +{ + make_grid(); + setup_system(); + assemble_system(); + solve(); +} + + +int +main(int argc, char **argv) +{ + initlog(); + + Utilities::MPI::MPI_InitFinalize mpi_initialization( + argc, argv, testing_max_num_threads()); + + try + { + Step4<2> test; + test.run(); + } + catch (const 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_tpetra/direct_solver_02.with_p4est=true.with_trilinos_with_amesos2=true.mpirun=1.output b/tests/trilinos_tpetra/direct_solver_02.with_p4est=true.with_trilinos_with_amesos2=true.mpirun=1.output new file mode 100644 index 0000000000..2f8d30713d --- /dev/null +++ b/tests/trilinos_tpetra/direct_solver_02.with_p4est=true.with_trilinos_with_amesos2=true.mpirun=1.output @@ -0,0 +1,11 @@ + +DEAL:cg::Starting value 21.8027 +DEAL:cg::Convergence step 109 value 0 +DEAL:cg::Starting value 0.0940512 +DEAL:cg::Convergence step 98 value 0 +DEAL:DirectKLU2::Starting value 0 +DEAL:DirectKLU2::Convergence step 0 value 0 +DEAL:DirectKLU2::Norm of error in direct solve 1: 0 +DEAL:DirectKLU2::Starting value 0 +DEAL:DirectKLU2::Convergence step 0 value 0 +DEAL:DirectKLU2::Norm of error in direct solve 2: 0