# 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
)
#
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 \
/* 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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 ¶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 <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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 ¶meter_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
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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+ };
+}
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+ };
+}
--- /dev/null
+
+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
projects / dealii.git / commitdiff