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/nonlinear_implicit_system.h"
      22             : #include "libmesh/diff_solver.h"
      23             : #include "libmesh/equation_systems.h"
      24             : #include "libmesh/libmesh_logging.h"
      25             : #include "libmesh/nonlinear_solver.h"
      26             : #include "libmesh/sparse_matrix.h"
      27             : #include "libmesh/static_condensation.h"
      28             : #include "libmesh/static_condensation_preconditioner.h"
      29             : 
      30             : namespace libMesh
      31             : {
      32             : 
      33        1470 : NonlinearImplicitSystem::NonlinearImplicitSystem (EquationSystems & es,
      34             :                                                   const std::string & name_in,
      35        1470 :                                                   const unsigned int number_in) :
      36             : 
      37             :   Parent                    (es, name_in, number_in),
      38        1386 :   nonlinear_solver          (NonlinearSolver<Number>::build(*this)),
      39        1386 :   diff_solver               (),
      40        1386 :   _n_nonlinear_iterations   (0),
      41        1470 :   _final_nonlinear_residual (1.e20)
      42             : {
      43             :   // Set default parameters
      44             :   // These were chosen to match the Petsc defaults
      45        1470 :   es.parameters.set<Real>        ("linear solver tolerance") = 1e-5;
      46        1470 :   es.parameters.set<Real>        ("linear solver minimum tolerance") = 1e-5;
      47        1470 :   es.parameters.set<unsigned int>("linear solver maximum iterations") = 10000;
      48             : 
      49        1470 :   es.parameters.set<unsigned int>("nonlinear solver maximum iterations") = 50;
      50        1470 :   es.parameters.set<unsigned int>("nonlinear solver maximum function evaluations") = 10000;
      51             : 
      52        1470 :   es.parameters.set<Real>("nonlinear solver absolute residual tolerance") = 1e-35;
      53        1470 :   es.parameters.set<Real>("nonlinear solver relative residual tolerance") = 1e-8;
      54        1470 :   es.parameters.set<Real>("nonlinear solver divergence tolerance") = 1e+4;
      55        1470 :   es.parameters.set<Real>("nonlinear solver absolute step tolerance") = 1e-8;
      56        1470 :   es.parameters.set<Real>("nonlinear solver relative step tolerance") = 1e-8;
      57             : 
      58        1470 :   es.parameters.set<bool>("reuse preconditioner") = false;
      59        1470 :   es.parameters.set<unsigned int>("reuse preconditioner maximum linear iterations") = 1;
      60             : 
      61        1470 :   if (this->has_static_condensation())
      62         280 :     this->setup_static_condensation_preconditioner(*nonlinear_solver);
      63        1470 : }
      64             : 
      65             : 
      66             : 
      67        2508 : NonlinearImplicitSystem::~NonlinearImplicitSystem () = default;
      68             : 
      69             : 
      70             : 
      71           0 : void NonlinearImplicitSystem::create_static_condensation()
      72             : {
      73           0 :   Parent::create_static_condensation();
      74           0 :   this->setup_static_condensation_preconditioner(*nonlinear_solver);
      75           0 : }
      76             : 
      77             : 
      78             : 
      79         280 : void NonlinearImplicitSystem::clear ()
      80             : {
      81             :   // clear the nonlinear solver
      82         280 :   nonlinear_solver->clear();
      83             : 
      84             :   // FIXME - this is necessary for petsc_auto_fieldsplit
      85             :   // nonlinear_solver->init_names(*this);
      86             : 
      87             :   // clear the parent data
      88         280 :   Parent::clear();
      89             : 
      90             :   // And restore any StaticCondensation to defaults
      91         280 :   if (this->has_static_condensation())
      92           0 :     this->setup_static_condensation_preconditioner(*nonlinear_solver);
      93         280 : }
      94             : 
      95             : 
      96             : 
      97           0 : void NonlinearImplicitSystem::reinit ()
      98             : {
      99             :   // re-initialize the nonlinear solver interface
     100           0 :   nonlinear_solver->clear();
     101             : 
     102             :   // force the solver to get a new preconditioner, in
     103             :   // case reuse was set
     104           0 :   nonlinear_solver->force_new_preconditioner();
     105             : 
     106             :   // FIXME - this is necessary for petsc_auto_fieldsplit
     107             :   // nonlinear_solver->init_names(*this);
     108             : 
     109           0 :   if (diff_solver.get())
     110           0 :     diff_solver->reinit();
     111             : 
     112             :   // initialize parent data
     113           0 :   Parent::reinit();
     114           0 : }
     115             : 
     116             : 
     117             : 
     118       29400 : void NonlinearImplicitSystem::set_solver_parameters ()
     119             : {
     120             :   // Get a reference to the EquationSystems
     121             :   const EquationSystems & es =
     122        1680 :     this->get_equation_systems();
     123             : 
     124             :   // Get the user-specified nonlinear solver tolerances
     125       57960 :   const unsigned int maxits = parameters.have_parameter<unsigned int>("nonlinear solver maximum iterations") ?
     126           0 :     parameters.get<unsigned int>("nonlinear solver maximum iterations") :
     127       57960 :     es.parameters.get<unsigned int>("nonlinear solver maximum iterations");
     128             : 
     129       57960 :   const unsigned int maxfuncs = parameters.have_parameter<unsigned int>("nonlinear solver maximum function evaluations") ?
     130           0 :     parameters.get<unsigned int>("nonlinear solver maximum function evaluations") :
     131       57960 :     es.parameters.get<unsigned int>("nonlinear solver maximum function evaluations");
     132             : 
     133       57960 :   const double abs_resid_tol = parameters.have_parameter<Real>("nonlinear solver absolute residual tolerance") ?
     134           0 :     double(parameters.get<Real>("nonlinear solver absolute residual tolerance")) :
     135       57960 :     double(es.parameters.get<Real>("nonlinear solver absolute residual tolerance"));
     136             : 
     137       57960 :   const double rel_resid_tol = parameters.have_parameter<Real>("nonlinear solver relative residual tolerance") ?
     138           0 :     double(parameters.get<Real>("nonlinear solver relative residual tolerance")) :
     139       57960 :     double(es.parameters.get<Real>("nonlinear solver relative residual tolerance"));
     140             : 
     141       57960 :   const double div_tol = parameters.have_parameter<Real>("nonlinear solver divergence tolerance") ?
     142           0 :     double(parameters.get<Real>("nonlinear solver divergence tolerance")) :
     143       57960 :     double(es.parameters.get<Real>("nonlinear solver divergence tolerance"));
     144             : 
     145       57960 :   const double abs_step_tol = parameters.have_parameter<Real>("nonlinear solver absolute step tolerance") ?
     146           0 :     double(parameters.get<Real>("nonlinear solver absolute step tolerance")) :
     147       57960 :     double(es.parameters.get<Real>("nonlinear solver absolute step tolerance"));
     148             : 
     149       57960 :   const double rel_step_tol = parameters.have_parameter<Real>("nonlinear solver relative step tolerance")?
     150           0 :     double(parameters.get<Real>("nonlinear solver relative step tolerance")) :
     151       29400 :     double(es.parameters.get<Real>("nonlinear solver relative step tolerance"));
     152             : 
     153             :   // Get the user-specified linear solver tolerances
     154       29400 :   const auto [maxlinearits, linear_tol] = this->Parent::get_linear_solve_parameters();
     155             : 
     156       57960 :   const double linear_min_tol = parameters.have_parameter<Real>("linear solver minimum tolerance") ?
     157           0 :     double(parameters.get<Real>("linear solver minimum tolerance")) :
     158       57960 :     double(es.parameters.get<Real>("linear solver minimum tolerance"));
     159             : 
     160       29400 :   const bool reuse_preconditioner = parameters.have_parameter<unsigned int>("reuse preconditioner") ?
     161           0 :     parameters.get<unsigned int>("reuse preconditioner") :
     162       29400 :       es.parameters.get<bool>("reuse preconditioner");
     163             :   const unsigned int reuse_preconditioner_max_linear_its =
     164       57960 :     parameters.have_parameter<unsigned int>("reuse preconditioner maximum linear iterations") ?
     165           0 :     parameters.get<unsigned int>("reuse preconditioner maximum linear iterations") :
     166       57960 :       es.parameters.get<unsigned int>("reuse preconditioner maximum linear iterations");
     167             : 
     168             :   // Set all the parameters on the NonlinearSolver
     169       29400 :   nonlinear_solver->max_nonlinear_iterations = maxits;
     170       29400 :   nonlinear_solver->max_function_evaluations = maxfuncs;
     171       29400 :   nonlinear_solver->absolute_residual_tolerance = abs_resid_tol;
     172       29400 :   nonlinear_solver->relative_residual_tolerance = rel_resid_tol;
     173       29400 :   nonlinear_solver->divergence_tolerance = div_tol;
     174       29400 :   nonlinear_solver->absolute_step_tolerance = abs_step_tol;
     175       29400 :   nonlinear_solver->relative_step_tolerance = rel_step_tol;
     176       29400 :   nonlinear_solver->max_linear_iterations = maxlinearits;
     177       29400 :   nonlinear_solver->initial_linear_tolerance = linear_tol;
     178       29400 :   nonlinear_solver->minimum_linear_tolerance = linear_min_tol;
     179       29400 :   nonlinear_solver->set_reuse_preconditioner(reuse_preconditioner);
     180       29400 :   nonlinear_solver->set_reuse_preconditioner_max_linear_its(reuse_preconditioner_max_linear_its);
     181             : 
     182       29400 :   if (diff_solver.get())
     183             :     {
     184           0 :       diff_solver->max_nonlinear_iterations = maxits;
     185           0 :       diff_solver->absolute_residual_tolerance = abs_resid_tol;
     186           0 :       diff_solver->relative_residual_tolerance = rel_resid_tol;
     187           0 :       diff_solver->absolute_step_tolerance = abs_step_tol;
     188           0 :       diff_solver->relative_step_tolerance = rel_step_tol;
     189           0 :       diff_solver->max_linear_iterations = maxlinearits;
     190           0 :       diff_solver->initial_linear_tolerance = linear_tol;
     191           0 :       diff_solver->minimum_linear_tolerance = linear_min_tol;
     192             :     }
     193       29400 : }
     194             : 
     195             : 
     196             : 
     197       29400 : void NonlinearImplicitSystem::solve ()
     198             : {
     199             :   // Log how long the nonlinear solve takes.
     200        1680 :   LOG_SCOPE("solve()", "System");
     201             : 
     202       29400 :   this->set_solver_parameters();
     203             : 
     204       29400 :   if (diff_solver.get())
     205             :     {
     206           0 :       diff_solver->solve();
     207             : 
     208             :       // Store the number of nonlinear iterations required to
     209             :       // solve and the final residual.
     210           0 :       _n_nonlinear_iterations   = diff_solver->total_outer_iterations();
     211           0 :       _final_nonlinear_residual = 0.; // FIXME - support this!
     212             :     }
     213             :   else
     214             :     {
     215       29400 :       if (this->prefix_with_name())
     216           0 :         nonlinear_solver->init(this->prefix().c_str());
     217             :       else
     218       29400 :         nonlinear_solver->init();
     219             : 
     220             :       // FIXME - this is necessary for petsc_auto_fieldsplit
     221             :       // nonlinear_solver->init_names(*this);
     222             : 
     223             :       // Solve the nonlinear system.
     224             :       // Store the number of nonlinear iterations required to
     225             :       // solve and the final residual.
     226       29400 :       std::tie(_n_nonlinear_iterations, _final_nonlinear_residual) =
     227       30240 :         nonlinear_solver->solve (*matrix, *solution, *rhs,
     228         840 :                                  nonlinear_solver->relative_residual_tolerance,
     229        3360 :                                  nonlinear_solver->max_linear_iterations);
     230             :     }
     231             : 
     232             :   // Update the system after the solve
     233       29400 :   this->update();
     234       29400 : }
     235             : 
     236             : 
     237             : 
     238           0 : std::pair<unsigned int, Real> NonlinearImplicitSystem::get_linear_solve_parameters() const
     239             : {
     240           0 :   if (diff_solver.get())
     241           0 :     return std::make_pair(this->diff_solver->max_linear_iterations,
     242           0 :                           this->diff_solver->relative_residual_tolerance);
     243           0 :   return std::make_pair(this->nonlinear_solver->max_linear_iterations,
     244           0 :                         this->nonlinear_solver->relative_residual_tolerance);
     245             : }
     246             : 
     247             : 
     248             : 
     249           0 : void NonlinearImplicitSystem::assembly(bool get_residual,
     250             :                                        bool get_jacobian,
     251             :                                        bool /*apply_heterogeneous_constraints*/,
     252             :                                        bool /*apply_no_constraints*/)
     253             : {
     254             :   // Get current_local_solution in sync
     255           0 :   this->update();
     256             : 
     257             :   //-----------------------------------------------------------------------------
     258             :   // if the user has provided both function pointers and objects only the pointer
     259             :   // will be used, so catch that as an error
     260           0 :   libmesh_error_msg_if(nonlinear_solver->jacobian && nonlinear_solver->jacobian_object,
     261             :                        "ERROR: cannot specify both a function and object to compute the Jacobian!");
     262             : 
     263           0 :   libmesh_error_msg_if(nonlinear_solver->residual && nonlinear_solver->residual_object,
     264             :                        "ERROR: cannot specify both a function and object to compute the Residual!");
     265             : 
     266           0 :   libmesh_error_msg_if(nonlinear_solver->matvec && nonlinear_solver->residual_and_jacobian_object,
     267             :                        "ERROR: cannot specify both a function and object to compute the combined Residual & Jacobian!");
     268             : 
     269             : 
     270           0 :   if (get_jacobian)
     271             :     {
     272           0 :       if (nonlinear_solver->jacobian != nullptr)
     273           0 :         nonlinear_solver->jacobian (*current_local_solution.get(), *matrix, *this);
     274             : 
     275           0 :       else if (nonlinear_solver->jacobian_object != nullptr)
     276           0 :         nonlinear_solver->jacobian_object->jacobian (*current_local_solution.get(), *matrix, *this);
     277             : 
     278           0 :       else if (nonlinear_solver->matvec != nullptr)
     279           0 :         nonlinear_solver->matvec (*current_local_solution.get(), get_residual ? rhs : nullptr, matrix, *this);
     280             : 
     281           0 :       else if (nonlinear_solver->residual_and_jacobian_object != nullptr)
     282           0 :         nonlinear_solver->residual_and_jacobian_object->residual_and_jacobian (*current_local_solution.get(), get_residual ? rhs : nullptr, matrix, *this);
     283             : 
     284             :       else
     285           0 :         libmesh_error_msg("Error! Unable to compute residual and/or Jacobian!");
     286             :     }
     287             : 
     288           0 :   if (get_residual)
     289             :     {
     290           0 :       if (nonlinear_solver->residual != nullptr)
     291           0 :         nonlinear_solver->residual (*current_local_solution.get(), *rhs, *this);
     292             : 
     293           0 :       else if (nonlinear_solver->residual_object != nullptr)
     294           0 :         nonlinear_solver->residual_object->residual (*current_local_solution.get(), *rhs, *this);
     295             : 
     296           0 :       else if (nonlinear_solver->matvec != nullptr)
     297             :         {
     298             :           // we might have already grabbed the residual and jacobian together
     299           0 :           if (!get_jacobian)
     300           0 :             nonlinear_solver->matvec (*current_local_solution.get(), rhs, nullptr, *this);
     301             :         }
     302             : 
     303           0 :       else if (nonlinear_solver->residual_and_jacobian_object != nullptr)
     304             :         {
     305             :           // we might have already grabbed the residual and jacobian together
     306           0 :           if (!get_jacobian)
     307           0 :             nonlinear_solver->residual_and_jacobian_object->residual_and_jacobian (*current_local_solution.get(), rhs, nullptr, *this);
     308             :         }
     309             : 
     310             :       else
     311           0 :         libmesh_error_msg("Error! Unable to compute residual and/or Jacobian!");
     312             :     }
     313             :   else
     314           0 :     libmesh_assert(get_jacobian);  // I can't believe you really wanted to assemble *nothing*
     315           0 : }
     316             : 
     317             : 
     318             : 
     319             : 
     320           0 : unsigned NonlinearImplicitSystem::get_current_nonlinear_iteration_number() const
     321             : {
     322           0 :   return nonlinear_solver->get_current_nonlinear_iteration_number();
     323             : }
     324             : 
     325             : 
     326             : 
     327             : } // namespace libMesh