Line data    Source code

       1             : // The libMesh Finite Element Library.
       2             : // Copyright (C) 2002-2026 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
       3             : 
       4             : // This library is free software; you can redistribute it and/or
       5             : // modify it under the terms of the GNU Lesser General Public
       6             : // License as published by the Free Software Foundation; either
       7             : // version 2.1 of the License, or (at your option) any later version.
       8             : 
       9             : // This library is distributed in the hope that it will be useful,
      10             : // but WITHOUT ANY WARRANTY; without even the implied warranty of
      11             : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
      12             : // Lesser General Public License for more details.
      13             : 
      14             : // You should have received a copy of the GNU Lesser General Public
      15             : // License along with this library; if not, write to the Free Software
      16             : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
      17             : 
      18             : 
      19             : 
      20             : // Local Includes
      21             : #include "libmesh/libmesh_logging.h"
      22             : #include "libmesh/linear_solver.h"
      23             : #include "libmesh/laspack_linear_solver.h"
      24             : #include "libmesh/eigen_sparse_linear_solver.h"
      25             : #include "libmesh/petsc_linear_solver.h"
      26             : #include "libmesh/trilinos_aztec_linear_solver.h"
      27             : #include "libmesh/preconditioner.h"
      28             : #include "libmesh/sparse_matrix.h"
      29             : #include "libmesh/enum_to_string.h"
      30             : #include "libmesh/solver_configuration.h"
      31             : #include "libmesh/enum_solver_package.h"
      32             : #include "libmesh/enum_preconditioner_type.h"
      33             : #include "libmesh/enum_solver_type.h"
      34             : 
      35             : // C++ Includes
      36             : #include <memory>
      37             : 
      38             : namespace libMesh
      39             : {
      40             : 
      41             : //------------------------------------------------------------------
      42             : // LinearSolver members
      43             : template <typename T>
      44       30622 : LinearSolver<T>::LinearSolver (const libMesh::Parallel::Communicator & comm_in) :
      45             :   ParallelObject       (comm_in),
      46       28914 :   _solver_type         (GMRES),
      47       28914 :   _preconditioner_type (ILU_PRECOND),
      48       28914 :   _is_initialized      (false),
      49       28914 :   _preconditioner      (nullptr),
      50       28914 :   same_preconditioner  (false),
      51       30622 :   _solver_configuration(nullptr)
      52             : {
      53       30622 : }
      54             : 
      55             : 
      56             : 
      57             : template <typename T>
      58             : std::unique_ptr<LinearSolver<T>>
      59       30409 : LinearSolver<T>::build(const libMesh::Parallel::Communicator & comm,
      60             :                        const SolverPackage solver_package)
      61             : {
      62             :   // Avoid unused parameter warnings when no solver packages are enabled.
      63         848 :   libmesh_ignore(comm);
      64             : 
      65             :   // Build the appropriate solver
      66       30409 :   switch (solver_package)
      67             :     {
      68             : #ifdef LIBMESH_HAVE_LASPACK
      69           0 :     case LASPACK_SOLVERS:
      70           0 :       return std::make_unique<LaspackLinearSolver<T>>(comm);
      71             : #endif
      72             : 
      73             : 
      74             : #ifdef LIBMESH_HAVE_PETSC
      75       30100 :     case PETSC_SOLVERS:
      76       30100 :       return std::make_unique<PetscLinearSolver<T>>(comm);
      77             : #endif
      78             : 
      79             : 
      80             : #ifdef LIBMESH_TRILINOS_HAVE_AZTECOO
      81             :     case TRILINOS_SOLVERS:
      82             :       return std::make_unique<AztecLinearSolver<T>>(comm);
      83             : #endif
      84             : 
      85             : 
      86             : #ifdef LIBMESH_HAVE_EIGEN
      87         309 :     case EIGEN_SOLVERS:
      88         309 :       return std::make_unique<EigenSparseLinearSolver<T>>(comm);
      89             : #endif
      90             : 
      91           0 :     default:
      92           0 :       libmesh_error_msg("ERROR:  Unrecognized solver package: " << solver_package);
      93             :     }
      94             : 
      95             :   return std::unique_ptr<LinearSolver<T>>();
      96             : }
      97             : 
      98             : template <typename T>
      99             : PreconditionerType
     100           0 : LinearSolver<T>::preconditioner_type () const
     101             : {
     102           0 :   if (_preconditioner)
     103           0 :     return _preconditioner->type();
     104             : 
     105           0 :   return _preconditioner_type;
     106             : }
     107             : 
     108             : template <typename T>
     109             : void
     110         994 : LinearSolver<T>::set_preconditioner_type (const PreconditionerType pct)
     111             : {
     112         994 :   if (_preconditioner)
     113           0 :     _preconditioner->set_type(pct);
     114             :   else
     115         994 :     _preconditioner_type = pct;
     116         994 : }
     117             : 
     118             : template <typename T>
     119             : void
     120         350 : LinearSolver<T>::attach_preconditioner(Preconditioner<T> * preconditioner)
     121             : {
     122         350 :   libmesh_error_msg_if(this->_is_initialized,
     123             :                        "Preconditioner must be attached before the solver is initialized!");
     124             : 
     125         350 :   _preconditioner_type = SHELL_PRECOND;
     126         350 :   _preconditioner = preconditioner;
     127         350 : }
     128             : 
     129             : template <typename T>
     130             : void
     131      288378 : LinearSolver<T>::reuse_preconditioner(bool reuse_flag)
     132             : {
     133      288378 :   same_preconditioner = reuse_flag;
     134      288378 : }
     135             : 
     136             : template <typename T>
     137             : void
     138           0 : LinearSolver<T>::restrict_solve_to(const std::vector<unsigned int> * const dofs,
     139             :                                    const SubsetSolveMode /*subset_solve_mode*/)
     140             : {
     141           0 :   if (dofs != nullptr)
     142           0 :     libmesh_not_implemented();
     143           0 : }
     144             : 
     145             : 
     146             : template <typename T>
     147           0 : std::pair<unsigned int, Real> LinearSolver<T>::adjoint_solve (SparseMatrix<T> & mat,
     148             :                                                               NumericVector<T> & sol,
     149             :                                                               NumericVector<T> & rhs,
     150             :                                                               const std::optional<double> tol,
     151             :                                                               const std::optional<unsigned int> n_its)
     152             : {
     153             :   // Log how long the linear solve takes.
     154           0 :   LOG_SCOPE("adjoint_solve()", "LinearSolver");
     155             : 
     156             :   // Take the discrete adjoint
     157           0 :   mat.close();
     158           0 :   mat.get_transpose(mat);
     159             : 
     160             :   // Call the solve function for the relevant linear algebra library and
     161             :   // solve the transpose matrix
     162           0 :   const std::pair<unsigned int, Real> totalrval =  this->solve (mat, sol, rhs, tol, n_its);
     163             : 
     164             :   // Now transpose back and restore the original matrix
     165             :   // by taking the discrete adjoint
     166           0 :   mat.get_transpose(mat);
     167             : 
     168           0 :   return totalrval;
     169             : }
     170             : 
     171             : template <typename T>
     172           0 : void LinearSolver<T>::print_converged_reason() const
     173             : {
     174           0 :   LinearConvergenceReason reason = this->get_converged_reason();
     175           0 :   libMesh::out << "Linear solver convergence/divergence reason: " << Utility::enum_to_string(reason) << std::endl;
     176           0 : }
     177             : 
     178             : template <typename T>
     179      210203 : SolverConfiguration * LinearSolver<T>::solver_configuration() const
     180             : {
     181      210203 :   return _solver_configuration;
     182             : }
     183             : 
     184             : template <typename T>
     185         709 : void LinearSolver<T>::set_solver_configuration(SolverConfiguration * solver_configuration)
     186             : {
     187         709 :   _solver_configuration = solver_configuration;
     188         709 : }
     189             : 
     190             : template <typename T>
     191      900020 : double LinearSolver<T>::get_real_solver_setting (const std::string & setting_name,
     192             :                                                  const std::optional<double> & setting,
     193             :                                                  const std::optional<double> default_value)
     194             : {
     195      900020 :   if (setting.has_value())
     196      451618 :     return setting.value();
     197             : 
     198      448402 :   if (_solver_configuration)
     199         567 :     if (const auto it = this->_solver_configuration->real_valued_data.find(setting_name);
     200          32 :         it != this->_solver_configuration->real_valued_data.end())
     201         355 :       return double(it->second);
     202             : 
     203      448047 :   if (default_value.has_value())
     204      447976 :     return default_value.value();
     205             : 
     206         209 :   libmesh_error_msg("Linear solver setting '"
     207             :                     << setting_name
     208             :                     << "' must be supplied through an input argument, a SolverConfiguration "
     209             :                     << "object, or a default value.");
     210             : 
     211             :   return 0.0;
     212             : }
     213             : 
     214             : template <typename T>
     215      451973 : int LinearSolver<T>::get_int_solver_setting (const std::string & setting_name,
     216             :                                              const std::optional<int> & setting,
     217             :                                              const std::optional<int> default_value)
     218             : {
     219      451973 :   if (setting.has_value())
     220      451618 :     return setting.value();
     221             : 
     222         355 :   if (_solver_configuration)
     223         355 :     if (const auto it = this->_solver_configuration->int_valued_data.find(setting_name);
     224          20 :         it != this->_solver_configuration->int_valued_data.end())
     225         213 :       return it->second;
     226             : 
     227         142 :   if (default_value.has_value())
     228          71 :     return default_value.value();
     229             : 
     230         209 :   libmesh_error_msg("Linear solver setting '"
     231             :                     << setting_name
     232             :                     << "' must be supplied through an input argument, a SolverConfiguration "
     233             :                     << "object, or a default value.");
     234             : 
     235             :   return 0;
     236             : }
     237             : 
     238             : //------------------------------------------------------------------
     239             : // Explicit instantiations
     240             : template class LIBMESH_EXPORT LinearSolver<Number>;
     241             : 
     242             : 
     243             : 
     244             : } // namespace libMesh