]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add initial support for Amesos2 16602/head
authorJan Philipp Thiele <thiele@wias-berlin.de>
Mon, 12 Feb 2024 19:53:17 +0000 (20:53 +0100)
committerJan Philipp Thiele <thiele@wias-berlin.de>
Sat, 17 Feb 2024 12:15:02 +0000 (13:15 +0100)
cmake/configure/configure_20_trilinos.cmake
doc/doxygen/options.dox.in
include/deal.II/base/config.h.in
include/deal.II/lac/trilinos_tpetra_solver_direct.h [new file with mode: 0644]
include/deal.II/lac/trilinos_tpetra_solver_direct.templates.h [new file with mode: 0644]
source/lac/CMakeLists.txt
source/lac/trilinos_tpetra_solver_direct.cc [new file with mode: 0644]
tests/trilinos_tpetra/direct_solver_01.cc [new file with mode: 0644]
tests/trilinos_tpetra/direct_solver_01.with_p4est=true.with_trilinos_with_amesos2=true.mpirun=1.output [new file with mode: 0644]
tests/trilinos_tpetra/direct_solver_02.cc [new file with mode: 0644]
tests/trilinos_tpetra/direct_solver_02.with_p4est=true.with_trilinos_with_amesos2=true.mpirun=1.output [new file with mode: 0644]

index f8ba7e57d71a8d1adf2c1c6f39169263adca92c0..d9f3c280f775d7be274c9f62481ed936a0d47ab8 100644 (file)
@@ -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 
 )
 
 #
index 0e8058e97005ac152e79ff66398f7145f11c15ff..82d336b9eeffc9b31dda45c6bc859e85bfa9b95a 100644 (file)
@@ -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 \
index c9943b3efea1954c76f824d6e5d509d244ec3fa5..25367eab2ab173265316e99ff9acdf9687585d01 100644 (file)
 
 /* 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 (file)
index 0000000..b422f43
--- /dev/null
@@ -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 <deal.II/base/config.h>
+
+#include "deal.II/grid/grid_generator.h"
+
+#include <string>
+
+#ifdef DEAL_II_TRILINOS_WITH_TPETRA
+#  ifdef DEAL_II_TRILINOS_WITH_AMESOS2
+
+#    include "deal.II/base/types.h"
+
+#    include <deal.II/lac/la_parallel_vector.h>
+#    include <deal.II/lac/solver_control.h>
+#    include <deal.II/lac/trilinos_tpetra_sparse_matrix.h>
+
+#    include <Amesos2.hpp>
+#    include <Teuchos_ConfigDefs.hpp>
+#    include <Teuchos_ParameterList.hpp>
+#    include <Teuchos_RCPDecl.hpp>
+
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace LinearAlgebra
+{
+  namespace TpetraWrappers
+  {
+    // forward declarations
+#    ifndef DOXYGEN
+    template <typename Number, typename MemorySpace>
+    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 <typename Number, typename MemorySpace = dealii::MemorySpace::Host>
+    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<Number, int, size_type, NodeType>;
+
+      /**
+       * Typedef for the MultiVector type
+       */
+      using MultiVector = Tpetra::MultiVector<Number, int, size_type, NodeType>;
+
+      /**
+       * Destructor.
+       */
+      virtual ~SolverDirectBase() = default;
+
+      /**
+       * Initializes the direct solver for the matrix <tt>A</tt> 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<Number, MemorySpace> &A);
+
+      /**
+       * Solve the linear system <tt>Ax=b</tt> based on the
+       * package set in initialize(). Note the matrix is not refactorized during
+       * this call.
+       */
+      void
+      solve(Vector<Number, MemorySpace>       &x,
+            const Vector<Number, MemorySpace> &b);
+
+      /**
+       * Solve the linear system <tt>Ax=b</tt>. 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<Number, MemorySpace> &A,
+            Vector<Number, MemorySpace>             &x,
+            const Vector<Number, MemorySpace>       &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<typename SparseMatrix<Number, MemorySpace>::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 <typename Number, typename MemorySpace = dealii::MemorySpace::Host>
+    class SolverDirect : public SolverDirectBase<Number, MemorySpace>
+    {
+    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:
+         * <ul>
+         * <li>  "Basker" </li>
+         * <li>  "Cholmod" </li>
+         * <li>  "cuSOLVER" </li>
+         * <li>  "KLU2" </li>
+         * <li>  "LAPACK" </li>
+         * <li>  "MUMPS" </li>
+         * <li>  "PardisoMKL" </li>
+         * <li>  "ShyLUBasker" </li>
+         * <li>  "SuperLU" </li>
+         * <li>  "SuperLU_DIST" </li>
+         * <li>  "SuperLU_MT" </li>
+         * <li>  "Tacho" </li>
+         * <li>  "UMFPACK" </li>
+         * </ul>
+         * 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 &parameter_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 <typename Number, typename MemorySpace = dealii::MemorySpace::Host>
+    class SolverDirectKLU2 : public SolverDirectBase<Number, MemorySpace>
+    {
+    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 (file)
index 0000000..238f7ef
--- /dev/null
@@ -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 <deal.II/base/config.h>
+
+#include "deal.II/lac/solver_control.h"
+
+#include <string>
+
+
+#ifdef DEAL_II_TRILINOS_WITH_TPETRA
+#  ifdef DEAL_II_TRILINOS_WITH_AMESOS2
+
+#    include "deal.II/base/template_constraints.h"
+#    include <deal.II/base/conditional_ostream.h>
+
+#    include <deal.II/lac/trilinos_tpetra_solver_direct.h>
+
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace LinearAlgebra
+{
+  namespace TpetraWrappers
+  {
+    /* ---------------------- SolverDirectBase ------------------------ */
+
+    template <typename Number, typename MemorySpace>
+    SolverDirectBase<Number, MemorySpace>::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 <typename Number, typename MemorySpace>
+    SolverControl &
+    SolverDirectBase<Number, MemorySpace>::control() const
+    {
+      return solver_control;
+    }
+
+
+
+    template <typename Number, typename MemorySpace>
+    void
+    SolverDirectBase<Number, MemorySpace>::initialize(
+      const SparseMatrix<Number, MemorySpace> &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<typename SparseMatrix<Number, MemorySpace>::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 <typename Number, typename MemorySpace>
+    void
+    SolverDirectBase<Number, MemorySpace>::solve(
+      Vector<Number, MemorySpace>       &x,
+      const Vector<Number, MemorySpace> &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 <typename Number, typename MemorySpace>
+    void
+    SolverDirectBase<Number, MemorySpace>::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 <typename Number, typename MemorySpace>
+    void
+    SolverDirectBase<Number, MemorySpace>::solve(
+      const SparseMatrix<Number, MemorySpace> &A,
+      Vector<Number, MemorySpace>             &x,
+      const Vector<Number, MemorySpace>       &b)
+    {
+      solver =
+        Amesos2::create<typename SparseMatrix<Number, MemorySpace>::MatrixType,
+                        MultiVector>(solver_type,
+                                     A.trilinos_rcp(),
+                                     x.trilinos_rcp(),
+                                     b.trilinos_rcp());
+      do_solve();
+    }
+
+
+
+    template <typename Number, typename MemorySpace>
+    SolverDirect<Number, MemorySpace>::AdditionalData::AdditionalData(
+      const std::string &solver_name,
+      const bool         output_solver_details)
+      : solver_name(solver_name)
+      , output_solver_details(output_solver_details)
+    {}
+
+
+
+    template <typename Number, typename MemorySpace>
+    SolverDirect<Number, MemorySpace>::SolverDirect(SolverControl        &cn,
+                                                    const AdditionalData &ad)
+      : SolverDirectBase<Number, MemorySpace>(cn,
+                                              ad.solver_name,
+                                              ad.output_solver_details)
+    {}
+
+
+
+    template <typename Number, typename MemorySpace>
+    void
+    SolverDirect<Number, MemorySpace>::set_pararameter_list(
+      Teuchos::ParameterList &parameter_list)
+    {
+      this->parameter_list.setParameters(parameter_list);
+    }
+
+
+
+    template <typename Number, typename MemorySpace>
+    SolverDirectKLU2<Number, MemorySpace>::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 <typename Number, typename MemorySpace>
+    SolverDirectKLU2<Number, MemorySpace>::SolverDirectKLU2(
+      SolverControl        &cn,
+      const AdditionalData &ad)
+      : SolverDirectBase<Number, MemorySpace>(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
index 207197951e4044817cae680c4134bc248017f563..f2548f3924abb1f1661d998fed3adc5cc7eb4b73 100644 (file)
@@ -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 (file)
index 0000000..512baf8
--- /dev/null
@@ -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<float>;
+    template class SolverDirect<float>;
+    template class SolverDirectKLU2<float>;
+#      endif
+
+#      ifdef HAVE_TPETRA_INST_DOUBLE
+    template class SolverDirectBase<double>;
+    template class SolverDirect<double>;
+    template class SolverDirectKLU2<double>;
+#      endif
+#      ifdef DEAL_II_WITH_COMPLEX_VALUES
+#        ifdef HAVE_TPETRA_INST_COMPLEX_FLOAT
+    template class SolverDirectBase<complex<float>>;
+    template class SolverDirect<complex<float>>;
+    template class SolverDirectKLU2<complex<float>>;
+#        endif
+#        ifdef HAVE_TPETRA_INST_COMPLEX_DOUBLE
+    template class SolverDirectBase<complex<double>>;
+    template class SolverDirect<complex<double>>;
+    template class SolverDirectKLU2<complex<double>>;
+#        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 (file)
index 0000000..8889251
--- /dev/null
@@ -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 <deal.II/base/conditional_ostream.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/index_set.h>
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include "deal.II/lac/exceptions.h"
+#include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/sparsity_tools.h>
+// #include <deal.II/lac/trilinos_precondition.h>
+#include <deal.II/lac/trilinos_tpetra_solver_direct.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include "../tests.h"
+
+
+
+template <int dim>
+class Step4
+{
+public:
+  Step4();
+  void
+  run();
+
+private:
+  void
+  make_grid();
+  void
+  setup_system();
+  void
+  assemble_system();
+  void
+  solve();
+
+  parallel::distributed::Triangulation<dim> triangulation;
+  FE_Q<dim>                                 fe;
+  DoFHandler<dim>                           dof_handler;
+
+  AffineConstraints<double> constraints;
+  SparsityPattern           sparsity_pattern;
+
+  LinearAlgebra::TpetraWrappers::SparseMatrix<double> system_matrix;
+
+  LinearAlgebra::TpetraWrappers::Vector<double> solution;
+  LinearAlgebra::TpetraWrappers::Vector<double> system_rhs;
+  LinearAlgebra::TpetraWrappers::Vector<double> system_rhs_two;
+};
+
+
+template <int dim>
+class RightHandSide : public Function<dim>
+{
+public:
+  RightHandSide()
+    : Function<dim>()
+  {}
+
+  virtual double
+  value(const Point<dim> &p, const unsigned int component = 0) const;
+};
+
+
+template <int dim>
+class RightHandSideTwo : public Function<dim>
+{
+public:
+  RightHandSideTwo()
+    : Function<dim>()
+  {}
+
+  virtual double
+  value(const Point<dim> &p, const unsigned int component = 0) const;
+};
+
+
+
+template <int dim>
+class BoundaryValues : public Function<dim>
+{
+public:
+  BoundaryValues()
+    : Function<dim>()
+  {}
+
+  virtual double
+  value(const Point<dim> &p, const unsigned int component = 0) const;
+};
+
+
+
+template <int dim>
+double
+RightHandSide<dim>::value(const Point<dim> &p,
+                          const unsigned int /*component*/) const
+{
+  double return_value = 0;
+  for (unsigned int i = 0; i < dim; ++i)
+    return_value += 2 * std::pow(p[i], 2);
+
+  return return_value;
+}
+
+
+template <int dim>
+double
+RightHandSideTwo<dim>::value(const Point<dim> &p,
+                             const unsigned int /*component*/) const
+{
+  double return_value = 0;
+  for (unsigned int i = 0; i < dim; ++i)
+    return_value += 4 * std::pow(p(i), 4);
+
+  return return_value;
+}
+
+
+template <int dim>
+double
+BoundaryValues<dim>::value(const Point<dim> &p,
+                           const unsigned int /*component*/) const
+{
+  return p.square();
+}
+
+
+
+template <int dim>
+Step4<dim>::Step4()
+  : triangulation(MPI_COMM_WORLD,
+                  typename Triangulation<dim>::MeshSmoothing(
+                    Triangulation<dim>::smoothing_on_refinement |
+                    Triangulation<dim>::smoothing_on_coarsening))
+  , fe(1)
+  , dof_handler(triangulation)
+{}
+
+
+template <int dim>
+void
+Step4<dim>::make_grid()
+{
+  GridGenerator::hyper_cube(triangulation, -1, 1);
+  triangulation.refine_global(6);
+}
+
+
+
+template <int dim>
+void
+Step4<dim>::setup_system()
+{
+  dof_handler.distribute_dofs(fe);
+
+  constraints.clear();
+  std::map<unsigned int, double> boundary_values;
+  VectorTools::interpolate_boundary_values(dof_handler,
+                                           0,
+                                           BoundaryValues<dim>(),
+                                           constraints);
+  constraints.close();
+
+  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 <int dim>
+void
+Step4<dim>::assemble_system()
+{
+  QGauss<dim> quadrature_formula(fe.degree + 1);
+
+  const RightHandSide<dim> right_hand_side;
+
+  FEValues<dim> fe_values(fe,
+                          quadrature_formula,
+                          update_values | update_gradients |
+                            update_quadrature_points | update_JxW_values);
+
+  const unsigned int dofs_per_cell = fe.dofs_per_cell;
+  const unsigned int n_q_points    = quadrature_formula.size();
+
+  FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
+  Vector<double>     cell_rhs(dofs_per_cell);
+  Vector<double>     cell_rhs_two(dofs_per_cell);
+
+  std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
+
+  typename DoFHandler<dim>::active_cell_iterator cell =
+                                                   dof_handler.begin_active(),
+                                                 endc = dof_handler.end();
+
+  for (; cell != endc; ++cell)
+    {
+      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 <int dim>
+void
+Step4<dim>::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<double> temp_solution(system_rhs);
+  temp_solution = 0;
+  SolverControl solver_control(1000, 1e-12);
+  SolverCG<LinearAlgebra::TpetraWrappers::Vector<double>> solver(
+    solver_control);
+
+  solver.solve(system_matrix,
+               temp_solution,
+               system_rhs,
+               PreconditionIdentity());
+
+  constraints.distribute(temp_solution);
+  solution = temp_solution;
+
+  LinearAlgebra::TpetraWrappers::Vector<double> 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<double> output_two(temp_solution);
+
+  // factorize matrix for direct solver
+  temp_solution = 0;
+  solution      = 0;
+
+  {
+    deallog.push("DirectKLU2");
+    LinearAlgebra::TpetraWrappers::SolverDirect<double>::AdditionalData data(
+      "KLU2", true);
+    LinearAlgebra::TpetraWrappers::SolverDirect<double> 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<double>::AdditionalData data(
+      "MEASLES", true);
+    try
+      {
+        LinearAlgebra::TpetraWrappers::SolverDirect<double> direct_solver(
+          solver_control, data);
+      }
+    catch (dealii::ExceptionBase &exc)
+      {
+        deallog << "Error: " << std::endl;
+        exc.print_info(deallog.get_file_stream());
+      }
+  }
+}
+
+
+
+template <int dim>
+void
+Step4<dim>::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 (file)
index 0000000..1facc1d
--- /dev/null
@@ -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 <MEASLES>
+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 (file)
index 0000000..6b39431
--- /dev/null
@@ -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 <deal.II/base/conditional_ostream.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/index_set.h>
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/sparsity_tools.h>
+// #include <deal.II/lac/trilinos_precondition.h>
+#include <deal.II/lac/trilinos_tpetra_solver_direct.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include "../tests.h"
+
+
+
+template <int dim>
+class Step4
+{
+public:
+  Step4();
+  void
+  run();
+
+private:
+  void
+  make_grid();
+  void
+  setup_system();
+  void
+  assemble_system();
+  void
+  solve();
+
+  parallel::distributed::Triangulation<dim> triangulation;
+  FE_Q<dim>                                 fe;
+  DoFHandler<dim>                           dof_handler;
+
+  AffineConstraints<double> constraints;
+  SparsityPattern           sparsity_pattern;
+
+  LinearAlgebra::TpetraWrappers::SparseMatrix<double> system_matrix;
+
+  LinearAlgebra::TpetraWrappers::Vector<double> solution;
+  LinearAlgebra::TpetraWrappers::Vector<double> system_rhs;
+  LinearAlgebra::TpetraWrappers::Vector<double> system_rhs_two;
+};
+
+
+template <int dim>
+class RightHandSide : public Function<dim>
+{
+public:
+  RightHandSide()
+    : Function<dim>()
+  {}
+
+  virtual double
+  value(const Point<dim> &p, const unsigned int component = 0) const;
+};
+
+
+template <int dim>
+class RightHandSideTwo : public Function<dim>
+{
+public:
+  RightHandSideTwo()
+    : Function<dim>()
+  {}
+
+  virtual double
+  value(const Point<dim> &p, const unsigned int component = 0) const;
+};
+
+
+
+template <int dim>
+class BoundaryValues : public Function<dim>
+{
+public:
+  BoundaryValues()
+    : Function<dim>()
+  {}
+
+  virtual double
+  value(const Point<dim> &p, const unsigned int component = 0) const;
+};
+
+
+
+template <int dim>
+double
+RightHandSide<dim>::value(const Point<dim> &p,
+                          const unsigned int /*component*/) const
+{
+  double return_value = 0;
+  for (unsigned int i = 0; i < dim; ++i)
+    return_value += 2 * std::pow(p[i], 2);
+
+  return return_value;
+}
+
+
+template <int dim>
+double
+RightHandSideTwo<dim>::value(const Point<dim> &p,
+                             const unsigned int /*component*/) const
+{
+  double return_value = 0;
+  for (unsigned int i = 0; i < dim; ++i)
+    return_value += 4 * std::pow(p(i), 4);
+
+  return return_value;
+}
+
+
+template <int dim>
+double
+BoundaryValues<dim>::value(const Point<dim> &p,
+                           const unsigned int /*component*/) const
+{
+  return p.square();
+}
+
+
+
+template <int dim>
+Step4<dim>::Step4()
+  : triangulation(MPI_COMM_WORLD,
+                  typename Triangulation<dim>::MeshSmoothing(
+                    Triangulation<dim>::smoothing_on_refinement |
+                    Triangulation<dim>::smoothing_on_coarsening))
+  , fe(1)
+  , dof_handler(triangulation)
+{}
+
+
+template <int dim>
+void
+Step4<dim>::make_grid()
+{
+  GridGenerator::hyper_cube(triangulation, -1, 1);
+  triangulation.refine_global(6);
+}
+
+
+
+template <int dim>
+void
+Step4<dim>::setup_system()
+{
+  dof_handler.distribute_dofs(fe);
+
+  constraints.clear();
+  std::map<unsigned int, double> boundary_values;
+  VectorTools::interpolate_boundary_values(dof_handler,
+                                           0,
+                                           BoundaryValues<dim>(),
+                                           constraints);
+  constraints.close();
+
+  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 <int dim>
+void
+Step4<dim>::assemble_system()
+{
+  QGauss<dim> quadrature_formula(fe.degree + 1);
+
+  const RightHandSide<dim> right_hand_side;
+
+  FEValues<dim> fe_values(fe,
+                          quadrature_formula,
+                          update_values | update_gradients |
+                            update_quadrature_points | update_JxW_values);
+
+  const unsigned int dofs_per_cell = fe.dofs_per_cell;
+  const unsigned int n_q_points    = quadrature_formula.size();
+
+  FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
+  Vector<double>     cell_rhs(dofs_per_cell);
+  Vector<double>     cell_rhs_two(dofs_per_cell);
+
+  std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
+
+  typename DoFHandler<dim>::active_cell_iterator cell =
+                                                   dof_handler.begin_active(),
+                                                 endc = dof_handler.end();
+
+  for (; cell != endc; ++cell)
+    {
+      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 <int dim>
+void
+Step4<dim>::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<double> temp_solution(system_rhs);
+  temp_solution = 0;
+  SolverControl solver_control(1000, 1e-12);
+  SolverCG<LinearAlgebra::TpetraWrappers::Vector<double>> solver(
+    solver_control);
+
+  solver.solve(system_matrix,
+               temp_solution,
+               system_rhs,
+               PreconditionIdentity());
+
+  constraints.distribute(temp_solution);
+  solution = temp_solution;
+
+  LinearAlgebra::TpetraWrappers::Vector<double> 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<double> output_two(temp_solution);
+
+  // factorize matrix for direct solver
+  temp_solution = 0;
+  solution      = 0;
+
+  deallog.push("DirectKLU2");
+  LinearAlgebra::TpetraWrappers::SolverDirectKLU2<double>::AdditionalData data;
+  LinearAlgebra::TpetraWrappers::SolverDirectKLU2<double> 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 <int dim>
+void
+Step4<dim>::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 (file)
index 0000000..2f8d307
--- /dev/null
@@ -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

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

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.