www.mooseframework.org
NonlinearSystem.C
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://www.mooseframework.org
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 
23 #include "libmesh/nonlinear_solver.h"
24 #include "libmesh/petsc_nonlinear_solver.h"
25 #include "libmesh/sparse_matrix.h"
26 #include "libmesh/petsc_matrix.h"
27 #include "libmesh/diagonal_matrix.h"
28 #include "libmesh/default_coupling.h"
29 
30 namespace Moose
31 {
32 void
34  SparseMatrix<Number> & jacobian,
35  NonlinearImplicitSystem & sys)
36 {
37  FEProblemBase * p =
38  sys.get_equation_systems().parameters.get<FEProblemBase *>("_fe_problem_base");
39  p->computeJacobianSys(sys, soln, jacobian);
40 }
41 
42 void
44  NumericVector<Number> & upper,
45  NonlinearImplicitSystem & sys)
46 {
47  FEProblemBase * p =
48  sys.get_equation_systems().parameters.get<FEProblemBase *>("_fe_problem_base");
49  p->computeBounds(sys, lower, upper);
50 }
51 
52 void
53 compute_nullspace(std::vector<NumericVector<Number> *> & sp, NonlinearImplicitSystem & sys)
54 {
55  FEProblemBase * p =
56  sys.get_equation_systems().parameters.get<FEProblemBase *>("_fe_problem_base");
57  p->computeNullSpace(sys, sp);
58 }
59 
60 void
62  NonlinearImplicitSystem & sys)
63 {
64  FEProblemBase * p =
65  sys.get_equation_systems().parameters.get<FEProblemBase *>("_fe_problem_base");
66  p->computeTransposeNullSpace(sys, sp);
67 }
68 
69 void
70 compute_nearnullspace(std::vector<NumericVector<Number> *> & sp, NonlinearImplicitSystem & sys)
71 {
72  FEProblemBase * p =
73  sys.get_equation_systems().parameters.get<FEProblemBase *>("_fe_problem_base");
74  p->computeNearNullSpace(sys, sp);
75 }
76 
77 void
79  NumericVector<Number> & search_direction,
80  NumericVector<Number> & new_soln,
81  bool & changed_search_direction,
82  bool & changed_new_soln,
83  NonlinearImplicitSystem & sys)
84 {
85  FEProblemBase * p =
86  sys.get_equation_systems().parameters.get<FEProblemBase *>("_fe_problem_base");
88  sys, old_soln, search_direction, new_soln, changed_search_direction, changed_new_soln);
89 }
90 } // namespace Moose
91 
92 NonlinearSystem::NonlinearSystem(FEProblemBase & fe_problem, const std::string & name)
94  fe_problem, fe_problem.es().add_system<NonlinearImplicitSystem>(name), name),
95  _nl_implicit_sys(fe_problem.es().get_system<NonlinearImplicitSystem>(name)),
96  _nl_residual_functor(_fe_problem),
97  _fd_residual_functor(_fe_problem),
98  _resid_and_jac_functor(_fe_problem),
99  _use_coloring_finite_difference(false),
100  _solution_is_invalid(false)
101 {
108 
109  PetscNonlinearSolver<Real> * petsc_solver =
110  static_cast<PetscNonlinearSolver<Real> *>(_nl_implicit_sys.nonlinear_solver.get());
111  if (petsc_solver)
112  {
113  petsc_solver->set_residual_zero_out(false);
114  petsc_solver->set_jacobian_zero_out(false);
115  petsc_solver->use_default_monitor(false);
116  }
117 }
118 
120 
121 void
123 {
125 
126  if (_automatic_scaling && _resid_vs_jac_scaling_param < 1. - TOLERANCE)
127  // Add diagonal matrix that will be used for computing scaling factors
128  _nl_implicit_sys.add_matrix<DiagonalMatrix>("scaling_matrix");
129 }
130 
131 void
133 {
134  // Only attach the postcheck function to the solver if we actually
135  // have dampers or if the FEProblemBase needs to update the solution,
136  // which is also done during the linesearch postcheck. It doesn't
137  // hurt to do this multiple times, it is just setting a pointer.
140  _nl_implicit_sys.nonlinear_solver->postcheck = Moose::compute_postcheck;
141 
143  {
144  TIME_SECTION("nlInitialResidual", 3, "Computing Initial Residual");
145  // Calculate the initial residual for use in the convergence criterion.
149  _nl_implicit_sys.rhs->close();
152  _console << "Initial residual before setting preset BCs: "
153  << _initial_residual_before_preset_bcs << std::endl;
154  }
155 
156  const bool presolve_succeeded = preSolve();
157  if (!presolve_succeeded)
158  return;
159 
161  {
162  _nl_implicit_sys.nonlinear_solver->fd_residual_object = &_fd_residual_functor;
164  }
165 
166  PetscNonlinearSolver<Real> & solver =
167  static_cast<PetscNonlinearSolver<Real> &>(*_nl_implicit_sys.nonlinear_solver);
168  solver.mffd_residual_object = &_fd_residual_functor;
169 
170  solver.set_snesmf_reuse_base(_fe_problem.useSNESMFReuseBase());
171 
172  if (_time_integrator)
173  {
174  // reset solution invalid counter for the time step
176  _time_integrator->solve();
177  _time_integrator->postSolve();
178  _n_iters = _time_integrator->getNumNonlinearIterations();
179  _n_linear_iters = _time_integrator->getNumLinearIterations();
180  // Accumulate only the occurence of solution invalid warnings for the current time step counters
182  }
183  else
184  {
185  system().solve();
186  _n_iters = _nl_implicit_sys.n_nonlinear_iterations();
187  _n_linear_iters = solver.get_total_linear_iterations();
188  }
189 
190  // store info about the solve
191  _final_residual = _nl_implicit_sys.final_nonlinear_residual();
192 
193  // determine whether solution invalid occurs in the converged solution
195 
196  // output the solution invalid summary
198  {
199  // sync all solution invalid counts to rank 0 process
201 
203  mooseWarning("The Solution Invalidity warnings are detected but silenced! "
204  "Use Problem/allow_invalid_solution=false to activate ");
205  else
206  // output the occurrence of solution invalid in a summary table
208  }
209 
211  MatFDColoringDestroy(&_fdcoloring);
212 }
213 
214 void
216 {
217  PetscNonlinearSolver<Real> & solver =
218  static_cast<PetscNonlinearSolver<Real> &>(*sys().nonlinear_solver);
219 
220  if (exec_flag == EXEC_LINEAR || exec_flag == EXEC_POSTCHECK)
221  {
222  SNESSetFunctionDomainError(solver.snes());
223 
224  // Clean up by getting vectors into a valid state for a
225  // (possible) subsequent solve. There may be more than just
226  // these...
227  _nl_implicit_sys.rhs->close();
228  if (_Re_time)
229  _Re_time->close();
230  _Re_non_time->close();
231  }
232  else if (exec_flag == EXEC_NONLINEAR)
233  SNESSetJacobianDomainError(solver.snes());
234  else
235  mooseError("Unsupported execute flag: ", Moose::stringify(exec_flag));
236 }
237 
238 void
240 {
241  std::shared_ptr<FiniteDifferencePreconditioner> fdp =
243  if (!fdp)
244  mooseError("Did not setup finite difference preconditioner, and please add a preconditioning "
245  "block with type = fdp");
246 
247  if (fdp->finiteDifferenceType() == "coloring")
248  {
251  }
252 
253  else if (fdp->finiteDifferenceType() == "standard")
254  {
257  }
258  else
259  mooseError("Unknown finite difference type");
260 }
261 
262 void
264 {
265  // Make sure that libMesh isn't going to override our preconditioner
266  _nl_implicit_sys.nonlinear_solver->jacobian = nullptr;
267 
268  PetscNonlinearSolver<Number> * petsc_nonlinear_solver =
269  static_cast<PetscNonlinearSolver<Number> *>(_nl_implicit_sys.nonlinear_solver.get());
270 
271  PetscMatrix<Number> * petsc_mat =
272  static_cast<PetscMatrix<Number> *>(&_nl_implicit_sys.get_system_matrix());
273 
274  SNESSetJacobian(petsc_nonlinear_solver->snes(),
275  petsc_mat->mat(),
276  petsc_mat->mat(),
277  SNESComputeJacobianDefault,
278  nullptr);
279 }
280 
281 void
283 {
284  // Make sure that libMesh isn't going to override our preconditioner
285  _nl_implicit_sys.nonlinear_solver->jacobian = nullptr;
286 
287  PetscNonlinearSolver<Number> & petsc_nonlinear_solver =
288  dynamic_cast<PetscNonlinearSolver<Number> &>(*_nl_implicit_sys.nonlinear_solver);
289 
290  // Pointer to underlying PetscMatrix type
291  PetscMatrix<Number> * petsc_mat =
292  dynamic_cast<PetscMatrix<Number> *>(&_nl_implicit_sys.get_system_matrix());
293 
294  Moose::compute_jacobian(*_nl_implicit_sys.current_local_solution, *petsc_mat, _nl_implicit_sys);
295 
296  if (!petsc_mat)
297  mooseError("Could not convert to Petsc matrix.");
298 
299  petsc_mat->close();
300 
301  PetscErrorCode ierr = 0;
302  ISColoring iscoloring;
303 
304  // PETSc 3.5.x
305  MatColoring matcoloring;
306  ierr = MatColoringCreate(petsc_mat->mat(), &matcoloring);
307  CHKERRABORT(_communicator.get(), ierr);
308  ierr = MatColoringSetType(matcoloring, MATCOLORINGLF);
309  CHKERRABORT(_communicator.get(), ierr);
310  ierr = MatColoringSetFromOptions(matcoloring);
311  CHKERRABORT(_communicator.get(), ierr);
312  ierr = MatColoringApply(matcoloring, &iscoloring);
313  CHKERRABORT(_communicator.get(), ierr);
314  ierr = MatColoringDestroy(&matcoloring);
315  CHKERRABORT(_communicator.get(), ierr);
316 
317  MatFDColoringCreate(petsc_mat->mat(), iscoloring, &_fdcoloring);
318  MatFDColoringSetFromOptions(_fdcoloring);
319  // clang-format off
320  MatFDColoringSetFunction(_fdcoloring,
321  (PetscErrorCode(*)(void))(void (*)(void)) &
323  &petsc_nonlinear_solver);
324  // clang-format on
325  MatFDColoringSetUp(petsc_mat->mat(), iscoloring, _fdcoloring);
326  SNESSetJacobian(petsc_nonlinear_solver.snes(),
327  petsc_mat->mat(),
328  petsc_mat->mat(),
329  SNESComputeJacobianDefaultColor,
330  _fdcoloring);
331  // PETSc >=3.3.0
332  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  return petsc_solver->snes();
374  else
375  mooseError("It is not a petsc nonlinear solver");
376 }
377 
378 void
380 {
382  mooseError(
383  "Evaluting the residual and Jacobian together does not make sense for a JFNK solve type in "
384  "which only function evaluations are required, e.g. there is no need to form a matrix");
385 
386  nonlinearSolver()->residual_object = nullptr;
387  nonlinearSolver()->jacobian = nullptr;
389 }
std::string name(const ElemQuality q)
void compute_jacobian(const NumericVector< Number > &soln, SparseMatrix< Number > &jacobian, NonlinearImplicitSystem &sys)
std::unique_ptr< DiagonalMatrix< Number > > _scaling_matrix
A diagonal matrix used for computing scaling.
NumericVector< Number > * _Re_time
residual vector for time contributions
void solutionInvalidAccumulationTimeStep()
Pass the number of solution invalid occurrences from current iteration to cumulative time iteration c...
virtual void computeJacobianSys(NonlinearImplicitSystem &sys, const NumericVector< Number > &soln, SparseMatrix< Number > &jacobian)
Form a Jacobian matrix.
virtual void solve() override
Solve the system (using libMesh magic)
SolverParams & solverParams()
Get the solver parameters.
NonlinearImplicitSystem::ComputeResidualandJacobian * residual_and_jacobian_object
PetscErrorCode libmesh_petsc_snes_fd_residual(SNES, Vec x, Vec r, void *ctx)
NumericVector< Number > * _Re_non_time
residual vector for non-time contributions
virtual void init() override
Initialize the system.
virtual void computeNullSpace(NonlinearImplicitSystem &sys, std::vector< NumericVector< Number > *> &sp)
void mooseError(Args &&... args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Definition: MooseError.h:284
std::vector< std::pair< R1, R2 > > get(const std::string &param1, const std::string &param2) const
Combine two vector parameters into a single vector of pairs.
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:296
NonlinearImplicitSystem & _nl_implicit_sys
void resetSolutionInvalidTimeStep()
Reset the number of solution invalid occurrences back to zero for the current time step...
void attach_preconditioner(Preconditioner< Number > *preconditioner)
if(subdm)
NonlinearSystem(FEProblemBase &problem, const std::string &name)
virtual void setupFiniteDifferencedPreconditioner() override
ComputeResidualFunctor _nl_residual_functor
virtual void computeBounds(NonlinearImplicitSystem &sys, NumericVector< Number > &lower, NumericVector< Number > &upper)
Solving a linear problem.
Definition: MooseTypes.h:761
void compute_bounds(NumericVector< Number > &lower, NumericVector< Number > &upper, NonlinearImplicitSystem &sys)
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.
const Parallel::Communicator & _communicator
void residualAndJacobianTogether() override
Call this method if you want the residual and Jacobian to be computed simultaneously.
NonlinearImplicitSystem::ComputeResidual * residual_object
void init() override
Initialize the system.
bool hasDampers()
Whether or not this system has dampers.
Specialization of SubProblem for solving nonlinear equations plus auxiliary equations.
void(* jacobian)(const NumericVector< Number > &X, SparseMatrix< Number > &J, sys_type &S)
ierr
virtual NonlinearSolver< Number > * nonlinearSolver() override
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)
Real _resid_vs_jac_scaling_param
The param that indicates the weighting of the residual vs the Jacobian in determining variable scalin...
FEProblemBase & _fe_problem
ComputeFDResidualFunctor _fd_residual_functor
virtual void attachPreconditioner(Preconditioner< Number > *preconditioner) override
Attach a customized preconditioner that requires physics knowledge.
virtual bool shouldUpdateSolution()
Check to see whether the problem should update the solution.
ComputeResidualAndJacobian _resid_and_jac_functor
bool _automatic_scaling
Whether to automatically scale the variables.
Definition: SystemBase.h:979
SolutionInvalidity & solutionInvalidity()
Get the SolutionInvalidity for this app.
Definition: MooseApp.h:149
std::shared_ptr< MoosePreconditioner > _preconditioner
Preconditioner.
Jacobian-Free Newton Krylov.
Definition: MooseTypes.h:758
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(* transpose_nullspace)(std::vector< NumericVector< Number > *> &sp, sys_type &S)
void setupColoringFiniteDifferencedPreconditioner()
According to the nonzero pattern provided in the matrix, a graph is constructed.
bool allowInvalidSolution() const
Whether or not the invalid solutions are allowed.
Moose::SolveType _type
Definition: SolverParams.h:19
virtual void computePostCheck(NonlinearImplicitSystem &sys, const NumericVector< Number > &old_soln, NumericVector< Number > &search_direction, NumericVector< Number > &new_soln, bool &changed_search_direction, bool &changed_new_soln)
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.
const ExecFlagType EXEC_LINEAR
Definition: Moose.C:29
std::string stringify(const T &t)
conversion to string
Definition: Conversion.h:62
void(* nullspace)(std::vector< NumericVector< Number > *> &sp, sys_type &S)
virtual NumericVector< Number > & RHS() override
void(* nearnullspace)(std::vector< NumericVector< Number > *> &sp, sys_type &S)
virtual void close()=0
virtual bool converged() override
Returns the convergence state.
bool solutionInvalid() const
Loop over all the tracked objects and determine whether solution invalid is detected.
virtual SNES getSNES() override
bool _compute_initial_residual_before_preset_bcs
const ExecFlagType EXEC_POSTCHECK
Definition: Moose.C:31
const ExecFlagType EXEC_NONLINEAR
Definition: Moose.C:30
MooseApp & _app
Definition: SystemBase.h:926
Finite difference preconditioner.
Class for containing MooseEnum item information.
Definition: MooseEnumItem.h:18
virtual System & system() override
Get the reference to the libMesh system.
void compute_nullspace(std::vector< NumericVector< Number > *> &sp, NonlinearImplicitSystem &sys)
std::shared_ptr< TimeIntegrator > _time_integrator
Time integrator.
Definition: SystemBase.h:973
const NumericVector< Number > * _current_solution
solution vector from nonlinear solver
void(* bounds)(NumericVector< Number > &XL, NumericVector< Number > &XU, sys_type &S)
bool getFailNextNonlinearConvergenceCheck() const
Whether it will skip further residual evaluations and fail the next nonlinear convergence check...
virtual NonlinearImplicitSystem & sys()
void print(const ConsoleStream &console) const
Print the summary table of Solution Invalid warnings.
const InputParameters & parameters() const
Get the parameters of the object.
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.
virtual void stopSolve(const ExecFlagType &exec_flag) override
Quit the current solve as soon as possible.
virtual void computeResidualSys(NonlinearImplicitSystem &sys, const NumericVector< Number > &soln, NumericVector< Number > &residual)
This function is called by Libmesh to form a residual.
virtual bool hasException()
Whether or not an exception has occurred.
virtual void computeTransposeNullSpace(NonlinearImplicitSystem &sys, std::vector< NumericVector< Number > *> &sp)
virtual void computeNearNullSpace(NonlinearImplicitSystem &sys, std::vector< NumericVector< Number > *> &sp)
virtual ~NonlinearSystem()