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             : #ifndef LIBMESH_DIFF_SOLVER_H
      21             : #define LIBMESH_DIFF_SOLVER_H
      22             : 
      23             : // Local includes
      24             : #include "libmesh/libmesh_common.h"
      25             : #include "libmesh/reference_counted_object.h"
      26             : #include "libmesh/parallel_object.h"
      27             : 
      28             : // C++ includes
      29             : #include <vector>
      30             : #include <memory>
      31             : 
      32             : namespace libMesh
      33             : {
      34             : 
      35             : // Forward Declarations
      36             : class ImplicitSystem;
      37             : template <typename T> class NumericVector;
      38             : 
      39             : /**
      40             :  * Functor for use as callback in solve of nonlinear solver.
      41             :  */
      42             : class LinearSolutionMonitor
      43             : {
      44             : public:
      45             :   virtual void operator() (const NumericVector<Number> & delta_u,
      46             :                            const Real & norm_delta_u,
      47             :                            const NumericVector<Number> & u,
      48             :                            const Real & norm_u,
      49             :                            const NumericVector<Number> & res,
      50             :                            const Real & norm_res,
      51             :                            const unsigned int iteration) = 0;
      52             :   virtual ~LinearSolutionMonitor() = default;
      53             : };
      54             : 
      55             : /**
      56             :  * This is a generic class that defines a solver to handle
      57             :  * ImplicitSystem classes, including NonlinearImplicitSystem and
      58             :  * DifferentiableSystem   A user can define a solver by
      59             :  * deriving from this class and implementing certain functions.
      60             :  *
      61             :  * This class is part of the new DifferentiableSystem framework,
      62             :  * which is still experimental.  Users of this framework should
      63             :  * beware of bugs and future API changes.
      64             :  *
      65             :  * \author Roy H. Stogner
      66             :  * \date 2006-2010
      67             :  */
      68             : class DiffSolver : public ReferenceCountedObject<DiffSolver>,
      69             :                    public ParallelObject
      70             : {
      71             : public:
      72             :   /**
      73             :    * The type of system
      74             :    */
      75             :   typedef ImplicitSystem sys_type;
      76             : 
      77             :   /**
      78             :    * Constructor. Requires a reference to the system
      79             :    * to be solved.
      80             :    */
      81             :   DiffSolver (sys_type & s);
      82             : 
      83             :   /**
      84             :    * Factory method.  Requires a reference to the system
      85             :    * to be solved.
      86             :    *
      87             :    * \returns A NewtonSolver by default.
      88             :    */
      89             :   static std::unique_ptr<DiffSolver> build(sys_type & s);
      90             : 
      91             :   /**
      92             :    * Destructor.
      93             :    */
      94        5395 :   virtual ~DiffSolver () = default;
      95             : 
      96             :   /**
      97             :    * The initialization function.  This method is used to
      98             :    * initialize internal data structures before a simulation begins.
      99             :    */
     100             :   virtual void init ();
     101             : 
     102             :   /**
     103             :    * The reinitialization function.  This method is used after
     104             :    * changes in the mesh.
     105             :    */
     106             :   virtual void reinit ();
     107             : 
     108             :   /**
     109             :    * This method performs a solve.  What occurs in
     110             :    * this method will depend on the type of solver.  See
     111             :    * the subclasses for more details.
     112             :    */
     113             :   virtual unsigned int solve () = 0;
     114             : 
     115             :   /**
     116             :    * \returns The number of "outer" (e.g. quasi-Newton) iterations
     117             :    * required by the last solve.
     118             :    */
     119           0 :   unsigned int total_outer_iterations() { return _outer_iterations; }
     120             : 
     121             :   /**
     122             :    * \returns The number of "inner" (e.g. Krylov) iterations
     123             :    * required by the last solve.
     124             :    */
     125             :   unsigned int total_inner_iterations() { return _inner_iterations; }
     126             : 
     127             :   /**
     128             :    * \returns The value of the SolveResult from the last solve.
     129             :    */
     130             :   unsigned int solve_result() { return _solve_result; }
     131             : 
     132             :   /**
     133             :    * \returns A constant reference to the system we are solving.
     134             :    */
     135             :   const sys_type & system () const { return _system; }
     136             : 
     137             :   /**
     138             :    * \returns A writable reference to the system we are solving.
     139             :    */
     140       19069 :   sys_type & system () { return _system; }
     141             : 
     142             :   /**
     143             :    * Each linear solver step should exit after \p max_linear_iterations
     144             :    * is exceeded.
     145             :    */
     146             :   unsigned int max_linear_iterations;
     147             : 
     148             :   /**
     149             :    * The DiffSolver should exit in failure if \p max_nonlinear_iterations
     150             :    * is exceeded and \p continue_after_max_iterations is false, or should
     151             :    * end the nonlinear solve if \p max_nonlinear_iterations is exceeded and \p
     152             :    * continue_after_max_iterations is true.
     153             :    */
     154             :   unsigned int max_nonlinear_iterations;
     155             : 
     156             :   /**
     157             :    * The DiffSolver should not print anything to libMesh::out
     158             :    * unless quiet is set to false; default is true.
     159             :    */
     160             :   bool quiet;
     161             : 
     162             :   /**
     163             :    * The DiffSolver may print a lot more to libMesh::out
     164             :    * if verbose is set to true; default is false.
     165             :    */
     166             :   bool verbose;
     167             : 
     168             :   /**
     169             :    * Defaults to true, telling the DiffSolver to continue rather than exit when
     170             :    * a solve has reached its maximum number of nonlinear iterations.
     171             :    */
     172             :   bool continue_after_max_iterations;
     173             : 
     174             :   /**
     175             :    * Defaults to false, telling the DiffSolver to throw an error when
     176             :    * the backtracking scheme fails to find a descent direction.
     177             :    */
     178             :   bool continue_after_backtrack_failure;
     179             : 
     180             :   /**
     181             :    * The DiffSolver should exit after the residual is
     182             :    * reduced to either less than absolute_residual_tolerance
     183             :    * or less than relative_residual_tolerance times the
     184             :    * initial residual.
     185             :    *
     186             :    * Users should increase any of these tolerances that they want to use for a
     187             :    * stopping condition.
     188             :    */
     189             :   Real absolute_residual_tolerance;
     190             :   Real relative_residual_tolerance;
     191             : 
     192             :   /**
     193             :    * The DiffSolver should exit after the full nonlinear step norm is
     194             :    * reduced to either less than absolute_step_tolerance
     195             :    * or less than relative_step_tolerance times the largest
     196             :    * nonlinear solution which has been seen so far.
     197             :    *
     198             :    * Users should increase any of these tolerances that they want to use for a
     199             :    * stopping condition.
     200             :    */
     201             :   Real absolute_step_tolerance;
     202             :   Real relative_step_tolerance;
     203             : 
     204             :   /**
     205             :    * Any required linear solves will at first be done with this tolerance;
     206             :    * the DiffSolver may tighten the tolerance for later solves.
     207             :    */
     208             :   double initial_linear_tolerance;
     209             : 
     210             :   /**
     211             :    * The tolerance for linear solves is kept above this minimum
     212             :    */
     213             :   double minimum_linear_tolerance;
     214             : 
     215             :   /**
     216             :    * Enumeration return type for the solve() function.  Multiple SolveResults
     217             :    * may be combined (OR'd) in the single return.  To test which ones are present,
     218             :    * just AND the return value with any of the SolveResult flags defined below.
     219             :    */
     220             :   enum SolveResult {
     221             :     /**
     222             :      * A default or invalid solve result.  This usually means
     223             :      * no solve has occurred yet.
     224             :      */
     225             :     INVALID_SOLVE_RESULT = 0,
     226             : 
     227             :     /**
     228             :      * The solver converged but no
     229             :      * particular reason is specified.
     230             :      */
     231             :     CONVERGED_NO_REASON = 1,
     232             : 
     233             :     /**
     234             :      * The DiffSolver achieved the desired
     235             :      * absolute residual tolerance.
     236             :      */
     237             :     CONVERGED_ABSOLUTE_RESIDUAL = 2,
     238             : 
     239             :     /**
     240             :      * The DiffSolver achieved the desired
     241             :      * relative residual tolerance.
     242             :      */
     243             :     CONVERGED_RELATIVE_RESIDUAL = 4,
     244             : 
     245             :     /**
     246             :      * The DiffSolver achieved the desired
     247             :      * absolute step size tolerance.
     248             :      */
     249             :     CONVERGED_ABSOLUTE_STEP = 8,
     250             : 
     251             :     /**
     252             :      * The DiffSolver achieved the desired
     253             :      * relative step size tolerance.
     254             :      */
     255             :     CONVERGED_RELATIVE_STEP = 16,
     256             : 
     257             :     /**
     258             :      * The DiffSolver diverged but no
     259             :      * particular reason is specified.
     260             :      */
     261             :     DIVERGED_NO_REASON = 32,
     262             : 
     263             :     /**
     264             :      * The DiffSolver reached the maximum allowed
     265             :      * number of nonlinear iterations before satisfying
     266             :      * any convergence tests.
     267             :      */
     268             :     DIVERGED_MAX_NONLINEAR_ITERATIONS = 64,
     269             : 
     270             :     /**
     271             :      * The DiffSolver failed to find a descent direction
     272             :      * by backtracking (See newton_solver.C)
     273             :      */
     274             :     DIVERGED_BACKTRACKING_FAILURE = 128,
     275             : 
     276             :     /**
     277             :      * The linear solver used by the DiffSolver failed to
     278             :      * find a solution.
     279             :      */
     280             :     DIVERGED_LINEAR_SOLVER_FAILURE = 256
     281             :   };
     282             : 
     283             :   /**
     284             :    * Pointer to functor which is called right after each linear solve
     285             :    */
     286             :   std::unique_ptr<LinearSolutionMonitor> linear_solution_monitor;
     287             : 
     288             :   /**
     289             :    * Enable (or disable; it is \p true by default) exact enforcement
     290             :    * of constraints at the solver level, correcting any constrained
     291             :    * DoF coefficients in \p current_local_solution as well as applying
     292             :    * nonlinear residual and Jacobian terms based on constraint
     293             :    * equations.
     294             :    *
     295             :    * This is probably only safe to disable if user code is setting
     296             :    * nonlinear residual and Jacobian terms based on constraint
     297             :    * equations at an element-by-element level, by combining the
     298             :    * \p asymmetric_constraint_rows option with the
     299             :    * \p residual_constrain_element_vector processing option in
     300             :    * \p DofMap.
     301             :    */
     302        1573 :   virtual void set_exact_constraint_enforcement(bool enable)
     303             :   {
     304        1573 :     _exact_constraint_enforcement = enable;
     305        1573 :   }
     306             : 
     307        1080 :   bool exact_constraint_enforcement()
     308             :   {
     309       37170 :     return _exact_constraint_enforcement;
     310             :   }
     311             : 
     312             : protected:
     313             : 
     314             :   /**
     315             :    * Whether we should enforce exact constraints globally during a
     316             :    * solve.
     317             :    */
     318             :   bool _exact_constraint_enforcement;
     319             : 
     320             :   /**
     321             :    * The largest solution norm which the DiffSolver has yet seen will be stored
     322             :    * here, to be used for stopping criteria based on relative_step_tolerance
     323             :    */
     324             :   Real max_solution_norm;
     325             : 
     326             :   /**
     327             :    * The largest nonlinear residual which the DiffSolver has yet seen will be
     328             :    * stored here, to be used for stopping criteria based on
     329             :    * relative_residual_tolerance
     330             :    */
     331             :   Real max_residual_norm;
     332             : 
     333             :   /**
     334             :    * The number of outer iterations used by the last solve
     335             :    */
     336             :   unsigned int _outer_iterations;
     337             : 
     338             :   /**
     339             :    * The number of inner iterations used by the last solve
     340             :    */
     341             :   unsigned int _inner_iterations;
     342             : 
     343             :   /**
     344             :    * A reference to the system we are solving.
     345             :    */
     346             :   sys_type & _system;
     347             : 
     348             :   /**
     349             :    * Initialized to zero.  solve_result is typically set internally in
     350             :    * the solve() function before it returns.  When non-zero,
     351             :    * solve_result tells the result of the latest solve.  See enum
     352             :    * definition for description.
     353             :    */
     354             :   unsigned int _solve_result;
     355             : };
     356             : 
     357             : 
     358             : } // namespace libMesh
     359             : 
     360             : 
     361             : #endif // LIBMESH_DIFF_SOLVER_H