Line data    Source code

       1             : // The libMesh Finite Element Library.
       2             : // Copyright (C) 2002-2025 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       23098 : LinearSolver<T>::LinearSolver (const libMesh::Parallel::Communicator & comm_in) :
      45             :   ParallelObject       (comm_in),
      46       21750 :   _solver_type         (GMRES),
      47       21750 :   _preconditioner_type (ILU_PRECOND),
      48       21750 :   _is_initialized      (false),
      49       21750 :   _preconditioner      (nullptr),
      50       21750 :   same_preconditioner  (false),
      51       23098 :   _solver_configuration(nullptr)
      52             : {
      53       23098 : }
      54             : 
      55             : 
      56             : 
      57             : template <typename T>
      58             : std::unique_ptr<LinearSolver<T>>
      59       23098 : 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         674 :   libmesh_ignore(comm);
      64             : 
      65             :   // Build the appropriate solver
      66       23098 :   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       22822 :     case PETSC_SOLVERS:
      76       22822 :       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         276 :     case EIGEN_SOLVERS:
      88         276 :       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          70 : LinearSolver<T>::attach_preconditioner(Preconditioner<T> * preconditioner)
     121             : {
     122          70 :   libmesh_error_msg_if(this->_is_initialized,
     123             :                        "Preconditioner must be attached before the solver is initialized!");
     124             : 
     125          70 :   _preconditioner_type = SHELL_PRECOND;
     126          70 :   _preconditioner = preconditioner;
     127          70 : }
     128             : 
     129             : template <typename T>
     130             : void
     131      288307 : LinearSolver<T>::reuse_preconditioner(bool reuse_flag)
     132             : {
     133      288307 :   same_preconditioner = reuse_flag;
     134      288307 : }
     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          70 : void LinearSolver<T>::set_solver_configuration(SolverConfiguration & solver_configuration)
     180             : {
     181          70 :   _solver_configuration = &solver_configuration;
     182          70 : }
     183             : 
     184             : template <typename T>
     185      859129 : double LinearSolver<T>::get_real_solver_setting (const std::string & setting_name,
     186             :                                                  const std::optional<double> & setting,
     187             :                                                  const std::optional<double> default_value)
     188             : {
     189      859129 :   if (setting.has_value())
     190      431348 :     return setting.value();
     191      427781 :   else if (_solver_configuration)
     192             :   {
     193          70 :     if (const auto it = this->_solver_configuration->real_valued_data.find(setting_name);
     194           4 :         it != this->_solver_configuration->real_valued_data.end())
     195           0 :       return double(it->second);
     196             :   }
     197      427711 :   else if (default_value.has_value())
     198      427711 :     return default_value.value();
     199             :   else
     200           0 :     libmesh_error_msg("Iteration configuration parameter to the linear solver should either be supplied through input arguments or a SolverConfiguration object!");
     201             : 
     202           2 :   return 0.0;
     203             : }
     204             : 
     205             : template <typename T>
     206      431348 : int LinearSolver<T>::get_int_solver_setting (const std::string & setting_name,
     207             :                                              const std::optional<int> & setting,
     208             :                                              const std::optional<int> default_value)
     209             : {
     210      431348 :   if (setting.has_value())
     211      431348 :     return setting.value();
     212           0 :   else if (_solver_configuration)
     213             :   {
     214           0 :     if (const auto it = this->_solver_configuration->int_valued_data.find(setting_name);
     215           0 :         it != this->_solver_configuration->int_valued_data.end())
     216           0 :       return it->second;
     217             :   }
     218           0 :   else if (default_value.has_value())
     219           0 :     return default_value.value();
     220             :   else
     221           0 :     libmesh_error_msg("Iteration configuration parameter to the linear solver should either be supplied through input arguments or a SolverConfiguration object!");
     222             : 
     223           0 :   return 0.0;
     224             : 
     225             : }
     226             : 
     227             : //------------------------------------------------------------------
     228             : // Explicit instantiations
     229             : template class LIBMESH_EXPORT LinearSolver<Number>;
     230             : 
     231             : 
     232             : 
     233             : } // namespace libMesh