Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://mooseframework.inl.gov
       3             : //*
       4             : //* All rights reserved, see COPYRIGHT for full restrictions
       5             : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
       6             : //*
       7             : //* Licensed under LGPL 2.1, please see LICENSE for details
       8             : //* https://www.gnu.org/licenses/lgpl-2.1.html
       9             : 
      10             : // moose includes
      11             : #include "NonlinearSystem.h"
      12             : #include "FEProblem.h"
      13             : #include "DisplacedProblem.h"
      14             : #include "TimeIntegrator.h"
      15             : #include "FiniteDifferencePreconditioner.h"
      16             : #include "PetscSupport.h"
      17             : #include "ComputeResidualFunctor.h"
      18             : #include "ComputeFDResidualFunctor.h"
      19             : #include "MooseVariableScalar.h"
      20             : #include "MooseTypes.h"
      21             : #include "AuxiliarySystem.h"
      22             : #include "Console.h"
      23             : 
      24             : #include "libmesh/nonlinear_solver.h"
      25             : #include "libmesh/petsc_nonlinear_solver.h"
      26             : #include "libmesh/sparse_matrix.h"
      27             : #include "libmesh/petsc_matrix.h"
      28             : #include "libmesh/diagonal_matrix.h"
      29             : #include "libmesh/default_coupling.h"
      30             : #include "libmesh/petsc_solver_exception.h"
      31             : 
      32             : using namespace libMesh;
      33             : 
      34             : namespace Moose
      35             : {
      36             : void
      37      470651 : compute_jacobian(const NumericVector<Number> & soln,
      38             :                  SparseMatrix<Number> & jacobian,
      39             :                  NonlinearImplicitSystem & sys)
      40             : {
      41             :   FEProblemBase * p =
      42      470651 :       sys.get_equation_systems().parameters.get<FEProblemBase *>("_fe_problem_base");
      43      470651 :   p->computeJacobianSys(sys, soln, jacobian);
      44      470638 : }
      45             : 
      46             : void
      47         714 : compute_bounds(NumericVector<Number> & lower,
      48             :                NumericVector<Number> & upper,
      49             :                NonlinearImplicitSystem & sys)
      50             : {
      51             :   FEProblemBase * p =
      52         714 :       sys.get_equation_systems().parameters.get<FEProblemBase *>("_fe_problem_base");
      53         714 :   p->computeBounds(sys, lower, upper);
      54         714 : }
      55             : 
      56             : void
      57      294844 : compute_nullspace(std::vector<NumericVector<Number> *> & sp, NonlinearImplicitSystem & sys)
      58             : {
      59             :   FEProblemBase * p =
      60      294844 :       sys.get_equation_systems().parameters.get<FEProblemBase *>("_fe_problem_base");
      61      294844 :   p->computeNullSpace(sys, sp);
      62      294844 : }
      63             : 
      64             : void
      65      294844 : compute_transpose_nullspace(std::vector<NumericVector<Number> *> & sp,
      66             :                             NonlinearImplicitSystem & sys)
      67             : {
      68             :   FEProblemBase * p =
      69      294844 :       sys.get_equation_systems().parameters.get<FEProblemBase *>("_fe_problem_base");
      70      294844 :   p->computeTransposeNullSpace(sys, sp);
      71      294844 : }
      72             : 
      73             : void
      74      294844 : compute_nearnullspace(std::vector<NumericVector<Number> *> & sp, NonlinearImplicitSystem & sys)
      75             : {
      76             :   FEProblemBase * p =
      77      294844 :       sys.get_equation_systems().parameters.get<FEProblemBase *>("_fe_problem_base");
      78      294844 :   p->computeNearNullSpace(sys, sp);
      79      294844 : }
      80             : 
      81             : void
      82        2128 : compute_postcheck(const NumericVector<Number> & old_soln,
      83             :                   NumericVector<Number> & search_direction,
      84             :                   NumericVector<Number> & new_soln,
      85             :                   bool & changed_search_direction,
      86             :                   bool & changed_new_soln,
      87             :                   NonlinearImplicitSystem & sys)
      88             : {
      89             :   FEProblemBase * p =
      90        2128 :       sys.get_equation_systems().parameters.get<FEProblemBase *>("_fe_problem_base");
      91        2128 :   p->computePostCheck(
      92             :       sys, old_soln, search_direction, new_soln, changed_search_direction, changed_new_soln);
      93        2128 : }
      94             : } // namespace Moose
      95             : 
      96       60299 : NonlinearSystem::NonlinearSystem(FEProblemBase & fe_problem, const std::string & name)
      97             :   : NonlinearSystemBase(
      98       60299 :         fe_problem, fe_problem.es().add_system<NonlinearImplicitSystem>(name), name),
      99       60299 :     _nl_implicit_sys(fe_problem.es().get_system<NonlinearImplicitSystem>(name)),
     100       60299 :     _nl_residual_functor(_fe_problem),
     101       60299 :     _fd_residual_functor(_fe_problem),
     102       60299 :     _resid_and_jac_functor(_fe_problem),
     103      120598 :     _use_coloring_finite_difference(false)
     104             : {
     105       60299 :   nonlinearSolver()->residual_object = &_nl_residual_functor;
     106       60299 :   nonlinearSolver()->jacobian = Moose::compute_jacobian;
     107       60299 :   nonlinearSolver()->bounds = Moose::compute_bounds;
     108       60299 :   nonlinearSolver()->nullspace = Moose::compute_nullspace;
     109       60299 :   nonlinearSolver()->transpose_nullspace = Moose::compute_transpose_nullspace;
     110       60299 :   nonlinearSolver()->nearnullspace = Moose::compute_nearnullspace;
     111             : 
     112             :   PetscNonlinearSolver<Real> * petsc_solver =
     113       60299 :       static_cast<PetscNonlinearSolver<Real> *>(_nl_implicit_sys.nonlinear_solver.get());
     114       60299 :   if (petsc_solver)
     115             :   {
     116       60299 :     petsc_solver->set_residual_zero_out(false);
     117       60299 :     petsc_solver->set_jacobian_zero_out(false);
     118       60299 :     petsc_solver->use_default_monitor(false);
     119             :   }
     120       60299 : }
     121             : 
     122       57212 : NonlinearSystem::~NonlinearSystem() {}
     123             : 
     124             : void
     125      292764 : NonlinearSystem::potentiallySetupFiniteDifferencing()
     126             : {
     127      292764 :   if (_use_finite_differenced_preconditioner)
     128             :   {
     129         182 :     _nl_implicit_sys.nonlinear_solver->fd_residual_object = &_fd_residual_functor;
     130         182 :     setupFiniteDifferencedPreconditioner();
     131             :   }
     132             : 
     133             :   PetscNonlinearSolver<Real> & solver =
     134      292764 :       static_cast<PetscNonlinearSolver<Real> &>(*_nl_implicit_sys.nonlinear_solver);
     135      292764 :   solver.mffd_residual_object = &_fd_residual_functor;
     136             : 
     137      292764 :   solver.set_snesmf_reuse_base(_fe_problem.useSNESMFReuseBase());
     138      292764 : }
     139             : 
     140             : void
     141      287715 : NonlinearSystem::solve()
     142             : {
     143             :   // Only attach the postcheck function to the solver if we actually
     144             :   // have dampers or if the FEProblemBase needs to update the solution,
     145             :   // which is also done during the linesearch postcheck.  It doesn't
     146             :   // hurt to do this multiple times, it is just setting a pointer.
     147      575027 :   if (_fe_problem.hasDampers() || _fe_problem.shouldUpdateSolution() ||
     148      287312 :       _fe_problem.needsPreviousNewtonIteration())
     149         502 :     _nl_implicit_sys.nonlinear_solver->postcheck = Moose::compute_postcheck;
     150             : 
     151      287715 :   if (shouldEvaluatePreSMOResidual())
     152             :   {
     153         415 :     TIME_SECTION("nlPreSMOResidual", 3, "Computing Pre-SMO Residual");
     154             :     // Calculate the pre-SMO residual for use in the convergence criterion.
     155          83 :     _computing_pre_smo_residual = true;
     156          83 :     _fe_problem.computeResidualSys(_nl_implicit_sys, *_current_solution, *_nl_implicit_sys.rhs);
     157          83 :     _computing_pre_smo_residual = false;
     158          83 :     _nl_implicit_sys.rhs->close();
     159          83 :     _pre_smo_residual = _nl_implicit_sys.rhs->l2_norm();
     160          83 :     _console << " * Nonlinear |R| = "
     161         166 :              << Console::outputNorm(std::numeric_limits<Real>::max(), _pre_smo_residual)
     162          83 :              << " (Before preset BCs, predictors, correctors, and constraints)\n";
     163          83 :     _console << std::flush;
     164          83 :   }
     165             : 
     166      287715 :   const bool presolve_succeeded = preSolve();
     167      287715 :   if (!presolve_succeeded)
     168           0 :     return;
     169             : 
     170      287715 :   potentiallySetupFiniteDifferencing();
     171             : 
     172      287715 :   const bool time_integrator_solve = std::any_of(_time_integrators.begin(),
     173             :                                                  _time_integrators.end(),
     174      257955 :                                                  [](auto & ti) { return ti->overridesSolve(); });
     175             :   if (time_integrator_solve)
     176             :     mooseAssert(_time_integrators.size() == 1,
     177             :                 "If solve is overridden, then there must be only one time integrator");
     178             : 
     179      287715 :   if (time_integrator_solve)
     180       11099 :     _time_integrators.front()->solve();
     181             :   else
     182      276616 :     system().solve();
     183             : 
     184      545604 :   for (auto & ti : _time_integrators)
     185             :   {
     186      257941 :     if (!ti->overridesSolve())
     187      246842 :       ti->setNumIterationsLastSolve();
     188      257941 :     ti->postSolve();
     189             :   }
     190             : 
     191      287663 :   if (!_time_integrators.empty())
     192             :   {
     193      257601 :     _n_iters = _time_integrators.front()->getNumNonlinearIterations();
     194      257601 :     _n_linear_iters = _time_integrators.front()->getNumLinearIterations();
     195             :   }
     196             :   else
     197             :   {
     198       30062 :     _n_iters = _nl_implicit_sys.n_nonlinear_iterations();
     199       30062 :     _n_linear_iters = _nl_implicit_sys.nonlinear_solver->get_total_linear_iterations();
     200             :   }
     201             : 
     202             :   // store info about the solve
     203      287663 :   _final_residual = _nl_implicit_sys.final_nonlinear_residual();
     204             : 
     205             :   // determine whether solution invalid occurs in the converged solution
     206      287663 :   checkInvalidSolution();
     207             : 
     208      287660 :   if (_use_coloring_finite_difference)
     209          63 :     LibmeshPetscCall(MatFDColoringDestroy(&_fdcoloring));
     210             : }
     211             : 
     212             : void
     213         304 : NonlinearSystem::stopSolve(const ExecFlagType & exec_flag,
     214             :                            const std::set<TagID> & vector_tags_to_close)
     215             : {
     216             :   PetscNonlinearSolver<Real> & solver =
     217         304 :       static_cast<PetscNonlinearSolver<Real> &>(*sys().nonlinear_solver);
     218             : 
     219         304 :   if (exec_flag == EXEC_LINEAR || exec_flag == EXEC_POSTCHECK)
     220             :   {
     221         241 :     LibmeshPetscCall(SNESSetFunctionDomainError(solver.snes()));
     222             : 
     223             :     // Clean up by getting vectors into a valid state for a
     224             :     // (possible) subsequent solve.
     225         241 :     closeTaggedVectors(vector_tags_to_close);
     226             :   }
     227          63 :   else if (exec_flag == EXEC_NONLINEAR)
     228          63 :     LibmeshPetscCall(SNESSetJacobianDomainError(solver.snes()));
     229             :   else
     230           0 :     mooseError("Unsupported execute flag: ", Moose::stringify(exec_flag));
     231         304 : }
     232             : 
     233             : void
     234         182 : NonlinearSystem::setupFiniteDifferencedPreconditioner()
     235             : {
     236             :   std::shared_ptr<FiniteDifferencePreconditioner> fdp =
     237         182 :       std::dynamic_pointer_cast<FiniteDifferencePreconditioner>(_preconditioner);
     238         182 :   if (!fdp)
     239           0 :     mooseError("Did not setup finite difference preconditioner, and please add a preconditioning "
     240             :                "block with type = fdp");
     241             : 
     242         182 :   if (fdp->finiteDifferenceType() == "coloring")
     243             :   {
     244         171 :     _use_coloring_finite_difference = true;
     245         171 :     setupColoringFiniteDifferencedPreconditioner();
     246             :   }
     247             : 
     248          11 :   else if (fdp->finiteDifferenceType() == "standard")
     249             :   {
     250          11 :     _use_coloring_finite_difference = false;
     251          11 :     setupStandardFiniteDifferencedPreconditioner();
     252             :   }
     253             :   else
     254           0 :     mooseError("Unknown finite difference type");
     255         182 : }
     256             : 
     257             : void
     258          11 : NonlinearSystem::setupStandardFiniteDifferencedPreconditioner()
     259             : {
     260             :   // Make sure that libMesh isn't going to override our preconditioner
     261          11 :   _nl_implicit_sys.nonlinear_solver->jacobian = nullptr;
     262             : 
     263             :   PetscNonlinearSolver<Number> * petsc_nonlinear_solver =
     264          11 :       static_cast<PetscNonlinearSolver<Number> *>(_nl_implicit_sys.nonlinear_solver.get());
     265             : 
     266             :   PetscMatrix<Number> * petsc_mat =
     267          11 :       static_cast<PetscMatrix<Number> *>(&_nl_implicit_sys.get_system_matrix());
     268             : 
     269          11 :   LibmeshPetscCall(SNESSetJacobian(petsc_nonlinear_solver->snes(),
     270             :                                    petsc_mat->mat(),
     271             :                                    petsc_mat->mat(),
     272             :                                    SNESComputeJacobianDefault,
     273             :                                    nullptr));
     274          11 : }
     275             : 
     276             : void
     277         171 : NonlinearSystem::setupColoringFiniteDifferencedPreconditioner()
     278             : {
     279             :   // Make sure that libMesh isn't going to override our preconditioner
     280         171 :   _nl_implicit_sys.nonlinear_solver->jacobian = nullptr;
     281             : 
     282             :   PetscNonlinearSolver<Number> & petsc_nonlinear_solver =
     283         171 :       dynamic_cast<PetscNonlinearSolver<Number> &>(*_nl_implicit_sys.nonlinear_solver);
     284             : 
     285             :   // Pointer to underlying PetscMatrix type
     286             :   PetscMatrix<Number> * petsc_mat =
     287         171 :       dynamic_cast<PetscMatrix<Number> *>(&_nl_implicit_sys.get_system_matrix());
     288             : 
     289         171 :   Moose::compute_jacobian(*_nl_implicit_sys.current_local_solution, *petsc_mat, _nl_implicit_sys);
     290             : 
     291         171 :   if (!petsc_mat)
     292           0 :     mooseError("Could not convert to Petsc matrix.");
     293             : 
     294         171 :   petsc_mat->close();
     295             : 
     296             :   ISColoring iscoloring;
     297             : 
     298             :   // PETSc 3.5.x
     299             :   MatColoring matcoloring;
     300         171 :   LibmeshPetscCallA(_communicator.get(), MatColoringCreate(petsc_mat->mat(), &matcoloring));
     301         171 :   LibmeshPetscCallA(_communicator.get(), MatColoringSetType(matcoloring, MATCOLORINGLF));
     302         171 :   LibmeshPetscCallA(_communicator.get(), MatColoringSetFromOptions(matcoloring));
     303         171 :   LibmeshPetscCallA(_communicator.get(), MatColoringApply(matcoloring, &iscoloring));
     304         171 :   LibmeshPetscCallA(_communicator.get(), MatColoringDestroy(&matcoloring));
     305             : 
     306         171 :   LibmeshPetscCallA(_communicator.get(),
     307             :                     MatFDColoringCreate(petsc_mat->mat(), iscoloring, &_fdcoloring));
     308         171 :   LibmeshPetscCallA(_communicator.get(), MatFDColoringSetFromOptions(_fdcoloring));
     309             :   // clang-format off
     310             : #if PETSC_VERSION_LESS_THAN(3, 24, 0)
     311             :   LibmeshPetscCallA(_communicator.get(),
     312             :                     MatFDColoringSetFunction(_fdcoloring,
     313             :                                              (PetscErrorCode(*)(void))(void (*)(void))
     314             :                                              &libMesh::libmesh_petsc_snes_fd_residual,
     315             :                                              &petsc_nonlinear_solver));
     316             : #else
     317         171 :   LibmeshPetscCallA(_communicator.get(),
     318             :                     MatFDColoringSetFunction(_fdcoloring,
     319             :                                              (MatFDColoringFn*)
     320             :                                              &libMesh::libmesh_petsc_snes_fd_residual,
     321             :                                              &petsc_nonlinear_solver));
     322             : #endif
     323             :   // clang-format on
     324         171 :   LibmeshPetscCallA(_communicator.get(),
     325             :                     MatFDColoringSetUp(petsc_mat->mat(), iscoloring, _fdcoloring));
     326         171 :   LibmeshPetscCallA(_communicator.get(),
     327             :                     SNESSetJacobian(petsc_nonlinear_solver.snes(),
     328             :                                     petsc_mat->mat(),
     329             :                                     petsc_mat->mat(),
     330             :                                     SNESComputeJacobianDefaultColor,
     331             :                                     _fdcoloring));
     332             :   // PETSc >=3.3.0
     333         171 :   LibmeshPetscCallA(_communicator.get(), ISColoringDestroy(&iscoloring));
     334         171 : }
     335             : 
     336             : bool
     337      308801 : NonlinearSystem::converged()
     338             : {
     339      308801 :   if (_fe_problem.hasException() || _fe_problem.getFailNextNonlinearConvergenceCheck())
     340         172 :     return false;
     341             :   // When not computing the residual (for example at the beginning of a time step),
     342             :   // we may be in the process of counting invalid solution warnings, so the call to
     343             :   // acceptInvalidSolution() would fail due to lack of parallel synchronization
     344             :   // TODO: think of a better solution
     345      308629 :   if (_app.solutionInvalidity().hasSynced() && !_fe_problem.acceptInvalidSolution())
     346             :   {
     347         462 :     mooseWarning("The solution is not converged due to the solution being invalid.");
     348         459 :     return false;
     349             :   }
     350      308167 :   return _nl_implicit_sys.nonlinear_solver->converged;
     351             : }
     352             : 
     353             : void
     354         107 : NonlinearSystem::attachPreconditioner(Preconditioner<Number> * preconditioner)
     355             : {
     356         107 :   nonlinearSolver()->attach_preconditioner(preconditioner);
     357         107 : }
     358             : 
     359             : void
     360         474 : NonlinearSystem::computeScalingJacobian()
     361             : {
     362         474 :   _fe_problem.computeJacobianSys(_nl_implicit_sys, *_current_solution, *_scaling_matrix);
     363         474 : }
     364             : 
     365             : void
     366          36 : NonlinearSystem::computeScalingResidual()
     367             : {
     368          36 :   _fe_problem.computeResidualSys(_nl_implicit_sys, *_current_solution, RHS());
     369          36 : }
     370             : 
     371             : SNES
     372     1084358 : NonlinearSystem::getSNES()
     373             : {
     374             :   PetscNonlinearSolver<Number> * petsc_solver =
     375     1084358 :       dynamic_cast<PetscNonlinearSolver<Number> *>(nonlinearSolver());
     376             : 
     377     1084358 :   if (petsc_solver)
     378             :   {
     379     1084358 :     const char * snes_prefix = nullptr;
     380     1084358 :     std::string snes_prefix_str;
     381     1084358 :     if (system().prefix_with_name())
     382             :     {
     383       39244 :       snes_prefix_str = system().prefix();
     384       39244 :       snes_prefix = snes_prefix_str.c_str();
     385             :     }
     386     2168716 :     return petsc_solver->snes(snes_prefix);
     387     1084358 :   }
     388             :   else
     389           0 :     mooseError("It is not a petsc nonlinear solver");
     390             : }
     391             : 
     392             : void
     393         460 : NonlinearSystem::residualAndJacobianTogether()
     394             : {
     395         460 :   if (_fe_problem.solverParams(number())._type == Moose::ST_JFNK)
     396           0 :     mooseError(
     397             :         "Evaluting the residual and Jacobian together does not make sense for a JFNK solve type in "
     398             :         "which only function evaluations are required, e.g. there is no need to form a matrix");
     399             : 
     400         460 :   nonlinearSolver()->residual_object = nullptr;
     401         460 :   nonlinearSolver()->jacobian = nullptr;
     402         460 :   nonlinearSolver()->residual_and_jacobian_object = &_resid_and_jac_functor;
     403         460 : }