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      524728 : compute_jacobian(const NumericVector<Number> & soln,
      39             :                  SparseMatrix<Number> & jacobian,
      40             :                  NonlinearImplicitSystem & sys)
      41             : {
      42             :   FEProblemBase * p =
      43      524728 :       sys.get_equation_systems().parameters.get<FEProblemBase *>("_fe_problem_base");
      44      524728 :   p->computeJacobianSys(sys, soln, jacobian);
      45      524708 : }
      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      324273 : compute_nullspace(std::vector<NumericVector<Number> *> & sp, NonlinearImplicitSystem & sys)
      59             : {
      60             :   FEProblemBase * p =
      61      324273 :       sys.get_equation_systems().parameters.get<FEProblemBase *>("_fe_problem_base");
      62      324273 :   p->computeNullSpace(sys, sp);
      63      324273 : }
      64             : 
      65             : void
      66      324273 : compute_transpose_nullspace(std::vector<NumericVector<Number> *> & sp,
      67             :                             NonlinearImplicitSystem & sys)
      68             : {
      69             :   FEProblemBase * p =
      70      324273 :       sys.get_equation_systems().parameters.get<FEProblemBase *>("_fe_problem_base");
      71      324273 :   p->computeTransposeNullSpace(sys, sp);
      72      324273 : }
      73             : 
      74             : void
      75      324273 : compute_nearnullspace(std::vector<NumericVector<Number> *> & sp, NonlinearImplicitSystem & sys)
      76             : {
      77             :   FEProblemBase * p =
      78      324273 :       sys.get_equation_systems().parameters.get<FEProblemBase *>("_fe_problem_base");
      79      324273 :   p->computeNearNullSpace(sys, sp);
      80      324273 : }
      81             : 
      82             : void
      83        1976 : 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        1976 :       sys.get_equation_systems().parameters.get<FEProblemBase *>("_fe_problem_base");
      92        1976 :   p->computePostCheck(
      93             :       sys, old_soln, search_direction, new_soln, changed_search_direction, changed_new_soln);
      94        1976 : }
      95             : } // namespace Moose
      96             : 
      97       62867 : NonlinearSystem::NonlinearSystem(FEProblemBase & fe_problem, const std::string & name)
      98             :   : NonlinearSystemBase(
      99       62867 :         fe_problem, fe_problem.es().add_system<NonlinearImplicitSystem>(name), name),
     100       62867 :     _nl_implicit_sys(fe_problem.es().get_system<NonlinearImplicitSystem>(name)),
     101       62867 :     _nl_residual_functor(_fe_problem),
     102       62867 :     _fd_residual_functor(_fe_problem),
     103       62867 :     _resid_and_jac_functor(_fe_problem),
     104      125734 :     _use_coloring_finite_difference(false)
     105             : {
     106       62867 :   nonlinearSolver()->residual_object = &_nl_residual_functor;
     107       62867 :   nonlinearSolver()->jacobian = Moose::compute_jacobian;
     108       62867 :   nonlinearSolver()->bounds = Moose::compute_bounds;
     109       62867 :   nonlinearSolver()->nullspace = Moose::compute_nullspace;
     110       62867 :   nonlinearSolver()->transpose_nullspace = Moose::compute_transpose_nullspace;
     111       62867 :   nonlinearSolver()->nearnullspace = Moose::compute_nearnullspace;
     112             : 
     113             :   PetscNonlinearSolver<Real> * petsc_solver =
     114       62867 :       static_cast<PetscNonlinearSolver<Real> *>(_nl_implicit_sys.nonlinear_solver.get());
     115       62867 :   if (petsc_solver)
     116             :   {
     117       62867 :     petsc_solver->set_residual_zero_out(false);
     118       62867 :     petsc_solver->set_jacobian_zero_out(false);
     119       62867 :     petsc_solver->use_default_monitor(false);
     120             :   }
     121       62867 : }
     122             : 
     123       58443 : NonlinearSystem::~NonlinearSystem() {}
     124             : 
     125             : void
     126      321948 : NonlinearSystem::potentiallySetupFiniteDifferencing()
     127             : {
     128      321948 :   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      321948 :       static_cast<PetscNonlinearSolver<Real> &>(*_nl_implicit_sys.nonlinear_solver);
     136      321948 :   solver.mffd_residual_object = &_fd_residual_functor;
     137             : 
     138      321948 :   solver.set_snesmf_reuse_base(_fe_problem.useSNESMFReuseBase());
     139      321948 : }
     140             : 
     141             : void
     142      316346 : 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      632291 :   if (_fe_problem.hasDampers() || _fe_problem.shouldUpdateSolution() ||
     149      315945 :       _fe_problem.needsPreviousNewtonIteration())
     150         513 :     _nl_implicit_sys.nonlinear_solver->postcheck = Moose::compute_postcheck;
     151             : 
     152             :   // reset solution invalid counter for the time step
     153      316346 :   if (!_time_integrators.empty())
     154      283664 :     _app.solutionInvalidity().resetSolutionInvalidTimeStep();
     155             : 
     156      316346 :   if (shouldEvaluatePreSMOResidual())
     157             :   {
     158         430 :     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      316346 :   const bool presolve_succeeded = preSolve();
     172      316346 :   if (!presolve_succeeded)
     173           0 :     return;
     174             : 
     175      316346 :   potentiallySetupFiniteDifferencing();
     176             : 
     177      316346 :   const bool time_integrator_solve = std::any_of(_time_integrators.begin(),
     178             :                                                  _time_integrators.end(),
     179      284382 :                                                  [](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      316346 :   if (time_integrator_solve)
     185       12208 :     _time_integrators.front()->solve();
     186             :   else
     187      304138 :     system().solve();
     188             : 
     189      600632 :   for (auto & ti : _time_integrators)
     190             :   {
     191      284362 :     if (!ti->overridesSolve())
     192      272154 :       ti->setNumIterationsLastSolve();
     193      284362 :     ti->postSolve();
     194             :   }
     195             : 
     196      316270 :   if (!_time_integrators.empty())
     197             :   {
     198      283644 :     _n_iters = _time_integrators.front()->getNumNonlinearIterations();
     199      283644 :     _n_linear_iters = _time_integrators.front()->getNumLinearIterations();
     200             :   }
     201             :   else
     202             :   {
     203       32626 :     _n_iters = _nl_implicit_sys.n_nonlinear_iterations();
     204       32626 :     _n_linear_iters = _nl_implicit_sys.nonlinear_solver->get_total_linear_iterations();
     205             :   }
     206             : 
     207             :   // store info about the solve
     208      316270 :   _final_residual = _nl_implicit_sys.final_nonlinear_residual();
     209             : 
     210             :   // Accumulate only the occurence of solution invalid warnings for each time step
     211      316270 :   _app.solutionInvalidity().solutionInvalidAccumulationTimeStep(_fe_problem.timeStep());
     212             : 
     213             :   // determine whether solution invalid occurs in the converged solution
     214      316270 :   checkInvalidSolution();
     215             : 
     216      316266 :   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             : #if PETSC_VERSION_LESS_THAN(3, 24, 0)
     319         194 :   LibmeshPetscCallA(_communicator.get(),
     320             :                     MatFDColoringSetFunction(_fdcoloring,
     321             :                                              (PetscErrorCode(*)(void))(void (*)(void))
     322             :                                              &libMesh::libmesh_petsc_snes_fd_residual,
     323             :                                              &petsc_nonlinear_solver));
     324             : #else
     325             :   LibmeshPetscCallA(_communicator.get(),
     326             :                     MatFDColoringSetFunction(_fdcoloring,
     327             :                                              (MatFDColoringFn*)
     328             :                                              &libMesh::libmesh_petsc_snes_fd_residual,
     329             :                                              &petsc_nonlinear_solver));
     330             : #endif
     331             :   // clang-format on
     332         194 :   LibmeshPetscCallA(_communicator.get(),
     333             :                     MatFDColoringSetUp(petsc_mat->mat(), iscoloring, _fdcoloring));
     334         194 :   LibmeshPetscCallA(_communicator.get(),
     335             :                     SNESSetJacobian(petsc_nonlinear_solver.snes(),
     336             :                                     petsc_mat->mat(),
     337             :                                     petsc_mat->mat(),
     338             :                                     SNESComputeJacobianDefaultColor,
     339             :                                     _fdcoloring));
     340             :   // PETSc >=3.3.0
     341         194 :   LibmeshPetscCallA(_communicator.get(), ISColoringDestroy(&iscoloring));
     342         194 : }
     343             : 
     344             : bool
     345      339527 : NonlinearSystem::converged()
     346             : {
     347      339527 :   if (_fe_problem.hasException() || _fe_problem.getFailNextNonlinearConvergenceCheck())
     348         223 :     return false;
     349      339304 :   if (!_fe_problem.acceptInvalidSolution())
     350             :   {
     351         515 :     mooseWarning("The solution is not converged due to the solution being invalid.");
     352         511 :     return false;
     353             :   }
     354      338789 :   return _nl_implicit_sys.nonlinear_solver->converged;
     355             : }
     356             : 
     357             : void
     358         127 : NonlinearSystem::attachPreconditioner(Preconditioner<Number> * preconditioner)
     359             : {
     360         127 :   nonlinearSolver()->attach_preconditioner(preconditioner);
     361         127 : }
     362             : 
     363             : void
     364         592 : NonlinearSystem::computeScalingJacobian()
     365             : {
     366         592 :   _fe_problem.computeJacobianSys(_nl_implicit_sys, *_current_solution, *_scaling_matrix);
     367         592 : }
     368             : 
     369             : void
     370          40 : NonlinearSystem::computeScalingResidual()
     371             : {
     372          40 :   _fe_problem.computeResidualSys(_nl_implicit_sys, *_current_solution, RHS());
     373          40 : }
     374             : 
     375             : SNES
     376     1195657 : NonlinearSystem::getSNES()
     377             : {
     378             :   PetscNonlinearSolver<Number> * petsc_solver =
     379     1195657 :       dynamic_cast<PetscNonlinearSolver<Number> *>(nonlinearSolver());
     380             : 
     381     1195657 :   if (petsc_solver)
     382             :   {
     383     1195657 :     const char * snes_prefix = nullptr;
     384     1195657 :     std::string snes_prefix_str;
     385     1195657 :     if (system().prefix_with_name())
     386             :     {
     387       42893 :       snes_prefix_str = system().prefix();
     388       42893 :       snes_prefix = snes_prefix_str.c_str();
     389             :     }
     390     2391314 :     return petsc_solver->snes(snes_prefix);
     391     1195657 :   }
     392             :   else
     393           0 :     mooseError("It is not a petsc nonlinear solver");
     394             : }
     395             : 
     396             : void
     397         340 : NonlinearSystem::residualAndJacobianTogether()
     398             : {
     399         340 :   if (_fe_problem.solverParams(number())._type == Moose::ST_JFNK)
     400           0 :     mooseError(
     401             :         "Evaluting the residual and Jacobian together does not make sense for a JFNK solve type in "
     402             :         "which only function evaluations are required, e.g. there is no need to form a matrix");
     403             : 
     404         340 :   nonlinearSolver()->residual_object = nullptr;
     405         340 :   nonlinearSolver()->jacobian = nullptr;
     406         340 :   nonlinearSolver()->residual_and_jacobian_object = &_resid_and_jac_functor;
     407         340 : }