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      475066 : compute_jacobian(const NumericVector<Number> & soln,
      39             :                  SparseMatrix<Number> & jacobian,
      40             :                  NonlinearImplicitSystem & sys)
      41             : {
      42             :   FEProblemBase * p =
      43      475066 :       sys.get_equation_systems().parameters.get<FEProblemBase *>("_fe_problem_base");
      44      475066 :   p->computeJacobianSys(sys, soln, jacobian);
      45      475046 : }
      46             : 
      47             : void
      48         704 : compute_bounds(NumericVector<Number> & lower,
      49             :                NumericVector<Number> & upper,
      50             :                NonlinearImplicitSystem & sys)
      51             : {
      52             :   FEProblemBase * p =
      53         704 :       sys.get_equation_systems().parameters.get<FEProblemBase *>("_fe_problem_base");
      54         704 :   p->computeBounds(sys, lower, upper);
      55         704 : }
      56             : 
      57             : void
      58      292689 : compute_nullspace(std::vector<NumericVector<Number> *> & sp, NonlinearImplicitSystem & sys)
      59             : {
      60             :   FEProblemBase * p =
      61      292689 :       sys.get_equation_systems().parameters.get<FEProblemBase *>("_fe_problem_base");
      62      292689 :   p->computeNullSpace(sys, sp);
      63      292689 : }
      64             : 
      65             : void
      66      292689 : compute_transpose_nullspace(std::vector<NumericVector<Number> *> & sp,
      67             :                             NonlinearImplicitSystem & sys)
      68             : {
      69             :   FEProblemBase * p =
      70      292689 :       sys.get_equation_systems().parameters.get<FEProblemBase *>("_fe_problem_base");
      71      292689 :   p->computeTransposeNullSpace(sys, sp);
      72      292689 : }
      73             : 
      74             : void
      75      292689 : compute_nearnullspace(std::vector<NumericVector<Number> *> & sp, NonlinearImplicitSystem & sys)
      76             : {
      77             :   FEProblemBase * p =
      78      292689 :       sys.get_equation_systems().parameters.get<FEProblemBase *>("_fe_problem_base");
      79      292689 :   p->computeNearNullSpace(sys, sp);
      80      292689 : }
      81             : 
      82             : void
      83        1900 : 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        1900 :       sys.get_equation_systems().parameters.get<FEProblemBase *>("_fe_problem_base");
      92        1900 :   p->computePostCheck(
      93             :       sys, old_soln, search_direction, new_soln, changed_search_direction, changed_new_soln);
      94        1900 : }
      95             : } // namespace Moose
      96             : 
      97       56513 : NonlinearSystem::NonlinearSystem(FEProblemBase & fe_problem, const std::string & name)
      98             :   : NonlinearSystemBase(
      99       56513 :         fe_problem, fe_problem.es().add_system<NonlinearImplicitSystem>(name), name),
     100       56513 :     _nl_implicit_sys(fe_problem.es().get_system<NonlinearImplicitSystem>(name)),
     101       56513 :     _nl_residual_functor(_fe_problem),
     102       56513 :     _fd_residual_functor(_fe_problem),
     103       56513 :     _resid_and_jac_functor(_fe_problem),
     104      113026 :     _use_coloring_finite_difference(false)
     105             : {
     106       56513 :   nonlinearSolver()->residual_object = &_nl_residual_functor;
     107       56513 :   nonlinearSolver()->jacobian = Moose::compute_jacobian;
     108       56513 :   nonlinearSolver()->bounds = Moose::compute_bounds;
     109       56513 :   nonlinearSolver()->nullspace = Moose::compute_nullspace;
     110       56513 :   nonlinearSolver()->transpose_nullspace = Moose::compute_transpose_nullspace;
     111       56513 :   nonlinearSolver()->nearnullspace = Moose::compute_nearnullspace;
     112             : 
     113             :   PetscNonlinearSolver<Real> * petsc_solver =
     114       56513 :       static_cast<PetscNonlinearSolver<Real> *>(_nl_implicit_sys.nonlinear_solver.get());
     115       56513 :   if (petsc_solver)
     116             :   {
     117       56513 :     petsc_solver->set_residual_zero_out(false);
     118       56513 :     petsc_solver->set_jacobian_zero_out(false);
     119       56513 :     petsc_solver->use_default_monitor(false);
     120             :   }
     121       56513 : }
     122             : 
     123       52240 : NonlinearSystem::~NonlinearSystem() {}
     124             : 
     125             : void
     126      290344 : NonlinearSystem::potentiallySetupFiniteDifferencing()
     127             : {
     128      290344 :   if (_use_finite_differenced_preconditioner)
     129             :   {
     130         184 :     _nl_implicit_sys.nonlinear_solver->fd_residual_object = &_fd_residual_functor;
     131         184 :     setupFiniteDifferencedPreconditioner();
     132             :   }
     133             : 
     134             :   PetscNonlinearSolver<Real> & solver =
     135      290344 :       static_cast<PetscNonlinearSolver<Real> &>(*_nl_implicit_sys.nonlinear_solver);
     136      290344 :   solver.mffd_residual_object = &_fd_residual_functor;
     137             : 
     138      290344 :   solver.set_snesmf_reuse_base(_fe_problem.useSNESMFReuseBase());
     139      290344 : }
     140             : 
     141             : void
     142      285236 : 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      570077 :   if (_fe_problem.hasDampers() || _fe_problem.shouldUpdateSolution() ||
     149      284841 :       _fe_problem.needsPreviousNewtonIteration())
     150         489 :     _nl_implicit_sys.nonlinear_solver->postcheck = Moose::compute_postcheck;
     151             : 
     152             :   // reset solution invalid counter for the time step
     153      285236 :   if (!_time_integrators.empty())
     154      255948 :     _app.solutionInvalidity().resetSolutionInvalidTimeStep();
     155             : 
     156      285236 :   if (shouldEvaluatePreSMOResidual())
     157             :   {
     158          77 :     TIME_SECTION("nlPreSMOResidual", 3, "Computing Pre-SMO Residual");
     159             :     // Calculate the pre-SMO residual for use in the convergence criterion.
     160          77 :     _computing_pre_smo_residual = true;
     161          77 :     _fe_problem.computeResidualSys(_nl_implicit_sys, *_current_solution, *_nl_implicit_sys.rhs);
     162          77 :     _computing_pre_smo_residual = false;
     163          77 :     _nl_implicit_sys.rhs->close();
     164          77 :     _pre_smo_residual = _nl_implicit_sys.rhs->l2_norm();
     165          77 :     _console << " * Nonlinear |R| = "
     166         154 :              << Console::outputNorm(std::numeric_limits<Real>::max(), _pre_smo_residual)
     167          77 :              << " (Before preset BCs, predictors, correctors, and constraints)\n";
     168          77 :     _console << std::flush;
     169          77 :   }
     170             : 
     171      285236 :   const bool presolve_succeeded = preSolve();
     172      285236 :   if (!presolve_succeeded)
     173           0 :     return;
     174             : 
     175      285236 :   potentiallySetupFiniteDifferencing();
     176             : 
     177      285236 :   const bool time_integrator_solve = std::any_of(_time_integrators.begin(),
     178             :                                                  _time_integrators.end(),
     179      256662 :                                                  [](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      285236 :   if (time_integrator_solve)
     185       11418 :     _time_integrators.front()->solve();
     186             :   else
     187      273818 :     system().solve();
     188             : 
     189      541804 :   for (auto & ti : _time_integrators)
     190             :   {
     191      256642 :     if (!ti->overridesSolve())
     192      245224 :       ti->setNumIterationsLastSolve();
     193      256642 :     ti->postSolve();
     194             :   }
     195             : 
     196      285162 :   if (!_time_integrators.empty())
     197             :   {
     198      255928 :     _n_iters = _time_integrators.front()->getNumNonlinearIterations();
     199      255928 :     _n_linear_iters = _time_integrators.front()->getNumLinearIterations();
     200             :   }
     201             :   else
     202             :   {
     203       29234 :     _n_iters = _nl_implicit_sys.n_nonlinear_iterations();
     204       29234 :     _n_linear_iters = _nl_implicit_sys.nonlinear_solver->get_total_linear_iterations();
     205             :   }
     206             : 
     207             :   // store info about the solve
     208      285162 :   _final_residual = _nl_implicit_sys.final_nonlinear_residual();
     209             : 
     210             :   // Accumulate only the occurence of solution invalid warnings for each time step
     211      285162 :   _app.solutionInvalidity().solutionInvalidAccumulationTimeStep(_fe_problem.timeStep());
     212             : 
     213             :   // determine whether solution invalid occurs in the converged solution
     214      285162 :   checkInvalidSolution();
     215             : 
     216      285158 :   if (_use_coloring_finite_difference)
     217          69 :     LibmeshPetscCall(MatFDColoringDestroy(&_fdcoloring));
     218             : }
     219             : 
     220             : void
     221         321 : NonlinearSystem::stopSolve(const ExecFlagType & exec_flag,
     222             :                            const std::set<TagID> & vector_tags_to_close)
     223             : {
     224             :   PetscNonlinearSolver<Real> & solver =
     225         321 :       static_cast<PetscNonlinearSolver<Real> &>(*sys().nonlinear_solver);
     226             : 
     227         321 :   if (exec_flag == EXEC_LINEAR || exec_flag == EXEC_POSTCHECK)
     228             :   {
     229         231 :     LibmeshPetscCall(SNESSetFunctionDomainError(solver.snes()));
     230             : 
     231             :     // Clean up by getting vectors into a valid state for a
     232             :     // (possible) subsequent solve.
     233         231 :     closeTaggedVectors(vector_tags_to_close);
     234             :   }
     235          90 :   else if (exec_flag == EXEC_NONLINEAR)
     236          90 :     LibmeshPetscCall(SNESSetJacobianDomainError(solver.snes()));
     237             :   else
     238           0 :     mooseError("Unsupported execute flag: ", Moose::stringify(exec_flag));
     239         321 : }
     240             : 
     241             : void
     242         184 : NonlinearSystem::setupFiniteDifferencedPreconditioner()
     243             : {
     244             :   std::shared_ptr<FiniteDifferencePreconditioner> fdp =
     245         184 :       std::dynamic_pointer_cast<FiniteDifferencePreconditioner>(_preconditioner);
     246         184 :   if (!fdp)
     247           0 :     mooseError("Did not setup finite difference preconditioner, and please add a preconditioning "
     248             :                "block with type = fdp");
     249             : 
     250         184 :   if (fdp->finiteDifferenceType() == "coloring")
     251             :   {
     252         168 :     _use_coloring_finite_difference = true;
     253         168 :     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         184 : }
     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         168 : NonlinearSystem::setupColoringFiniteDifferencedPreconditioner()
     286             : {
     287             :   // Make sure that libMesh isn't going to override our preconditioner
     288         168 :   _nl_implicit_sys.nonlinear_solver->jacobian = nullptr;
     289             : 
     290             :   PetscNonlinearSolver<Number> & petsc_nonlinear_solver =
     291         168 :       dynamic_cast<PetscNonlinearSolver<Number> &>(*_nl_implicit_sys.nonlinear_solver);
     292             : 
     293             :   // Pointer to underlying PetscMatrix type
     294             :   PetscMatrix<Number> * petsc_mat =
     295         168 :       dynamic_cast<PetscMatrix<Number> *>(&_nl_implicit_sys.get_system_matrix());
     296             : 
     297         168 :   Moose::compute_jacobian(*_nl_implicit_sys.current_local_solution, *petsc_mat, _nl_implicit_sys);
     298             : 
     299         168 :   if (!petsc_mat)
     300           0 :     mooseError("Could not convert to Petsc matrix.");
     301             : 
     302         168 :   petsc_mat->close();
     303             : 
     304             :   ISColoring iscoloring;
     305             : 
     306             :   // PETSc 3.5.x
     307             :   MatColoring matcoloring;
     308         168 :   LibmeshPetscCallA(_communicator.get(), MatColoringCreate(petsc_mat->mat(), &matcoloring));
     309         168 :   LibmeshPetscCallA(_communicator.get(), MatColoringSetType(matcoloring, MATCOLORINGLF));
     310         168 :   LibmeshPetscCallA(_communicator.get(), MatColoringSetFromOptions(matcoloring));
     311         168 :   LibmeshPetscCallA(_communicator.get(), MatColoringApply(matcoloring, &iscoloring));
     312         168 :   LibmeshPetscCallA(_communicator.get(), MatColoringDestroy(&matcoloring));
     313             : 
     314         168 :   LibmeshPetscCallA(_communicator.get(),
     315             :                     MatFDColoringCreate(petsc_mat->mat(), iscoloring, &_fdcoloring));
     316         168 :   LibmeshPetscCallA(_communicator.get(), MatFDColoringSetFromOptions(_fdcoloring));
     317             :   // clang-format off
     318         168 :   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         168 :   LibmeshPetscCallA(_communicator.get(),
     324             :                     MatFDColoringSetUp(petsc_mat->mat(), iscoloring, _fdcoloring));
     325         168 :   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         168 :   LibmeshPetscCallA(_communicator.get(), ISColoringDestroy(&iscoloring));
     333         168 : }
     334             : 
     335             : bool
     336      304097 : NonlinearSystem::converged()
     337             : {
     338      304097 :   if (_fe_problem.hasException() || _fe_problem.getFailNextNonlinearConvergenceCheck())
     339         207 :     return false;
     340      303890 :   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      303375 :   return _nl_implicit_sys.nonlinear_solver->converged;
     346             : }
     347             : 
     348             : void
     349         106 : NonlinearSystem::attachPreconditioner(Preconditioner<Number> * preconditioner)
     350             : {
     351         106 :   nonlinearSolver()->attach_preconditioner(preconditioner);
     352         106 : }
     353             : 
     354             : void
     355         584 : NonlinearSystem::computeScalingJacobian()
     356             : {
     357         584 :   _fe_problem.computeJacobianSys(_nl_implicit_sys, *_current_solution, *_scaling_matrix);
     358         584 : }
     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     1079804 : NonlinearSystem::getSNES()
     368             : {
     369             :   PetscNonlinearSolver<Number> * petsc_solver =
     370     1079804 :       dynamic_cast<PetscNonlinearSolver<Number> *>(nonlinearSolver());
     371             : 
     372     1079804 :   if (petsc_solver)
     373             :   {
     374     1079804 :     const char * snes_prefix = nullptr;
     375     1079804 :     std::string snes_prefix_str;
     376     1079804 :     if (system().prefix_with_name())
     377             :     {
     378       38587 :       snes_prefix_str = system().prefix();
     379       38587 :       snes_prefix = snes_prefix_str.c_str();
     380             :     }
     381     2159608 :     return petsc_solver->snes(snes_prefix);
     382     1079804 :   }
     383             :   else
     384           0 :     mooseError("It is not a petsc nonlinear solver");
     385             : }
     386             : 
     387             : void
     388         319 : NonlinearSystem::residualAndJacobianTogether()
     389             : {
     390         319 :   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         319 :   nonlinearSolver()->residual_object = nullptr;
     396         319 :   nonlinearSolver()->jacobian = nullptr;
     397         319 :   nonlinearSolver()->residual_and_jacobian_object = &_resid_and_jac_functor;
     398         319 : }