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