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