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 "TimeIntegrator.h"
15 #include "PetscSupport.h"
16 #include "ComputeResidualFunctor.h"
18 
19 #include "libmesh/nonlinear_solver.h"
20 #include "libmesh/petsc_nonlinear_solver.h"
21 #include "libmesh/sparse_matrix.h"
22 #include "libmesh/petsc_matrix.h"
23 #include "libmesh/default_coupling.h"
24 
25 namespace Moose
26 {
27 void
28 compute_jacobian(const NumericVector<Number> & soln,
29  SparseMatrix<Number> & jacobian,
30  NonlinearImplicitSystem & sys)
31 {
32  FEProblemBase * p =
33  sys.get_equation_systems().parameters.get<FEProblemBase *>("_fe_problem_base");
34  p->computeJacobianSys(sys, soln, jacobian);
35 }
36 
37 void
38 compute_bounds(NumericVector<Number> & lower,
39  NumericVector<Number> & upper,
40  NonlinearImplicitSystem & sys)
41 {
42  FEProblemBase * p =
43  sys.get_equation_systems().parameters.get<FEProblemBase *>("_fe_problem_base");
44  p->computeBounds(sys, lower, upper);
45 }
46 
47 void
48 compute_nullspace(std::vector<NumericVector<Number> *> & sp, NonlinearImplicitSystem & sys)
49 {
50  FEProblemBase * p =
51  sys.get_equation_systems().parameters.get<FEProblemBase *>("_fe_problem_base");
52  p->computeNullSpace(sys, sp);
53 }
54 
55 void
56 compute_transpose_nullspace(std::vector<NumericVector<Number> *> & sp,
57  NonlinearImplicitSystem & sys)
58 {
59  FEProblemBase * p =
60  sys.get_equation_systems().parameters.get<FEProblemBase *>("_fe_problem_base");
61  p->computeTransposeNullSpace(sys, sp);
62 }
63 
64 void
65 compute_nearnullspace(std::vector<NumericVector<Number> *> & sp, NonlinearImplicitSystem & sys)
66 {
67  FEProblemBase * p =
68  sys.get_equation_systems().parameters.get<FEProblemBase *>("_fe_problem_base");
69  p->computeNearNullSpace(sys, sp);
70 }
71 
72 void
73 compute_postcheck(const NumericVector<Number> & old_soln,
74  NumericVector<Number> & search_direction,
75  NumericVector<Number> & new_soln,
76  bool & changed_search_direction,
77  bool & changed_new_soln,
78  NonlinearImplicitSystem & sys)
79 {
80  FEProblemBase * p =
81  sys.get_equation_systems().parameters.get<FEProblemBase *>("_fe_problem_base");
83  sys, old_soln, search_direction, new_soln, changed_search_direction, changed_new_soln);
84 }
85 } // namespace Moose
86 
87 NonlinearSystem::NonlinearSystem(FEProblemBase & fe_problem, const std::string & name)
89  fe_problem, fe_problem.es().add_system<TransientNonlinearImplicitSystem>(name), name),
90  _transient_sys(fe_problem.es().get_system<TransientNonlinearImplicitSystem>(name)),
91  _nl_residual_functor(_fe_problem),
92  _fd_residual_functor(_fe_problem),
93  _use_coloring_finite_difference(false)
94 {
95  _sys.get_dof_map().add_coupling_functor(
96  _sys.get_dof_map().default_coupling(),
97  false); // The false keeps it from getting added to the mesh
98 
99  nonlinearSolver()->residual_object = &_nl_residual_functor;
103  nonlinearSolver()->transpose_nullspace = Moose::compute_transpose_nullspace;
105 
106 #ifdef LIBMESH_HAVE_PETSC
107  PetscNonlinearSolver<Real> * petsc_solver =
108  static_cast<PetscNonlinearSolver<Real> *>(_transient_sys.nonlinear_solver.get());
109  if (petsc_solver)
110  {
111  petsc_solver->set_residual_zero_out(false);
112  petsc_solver->set_jacobian_zero_out(false);
113  petsc_solver->use_default_monitor(false);
114  }
115 #endif
116 }
117 
119 
120 SparseMatrix<Number> &
122 {
123  if (!_subproblem.matrixTagExists(tag))
124  mooseError("Cannot add a tagged matrix with matrix_tag, ",
125  tag,
126  ", that tag does not exist in System ",
127  name());
128 
129  if (hasMatrix(tag))
130  return getMatrix(tag);
131 
132  auto matrix_name = _subproblem.matrixTagName(tag);
133 
134  SparseMatrix<Number> * mat = &_transient_sys.add_matrix(matrix_name);
135 
136  if (_tagged_matrices.size() < tag + 1)
137  _tagged_matrices.resize(tag + 1);
138 
139  _tagged_matrices[tag] = mat;
140 
141  return *mat;
142 }
143 
144 void
146 {
147  // Only attach the postcheck function to the solver if we actually
148  // have dampers or if the FEProblemBase needs to update the solution,
149  // which is also done during the linesearch postcheck. It doesn't
150  // hurt to do this multiple times, it is just setting a pointer.
153  _transient_sys.nonlinear_solver->postcheck = Moose::compute_postcheck;
154 
156  {
157  // Calculate the initial residual for use in the convergence criterion.
161  _transient_sys.rhs->close();
164  _console << "Initial residual before setting preset BCs: "
166  }
167 
168  // Clear the iteration counters
169  _current_l_its.clear();
170  _current_nl_its = 0;
171 
172  // Initialize the solution vector using a predictor and known values from nodal bcs
174 
176  {
177  _transient_sys.nonlinear_solver->fd_residual_object = &_fd_residual_functor;
179  }
180 
181 #ifdef LIBMESH_HAVE_PETSC
182  PetscNonlinearSolver<Real> & solver =
183  static_cast<PetscNonlinearSolver<Real> &>(*_transient_sys.nonlinear_solver);
184  solver.mffd_residual_object = &_fd_residual_functor;
185 
186  solver.set_snesmf_reuse_base(_fe_problem.useSNESMFReuseBase());
187 #endif
188 
189  if (_time_integrator)
190  {
191  _time_integrator->solve();
192  _time_integrator->postSolve();
193  _n_iters = _time_integrator->getNumNonlinearIterations();
194  _n_linear_iters = _time_integrator->getNumLinearIterations();
195  }
196  else
197  {
198  system().solve();
199  _n_iters = _transient_sys.n_nonlinear_iterations();
200 #ifdef LIBMESH_HAVE_PETSC
201  _n_linear_iters = solver.get_total_linear_iterations();
202 #endif
203  }
204 
205  // store info about the solve
206  _final_residual = _transient_sys.final_nonlinear_residual();
207 
208 #ifdef LIBMESH_HAVE_PETSC
210 #if PETSC_VERSION_LESS_THAN(3, 2, 0)
211  MatFDColoringDestroy(_fdcoloring);
212 #else
213  MatFDColoringDestroy(&_fdcoloring);
214 #endif
215 #endif
216 }
217 
218 void
220 {
221 #ifdef LIBMESH_HAVE_PETSC
222 #if PETSC_VERSION_LESS_THAN(3, 0, 0)
223 #else
224  PetscNonlinearSolver<Real> & solver =
225  static_cast<PetscNonlinearSolver<Real> &>(*sys().nonlinear_solver);
226  SNESSetFunctionDomainError(solver.snes());
227 #endif
228 #endif
229 
230  // Insert a NaN into the residual vector. As of PETSc-3.6, this
231  // should make PETSc return DIVERGED_NANORINF the next time it does
232  // a reduction. We'll write to the first local dof on every
233  // processor that has any dofs.
234  _transient_sys.rhs->close();
235 
236  if (_transient_sys.rhs->local_size())
237  _transient_sys.rhs->set(_transient_sys.rhs->first_local_index(),
238  std::numeric_limits<Real>::quiet_NaN());
239  _transient_sys.rhs->close();
240 
241  // Clean up by getting other vectors into a valid state for a
242  // (possible) subsequent solve. There may be more than just
243  // these...
244  if (_Re_time)
245  _Re_time->close();
246  _Re_non_time->close();
247 }
248 
249 void
251 {
252  std::shared_ptr<FiniteDifferencePreconditioner> fdp =
254  if (!fdp)
255  mooseError("Did not setup finite difference preconditioner, and please add a preconditioning "
256  "block with type = fdp");
257 
258  if (fdp->finiteDifferenceType() == "coloring")
259  {
262  }
263 
264  else if (fdp->finiteDifferenceType() == "standard")
265  {
268  }
269  else
270  mooseError("Unknown finite difference type");
271 }
272 
273 void
275 {
276 #if LIBMESH_HAVE_PETSC
277  // Make sure that libMesh isn't going to override our preconditioner
278  _transient_sys.nonlinear_solver->jacobian = nullptr;
279 
280  PetscNonlinearSolver<Number> * petsc_nonlinear_solver =
281  static_cast<PetscNonlinearSolver<Number> *>(_transient_sys.nonlinear_solver.get());
282 
283  PetscMatrix<Number> * petsc_mat = static_cast<PetscMatrix<Number> *>(_transient_sys.matrix);
284 
285  SNESSetJacobian(petsc_nonlinear_solver->snes(),
286  petsc_mat->mat(),
287  petsc_mat->mat(),
288 #if PETSC_VERSION_LESS_THAN(3, 4, 0)
289  SNESDefaultComputeJacobian,
290 #else
291  SNESComputeJacobianDefault,
292 #endif
293  nullptr);
294 #endif
295 }
296 
297 void
299 {
300 #ifdef LIBMESH_HAVE_PETSC
301  // Make sure that libMesh isn't going to override our preconditioner
302  _transient_sys.nonlinear_solver->jacobian = nullptr;
303 
304  PetscNonlinearSolver<Number> & petsc_nonlinear_solver =
305  dynamic_cast<PetscNonlinearSolver<Number> &>(*_transient_sys.nonlinear_solver);
306 
307  // Pointer to underlying PetscMatrix type
308  PetscMatrix<Number> * petsc_mat = dynamic_cast<PetscMatrix<Number> *>(_transient_sys.matrix);
309 
310 #if PETSC_VERSION_LESS_THAN(3, 2, 0)
311  // This variable is only needed for PETSC < 3.2.0
312  PetscVector<Number> * petsc_vec =
313  dynamic_cast<PetscVector<Number> *>(_transient_sys.solution.get());
314 #endif
315 
316  Moose::compute_jacobian(*_transient_sys.current_local_solution, *petsc_mat, _transient_sys);
317 
318  if (!petsc_mat)
319  mooseError("Could not convert to Petsc matrix.");
320 
321  petsc_mat->close();
322 
323  PetscErrorCode ierr = 0;
324  ISColoring iscoloring;
325 
326 #if PETSC_VERSION_LESS_THAN(3, 2, 0)
327  // PETSc 3.2.x
328  ierr = MatGetColoring(petsc_mat->mat(), MATCOLORING_LF, &iscoloring);
329  CHKERRABORT(libMesh::COMM_WORLD, ierr);
330 #elif PETSC_VERSION_LESS_THAN(3, 5, 0)
331  // PETSc 3.3.x, 3.4.x
332  ierr = MatGetColoring(petsc_mat->mat(), MATCOLORINGLF, &iscoloring);
333  CHKERRABORT(_communicator.get(), ierr);
334 #else
335  // PETSc 3.5.x
336  MatColoring matcoloring;
337  ierr = MatColoringCreate(petsc_mat->mat(), &matcoloring);
338  CHKERRABORT(_communicator.get(), ierr);
339  ierr = MatColoringSetType(matcoloring, MATCOLORINGLF);
340  CHKERRABORT(_communicator.get(), ierr);
341  ierr = MatColoringSetFromOptions(matcoloring);
342  CHKERRABORT(_communicator.get(), ierr);
343  ierr = MatColoringApply(matcoloring, &iscoloring);
344  CHKERRABORT(_communicator.get(), ierr);
345  ierr = MatColoringDestroy(&matcoloring);
346  CHKERRABORT(_communicator.get(), ierr);
347 #endif
348 
349  MatFDColoringCreate(petsc_mat->mat(), iscoloring, &_fdcoloring);
350  MatFDColoringSetFromOptions(_fdcoloring);
351  // clang-format off
352  MatFDColoringSetFunction(_fdcoloring,
353  (PetscErrorCode(*)(void))(void (*)(void)) &
354  libMesh::libmesh_petsc_snes_fd_residual,
355  &petsc_nonlinear_solver);
356  // clang-format on
357 #if !PETSC_RELEASE_LESS_THAN(3, 5, 0)
358  MatFDColoringSetUp(petsc_mat->mat(), iscoloring, _fdcoloring);
359 #endif
360 #if PETSC_VERSION_LESS_THAN(3, 4, 0)
361  SNESSetJacobian(petsc_nonlinear_solver.snes(),
362  petsc_mat->mat(),
363  petsc_mat->mat(),
364  SNESDefaultComputeJacobianColor,
365  _fdcoloring);
366 #else
367  SNESSetJacobian(petsc_nonlinear_solver.snes(),
368  petsc_mat->mat(),
369  petsc_mat->mat(),
370  SNESComputeJacobianDefaultColor,
371  _fdcoloring);
372 #endif
373 #if PETSC_VERSION_LESS_THAN(3, 2, 0)
374  Mat my_mat = petsc_mat->mat();
375  MatStructure my_struct;
376 
377  SNESComputeJacobian(
378  petsc_nonlinear_solver.snes(), petsc_vec->vec(), &my_mat, &my_mat, &my_struct);
379 #endif
380 
381 #if PETSC_VERSION_LESS_THAN(3, 2, 0)
382  ISColoringDestroy(iscoloring);
383 #else
384  // PETSc 3.3.0
385  ISColoringDestroy(&iscoloring);
386 #endif
387 
388 #endif
389 }
390 
391 bool
393 {
395  return false;
396 
397  return _transient_sys.nonlinear_solver->converged;
398 }
void compute_jacobian(const NumericVector< Number > &soln, SparseMatrix< Number > &jacobian, NonlinearImplicitSystem &sys)
virtual bool hasMatrix(TagID tag) const
Check if the tagged matrix exists in the system.
Definition: SystemBase.C:805
virtual SparseMatrix< Number > & addMatrix(TagID tag) override
Adds a jacobian sized vector.
NumericVector< Number > * _Re_time
residual vector for time contributions
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.
NumericVector< Number > * _Re_non_time
residual vector for non-time contributions
unsigned int TagID
Definition: MooseTypes.h:162
if(nl->nonlinearSolver() ->matvec &&nl->nonlinearSolver() ->residual_and_jacobian_object)
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:207
void compute_nearnullspace(std::vector< NumericVector< Number > *> &sp, NonlinearImplicitSystem &sys)
std::vector< SparseMatrix< Number > * > _tagged_matrices
Tagged matrices (pointer)
Definition: SystemBase.h:758
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:595
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.
virtual void stopSolve() override
Quit the current solve as soon as possible.
bool hasDampers()
Whether or not this system has dampers.
Specialization of SubProblem for solving nonlinear equations plus auxiliary equations.
virtual NonlinearSolver< Number > * nonlinearSolver() override
virtual bool matrixTagExists(const TagName &tag_name)
Check to see if a particular Tag exists.
Definition: SubProblem.C:129
Nonlinear system to be solved.
virtual const std::string & name() const
Definition: SystemBase.C:1088
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)
FEProblemBase & _fe_problem
std::vector< unsigned int > _current_l_its
ComputeFDResidualFunctor _fd_residual_functor
virtual bool shouldUpdateSolution()
Check to see whether the problem should update the solution.
const InputParameters & parameters() const
Get the parameters of the object.
Definition: MooseObject.h:65
std::shared_ptr< MoosePreconditioner > _preconditioner
Preconditioner.
SubProblem & _subproblem
Definition: SystemBase.h:735
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.
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
TransientNonlinearImplicitSystem & _transient_sys
virtual bool converged() override
Returns the convergence state.
bool _compute_initial_residual_before_preset_bcs
Finite difference preconditioner.
virtual System & system() override
Get the reference to the libMesh system.
virtual SparseMatrix< Number > & getMatrix(TagID tag)
Get a raw SparseMatrix.
Definition: SystemBase.C:811
void compute_nullspace(std::vector< NumericVector< Number > *> &sp, NonlinearImplicitSystem &sys)
virtual TransientNonlinearImplicitSystem & sys()
std::shared_ptr< TimeIntegrator > _time_integrator
Time integrator.
Definition: SystemBase.h:782
const NumericVector< Number > * _current_solution
solution vector from nonlinear solver
void setupStandardFiniteDifferencedPreconditioner()
Form preconditioning matrix via a standard finite difference method column-by-column.
bool useSNESMFReuseBase()
Return a flag that indicates if we are reusing the vector base.
ierr
Definition: Moose.h:112
const ConsoleStream _console
An instance of helper class to write streams to the Console objects.
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 TagName matrixTagName(TagID tag)
Retrieve the name associated with a TagID.
Definition: SubProblem.C:158
virtual ~NonlinearSystem()