https://mooseframework.inl.gov
NonlinearSystem.C
Go to the documentation of this file.
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"
16 #include "PetscSupport.h"
17 #include "ComputeResidualFunctor.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
39  SparseMatrix<Number> & jacobian,
41 {
42  FEProblemBase * p =
43  sys.get_equation_systems().parameters.get<FEProblemBase *>("_fe_problem_base");
44  p->computeJacobianSys(sys, soln, jacobian);
45 }
46 
47 void
49  NumericVector<Number> & upper,
51 {
52  FEProblemBase * p =
53  sys.get_equation_systems().parameters.get<FEProblemBase *>("_fe_problem_base");
54  p->computeBounds(sys, lower, upper);
55 }
56 
57 void
59 {
60  FEProblemBase * p =
61  sys.get_equation_systems().parameters.get<FEProblemBase *>("_fe_problem_base");
62  p->computeNullSpace(sys, sp);
63 }
64 
65 void
68 {
69  FEProblemBase * p =
70  sys.get_equation_systems().parameters.get<FEProblemBase *>("_fe_problem_base");
71  p->computeTransposeNullSpace(sys, sp);
72 }
73 
74 void
76 {
77  FEProblemBase * p =
78  sys.get_equation_systems().parameters.get<FEProblemBase *>("_fe_problem_base");
79  p->computeNearNullSpace(sys, sp);
80 }
81 
82 void
84  NumericVector<Number> & search_direction,
85  NumericVector<Number> & new_soln,
86  bool & changed_search_direction,
87  bool & changed_new_soln,
89 {
90  FEProblemBase * p =
91  sys.get_equation_systems().parameters.get<FEProblemBase *>("_fe_problem_base");
93  sys, old_soln, search_direction, new_soln, changed_search_direction, changed_new_soln);
94 }
95 } // namespace Moose
96 
97 NonlinearSystem::NonlinearSystem(FEProblemBase & fe_problem, const std::string & name)
99  fe_problem, fe_problem.es().add_system<NonlinearImplicitSystem>(name), name),
100  _nl_implicit_sys(fe_problem.es().get_system<NonlinearImplicitSystem>(name)),
101  _nl_residual_functor(_fe_problem),
102  _fd_residual_functor(_fe_problem),
103  _resid_and_jac_functor(_fe_problem),
104  _use_coloring_finite_difference(false)
105 {
106  nonlinearSolver()->residual_object = &_nl_residual_functor;
110  nonlinearSolver()->transpose_nullspace = Moose::compute_transpose_nullspace;
112 
113  PetscNonlinearSolver<Real> * petsc_solver =
115  if (petsc_solver)
116  {
117  petsc_solver->set_residual_zero_out(false);
118  petsc_solver->set_jacobian_zero_out(false);
119  petsc_solver->use_default_monitor(false);
120  }
121 }
122 
124 
125 void
127 {
129  {
132  }
133 
134  PetscNonlinearSolver<Real> & solver =
137 
139 }
140 
141 void
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.
151 
152  // reset solution invalid counter for the time step
153  if (!_time_integrators.empty())
155 
157  {
158  TIME_SECTION("nlPreSMOResidual", 3, "Computing Pre-SMO Residual");
159  // Calculate the pre-SMO residual for use in the convergence criterion.
165  _console << " * Nonlinear |R| = "
167  << " (Before preset BCs, predictors, correctors, and constraints)\n";
168  _console << std::flush;
169  }
170 
171  const bool presolve_succeeded = preSolve();
172  if (!presolve_succeeded)
173  return;
174 
176 
177  const bool time_integrator_solve = std::any_of(_time_integrators.begin(),
178  _time_integrators.end(),
179  [](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  if (time_integrator_solve)
185  _time_integrators.front()->solve();
186  else
187  system().solve();
188 
189  for (auto & ti : _time_integrators)
190  {
191  if (!ti->overridesSolve())
192  ti->setNumIterationsLastSolve();
193  ti->postSolve();
194  }
195 
196  if (!_time_integrators.empty())
197  {
198  _n_iters = _time_integrators.front()->getNumNonlinearIterations();
199  _n_linear_iters = _time_integrators.front()->getNumLinearIterations();
200  }
201  else
202  {
204  _n_linear_iters = _nl_implicit_sys.nonlinear_solver->get_total_linear_iterations();
205  }
206 
207  // store info about the solve
209 
210  // Accumulate only the occurence of solution invalid warnings for each time step
212 
213  // determine whether solution invalid occurs in the converged solution
215 
217  LibmeshPetscCall(MatFDColoringDestroy(&_fdcoloring));
218 }
219 
220 void
222  const std::set<TagID> & vector_tags_to_close)
223 {
224  PetscNonlinearSolver<Real> & solver =
226 
227  if (exec_flag == EXEC_LINEAR || exec_flag == EXEC_POSTCHECK)
228  {
229  LibmeshPetscCall(SNESSetFunctionDomainError(solver.snes()));
230 
231  // Clean up by getting vectors into a valid state for a
232  // (possible) subsequent solve.
233  closeTaggedVectors(vector_tags_to_close);
234  }
235  else if (exec_flag == EXEC_NONLINEAR)
236  LibmeshPetscCall(SNESSetJacobianDomainError(solver.snes()));
237  else
238  mooseError("Unsupported execute flag: ", Moose::stringify(exec_flag));
239 }
240 
241 void
243 {
244  std::shared_ptr<FiniteDifferencePreconditioner> fdp =
246  if (!fdp)
247  mooseError("Did not setup finite difference preconditioner, and please add a preconditioning "
248  "block with type = fdp");
249 
250  if (fdp->finiteDifferenceType() == "coloring")
251  {
254  }
255 
256  else if (fdp->finiteDifferenceType() == "standard")
257  {
260  }
261  else
262  mooseError("Unknown finite difference type");
263 }
264 
265 void
267 {
268  // Make sure that libMesh isn't going to override our preconditioner
269  _nl_implicit_sys.nonlinear_solver->jacobian = nullptr;
270 
271  PetscNonlinearSolver<Number> * petsc_nonlinear_solver =
273 
274  PetscMatrix<Number> * petsc_mat =
276 
277  LibmeshPetscCall(SNESSetJacobian(petsc_nonlinear_solver->snes(),
278  petsc_mat->mat(),
279  petsc_mat->mat(),
280  SNESComputeJacobianDefault,
281  nullptr));
282 }
283 
284 void
286 {
287  // Make sure that libMesh isn't going to override our preconditioner
288  _nl_implicit_sys.nonlinear_solver->jacobian = nullptr;
289 
290  PetscNonlinearSolver<Number> & petsc_nonlinear_solver =
292 
293  // Pointer to underlying PetscMatrix type
294  PetscMatrix<Number> * petsc_mat =
296 
298 
299  if (!petsc_mat)
300  mooseError("Could not convert to Petsc matrix.");
301 
302  petsc_mat->close();
303 
304  ISColoring iscoloring;
305 
306  // PETSc 3.5.x
307  MatColoring matcoloring;
308  LibmeshPetscCallA(_communicator.get(), MatColoringCreate(petsc_mat->mat(), &matcoloring));
309  LibmeshPetscCallA(_communicator.get(), MatColoringSetType(matcoloring, MATCOLORINGLF));
310  LibmeshPetscCallA(_communicator.get(), MatColoringSetFromOptions(matcoloring));
311  LibmeshPetscCallA(_communicator.get(), MatColoringApply(matcoloring, &iscoloring));
312  LibmeshPetscCallA(_communicator.get(), MatColoringDestroy(&matcoloring));
313 
314  LibmeshPetscCallA(_communicator.get(),
315  MatFDColoringCreate(petsc_mat->mat(), iscoloring, &_fdcoloring));
316  LibmeshPetscCallA(_communicator.get(), MatFDColoringSetFromOptions(_fdcoloring));
317  // clang-format off
318  LibmeshPetscCallA(_communicator.get(), MatFDColoringSetFunction(_fdcoloring,
319  (PetscErrorCode(*)(void))(void (*)(void)) &
321  &petsc_nonlinear_solver));
322  // clang-format on
323  LibmeshPetscCallA(_communicator.get(),
324  MatFDColoringSetUp(petsc_mat->mat(), iscoloring, _fdcoloring));
325  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  LibmeshPetscCallA(_communicator.get(), ISColoringDestroy(&iscoloring));
333 }
334 
335 bool
337 {
339  return false;
341  {
342  mooseWarning("The solution is not converged due to the solution being invalid.");
343  return false;
344  }
345  return _nl_implicit_sys.nonlinear_solver->converged;
346 }
347 
348 void
350 {
351  nonlinearSolver()->attach_preconditioner(preconditioner);
352 }
353 
354 void
356 {
358 }
359 
360 void
362 {
364 }
365 
366 SNES
368 {
369  PetscNonlinearSolver<Number> * petsc_solver =
371 
372  if (petsc_solver)
373  {
374  const char * snes_prefix = nullptr;
375  std::string snes_prefix_str;
376  if (system().prefix_with_name())
377  {
378  snes_prefix_str = system().prefix();
379  snes_prefix = snes_prefix_str.c_str();
380  }
381  return petsc_solver->snes(snes_prefix);
382  }
383  else
384  mooseError("It is not a petsc nonlinear solver");
385 }
386 
387 void
389 {
391  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  nonlinearSolver()->residual_object = nullptr;
396  nonlinearSolver()->jacobian = nullptr;
397  nonlinearSolver()->residual_and_jacobian_object = &_resid_and_jac_functor;
398 }
std::string name(const ElemQuality q)
std::vector< std::shared_ptr< TimeIntegrator > > _time_integrators
Time integrator.
Definition: SystemBase.h:1041
void compute_jacobian(const NumericVector< Number > &soln, SparseMatrix< Number > &jacobian, NonlinearImplicitSystem &sys)
SNES snes(const char *name=nullptr)
virtual libMesh::NonlinearSolver< Number > * nonlinearSolver() override
std::unique_ptr< NonlinearSolver< Number > > nonlinear_solver
virtual void solve() override
Solve the system (using libMesh magic)
void solutionInvalidAccumulationTimeStep(const unsigned int timestep_index)
Pass the number of solution invalid occurrences from current iteration to cumulative time iteration c...
PetscErrorCode libmesh_petsc_snes_fd_residual(SNES, Vec x, Vec r, void *ctx)
libMesh::NonlinearImplicitSystem & _nl_implicit_sys
void use_default_monitor(bool state)
void checkInvalidSolution()
Definition: SolverSystem.C:111
void mooseError(Args &&... args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Definition: MooseError.h:302
void compute_nearnullspace(std::vector< NumericVector< Number > *> &sp, NonlinearImplicitSystem &sys)
void mooseWarning(Args &&... args)
Emit a warning message with the given stringified, concatenated args.
Definition: MooseError.h:336
const EquationSystems & get_equation_systems() const
void resetSolutionInvalidTimeStep()
Reset the number of solution invalid occurrences back to zero for the current time step...
NonlinearSystem(FEProblemBase &problem, const std::string &name)
virtual void setupFiniteDifferencedPreconditioner() override
NumericVector< Number > * rhs
ComputeResidualFunctor _nl_residual_functor
void compute_bounds(NumericVector< Number > &lower, NumericVector< Number > &upper, NonlinearImplicitSystem &sys)
virtual void computeNearNullSpace(libMesh::NonlinearImplicitSystem &sys, std::vector< NumericVector< libMesh::Number > *> &sp)
std::unique_ptr< T_DEST, T_DELETER > dynamic_pointer_cast(std::unique_ptr< T_SRC, T_DELETER > &src)
These are reworked from https://stackoverflow.com/a/11003103.
unsigned int n_nonlinear_iterations() const
std::unique_ptr< libMesh::DiagonalMatrix< Number > > _scaling_matrix
A diagonal matrix used for computing scaling.
const Parallel::Communicator & _communicator
virtual void residualAndJacobianTogether() override
Call this method if you want the residual and Jacobian to be computed simultaneously.
The following methods are specializations for using the libMesh::Parallel::packed_range_* routines fo...
Real _pre_smo_residual
The pre-SMO residual, see setPreSMOResidual for a detailed explanation.
bool hasDampers()
Whether or not this system has dampers.
Specialization of SubProblem for solving nonlinear equations plus auxiliary equations.
virtual void computeJacobianSys(libMesh::NonlinearImplicitSystem &sys, const NumericVector< libMesh::Number > &soln, libMesh::SparseMatrix< libMesh::Number > &jacobian)
Form a Jacobian matrix.
auto max(const L &left, const R &right)
Nonlinear system to be solved.
void compute_postcheck(const NumericVector< Number > &old_soln, NumericVector< Number > &search_direction, NumericVector< Number > &new_soln, bool &changed_search_direction, bool &changed_new_soln, NonlinearImplicitSystem &sys)
virtual void stopSolve(const ExecFlagType &exec_flag, const std::set< TagID > &vector_tags_to_close) override
Quit the current solve as soon as possible.
virtual void computeResidualSys(libMesh::NonlinearImplicitSystem &sys, const NumericVector< libMesh::Number > &soln, NumericVector< libMesh::Number > &residual)
This function is called by Libmesh to form a residual.
const T & get(std::string_view) const
virtual Real l2_norm() const=0
ComputeFDResidualFunctor _fd_residual_functor
virtual bool shouldUpdateSolution()
Check to see whether the problem should update the solution.
std::string prefix() const
ComputeResidualAndJacobian _resid_and_jac_functor
SolutionInvalidity & solutionInvalidity()
Get the SolutionInvalidity for this app.
Definition: MooseApp.h:178
std::shared_ptr< MoosePreconditioner > _preconditioner
Preconditioner.
Jacobian-Free Newton Krylov.
Definition: MooseTypes.h:845
bool _use_finite_differenced_preconditioner
Whether or not to use a finite differenced preconditioner.
void needsPreviousNewtonIteration(bool state)
Set a flag that indicated that user required values for the previous Newton iterate.
void setupColoringFiniteDifferencedPreconditioner()
According to the nonzero pattern provided in the matrix, a graph is constructed.
bool shouldEvaluatePreSMOResidual() const
We offer the option to check convergence against the pre-SMO residual.
Moose::SolveType _type
Definition: SolverParams.h:19
virtual void computeNullSpace(libMesh::NonlinearImplicitSystem &sys, std::vector< NumericVector< libMesh::Number > *> &sp)
void compute_transpose_nullspace(std::vector< NumericVector< Number > *> &sp, NonlinearImplicitSystem &sys)
bool _use_coloring_finite_difference
void computeScalingResidual() override
Compute a "residual" for automatic scaling purposes.
virtual libMesh::NonlinearImplicitSystem & sys()
unsigned int number() const
Gets the number of this system.
Definition: SystemBase.C:1159
const ExecFlagType EXEC_LINEAR
Definition: Moose.C:29
std::string stringify(const T &t)
conversion to string
Definition: Conversion.h:64
void closeTaggedVectors(const std::set< TagID > &tags)
Close all vectors for given tags.
Definition: SystemBase.C:659
virtual NumericVector< Number > & RHS() override
virtual void solve()
virtual void close()=0
virtual int & timeStep() const
virtual bool converged() override
Returns the convergence state.
const NumericVector< Number > * _current_solution
solution vector from solver
Definition: SolverSystem.h:105
void set_residual_zero_out(bool state)
void set_jacobian_zero_out(bool state)
virtual SNES getSNES() override
const ExecFlagType EXEC_POSTCHECK
Definition: Moose.C:33
const ExecFlagType EXEC_NONLINEAR
Definition: Moose.C:31
MooseApp & _app
Definition: SystemBase.h:982
FEProblemBase & _fe_problem
the governing finite element/volume problem
Definition: SystemBase.h:980
Finite difference preconditioner.
virtual void attachPreconditioner(libMesh::Preconditioner< Number > *preconditioner) override
Attach a customized preconditioner that requires physics knowledge.
virtual void potentiallySetupFiniteDifferencing() override
Create finite differencing contexts for assembly of the Jacobian and/or approximating the action of t...
Class for containing MooseEnum item information.
Definition: MooseEnumItem.h:18
static std::string outputNorm(const Real &old_norm, const Real &norm, const unsigned int precision=6)
A helper function for outputting norms in color.
Definition: Console.C:619
bool acceptInvalidSolution() const
Whether or not to accept the solution based on its invalidity.
void compute_nullspace(std::vector< NumericVector< Number > *> &sp, NonlinearImplicitSystem &sys)
if(!dmm->_nl) SETERRQ(PETSC_COMM_WORLD
virtual void computeBounds(libMesh::NonlinearImplicitSystem &sys, NumericVector< libMesh::Number > &lower, NumericVector< libMesh::Number > &upper)
virtual void computePostCheck(libMesh::NonlinearImplicitSystem &sys, const NumericVector< libMesh::Number > &old_soln, NumericVector< libMesh::Number > &search_direction, NumericVector< libMesh::Number > &new_soln, bool &changed_search_direction, bool &changed_new_soln)
void set_snesmf_reuse_base(bool state)
bool getFailNextNonlinearConvergenceCheck() const
Whether it will skip further residual evaluations and fail the next nonlinear convergence check(s) ...
NonlinearImplicitSystem::ComputeResidual * mffd_residual_object
SolverParams & solverParams(unsigned int solver_sys_num=0)
Get the solver parameters.
std::unique_ptr< NumericVector< Number > > current_local_solution
virtual void computeTransposeNullSpace(libMesh::NonlinearImplicitSystem &sys, std::vector< NumericVector< libMesh::Number > *> &sp)
void setupStandardFiniteDifferencedPreconditioner()
Form preconditioning matrix via a standard finite difference method column-by-column.
void computeScalingJacobian() override
Compute a "Jacobian" for automatic scaling purposes.
bool useSNESMFReuseBase()
Return a flag that indicates if we are reusing the vector base.
MOOSE now contains C++17 code, so give a reasonable error message stating what the user can do to add...
const ConsoleStream _console
An instance of helper class to write streams to the Console objects.
bool preSolve()
Perform some steps to get ready for the solver.
void prefix_with_name(bool value)
virtual bool hasException()
Whether or not an exception has occurred.
const SparseMatrix< Number > & get_system_matrix() const
virtual ~NonlinearSystem()
virtual libMesh::System & system() override
Get the reference to the libMesh system.