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