www.mooseframework.org
NonlinearSystem.C
Go to the documentation of this file.
1 /****************************************************************/
2 /* DO NOT MODIFY THIS HEADER */
3 /* MOOSE - Multiphysics Object Oriented Simulation Environment */
4 /* */
5 /* (c) 2010 Battelle Energy Alliance, LLC */
6 /* ALL RIGHTS RESERVED */
7 /* */
8 /* Prepared by Battelle Energy Alliance, LLC */
9 /* Under Contract No. DE-AC07-05ID14517 */
10 /* With the U. S. Department of Energy */
11 /* */
12 /* See COPYRIGHT for full restrictions */
13 /****************************************************************/
14 
15 // moose includes
16 #include "NonlinearSystem.h"
17 #include "FEProblem.h"
18 #include "TimeIntegrator.h"
20 
21 #include "libmesh/nonlinear_solver.h"
22 #include "libmesh/petsc_nonlinear_solver.h"
23 #include "libmesh/sparse_matrix.h"
24 #include "libmesh/petsc_matrix.h"
25 
26 namespace Moose
27 {
28 void
29 compute_jacobian(const NumericVector<Number> & soln,
30  SparseMatrix<Number> & jacobian,
31  NonlinearImplicitSystem & sys)
32 {
33  FEProblemBase * p =
34  sys.get_equation_systems().parameters.get<FEProblemBase *>("_fe_problem_base");
35  p->computeJacobian(sys, soln, jacobian);
36 }
37 
38 void
39 compute_residual(const NumericVector<Number> & soln,
40  NumericVector<Number> & residual,
41  NonlinearImplicitSystem & sys)
42 {
43  FEProblemBase * p =
44  sys.get_equation_systems().parameters.get<FEProblemBase *>("_fe_problem_base");
45  p->computeResidual(sys, soln, residual);
46 }
47 
48 void
49 compute_bounds(NumericVector<Number> & lower,
50  NumericVector<Number> & upper,
51  NonlinearImplicitSystem & sys)
52 {
53  FEProblemBase * p =
54  sys.get_equation_systems().parameters.get<FEProblemBase *>("_fe_problem_base");
55  p->computeBounds(sys, lower, upper);
56 }
57 
58 void
59 compute_nullspace(std::vector<NumericVector<Number> *> & sp, NonlinearImplicitSystem & sys)
60 {
61  FEProblemBase * p =
62  sys.get_equation_systems().parameters.get<FEProblemBase *>("_fe_problem_base");
63  p->computeNullSpace(sys, sp);
64 }
65 
66 void
67 compute_transpose_nullspace(std::vector<NumericVector<Number> *> & sp,
68  NonlinearImplicitSystem & sys)
69 {
70  FEProblemBase * p =
71  sys.get_equation_systems().parameters.get<FEProblemBase *>("_fe_problem_base");
72  p->computeTransposeNullSpace(sys, sp);
73 }
74 
75 void
76 compute_nearnullspace(std::vector<NumericVector<Number> *> & sp, NonlinearImplicitSystem & sys)
77 {
78  FEProblemBase * p =
79  sys.get_equation_systems().parameters.get<FEProblemBase *>("_fe_problem_base");
80  p->computeNearNullSpace(sys, sp);
81 }
82 
83 void
84 compute_postcheck(const NumericVector<Number> & old_soln,
85  NumericVector<Number> & search_direction,
86  NumericVector<Number> & new_soln,
87  bool & changed_search_direction,
88  bool & changed_new_soln,
89  NonlinearImplicitSystem & sys)
90 {
91  FEProblemBase * p =
92  sys.get_equation_systems().parameters.get<FEProblemBase *>("_fe_problem_base");
94  sys, old_soln, search_direction, new_soln, changed_search_direction, changed_new_soln);
95 }
96 } // namespace Moose
97 
98 NonlinearSystem::NonlinearSystem(FEProblemBase & fe_problem, const std::string & name)
100  fe_problem, fe_problem.es().add_system<TransientNonlinearImplicitSystem>(name), name),
101  _transient_sys(fe_problem.es().get_system<TransientNonlinearImplicitSystem>(name)),
102  _use_coloring_finite_difference(false)
103 {
108  nonlinearSolver()->transpose_nullspace = Moose::compute_transpose_nullspace;
110 
111 #ifdef LIBMESH_HAVE_PETSC
112  PetscNonlinearSolver<Real> * petsc_solver =
113  static_cast<PetscNonlinearSolver<Real> *>(_transient_sys.nonlinear_solver.get());
114  if (petsc_solver)
115  {
116  petsc_solver->set_residual_zero_out(false);
117  petsc_solver->set_jacobian_zero_out(false);
118  petsc_solver->use_default_monitor(false);
119  }
120 #endif
121 }
122 
124 
125 void
127 {
128  // Only attach the postcheck function to the solver if we actually
129  // have dampers or if the FEProblemBase needs to update the solution,
130  // which is also done during the linesearch postcheck. It doesn't
131  // hurt to do this multiple times, it is just setting a pointer.
134  _transient_sys.nonlinear_solver->postcheck = Moose::compute_postcheck;
135 
137  {
138  // Calculate the initial residual for use in the convergence criterion.
142  _transient_sys.rhs->close();
145  _console << "Initial residual before setting preset BCs: "
147  }
148 
149  // Clear the iteration counters
150  _current_l_its.clear();
151  _current_nl_its = 0;
152 
153  // Initialize the solution vector using a predictor and known values from nodal bcs
155 
158 
159  if (_time_integrator)
160  {
161  _time_integrator->solve();
162  _time_integrator->postSolve();
163  }
164  else
165  system().solve();
166 
167  // store info about the solve
168  _n_iters = _transient_sys.n_nonlinear_iterations();
169  _final_residual = _transient_sys.final_nonlinear_residual();
170 
171 #ifdef LIBMESH_HAVE_PETSC
172  _n_linear_iters = static_cast<PetscNonlinearSolver<Real> &>(*_transient_sys.nonlinear_solver)
173  .get_total_linear_iterations();
174 #endif
175 
176 #ifdef LIBMESH_HAVE_PETSC
178 #if PETSC_VERSION_LESS_THAN(3, 2, 0)
179  MatFDColoringDestroy(_fdcoloring);
180 #else
181  MatFDColoringDestroy(&_fdcoloring);
182 #endif
183 #endif
184 }
185 
186 void
188 {
189 #ifdef LIBMESH_HAVE_PETSC
190 #if PETSC_VERSION_LESS_THAN(3, 0, 0)
191 #else
192  PetscNonlinearSolver<Real> & solver =
193  static_cast<PetscNonlinearSolver<Real> &>(*sys().nonlinear_solver);
194  SNESSetFunctionDomainError(solver.snes());
195 #endif
196 #endif
197 
198  // Insert a NaN into the residual vector. As of PETSc-3.6, this
199  // should make PETSc return DIVERGED_NANORINF the next time it does
200  // a reduction. We'll write to the first local dof on every
201  // processor that has any dofs.
202  if (_transient_sys.rhs->local_size())
203  _transient_sys.rhs->set(_transient_sys.rhs->first_local_index(),
204  std::numeric_limits<Real>::quiet_NaN());
205  _transient_sys.rhs->close();
206 
207  // Clean up by getting other vectors into a valid state for a
208  // (possible) subsequent solve. There may be more than just
209  // these...
210  if (_Re_time)
211  _Re_time->close();
212  _Re_non_time->close();
213 }
214 
215 void
217 {
218  std::shared_ptr<FiniteDifferencePreconditioner> fdp =
219  std::dynamic_pointer_cast<FiniteDifferencePreconditioner>(_preconditioner);
220  if (!fdp)
221  mooseError("Did not setup finite difference preconditioner, and please add a preconditioning "
222  "block with type = fdp");
223 
224  if (fdp->finiteDifferenceType() == "coloring")
225  {
228  }
229 
230  else if (fdp->finiteDifferenceType() == "standard")
231  {
234  }
235  else
236  mooseError("Unknown finite difference type");
237 }
238 
239 void
241 {
242 #if LIBMESH_HAVE_PETSC
243  PetscNonlinearSolver<Number> * petsc_nonlinear_solver =
244  static_cast<PetscNonlinearSolver<Number> *>(_transient_sys.nonlinear_solver.get());
245 
246  PetscMatrix<Number> * petsc_mat = static_cast<PetscMatrix<Number> *>(_transient_sys.matrix);
247 
248  SNESSetJacobian(petsc_nonlinear_solver->snes(),
249  petsc_mat->mat(),
250  petsc_mat->mat(),
251 #if PETSC_VERSION_LESS_THAN(3, 4, 0)
252  SNESDefaultComputeJacobian,
253 #else
254  SNESComputeJacobianDefault,
255 #endif
256  nullptr);
257 #endif
258 }
259 
260 void
262 {
263 #ifdef LIBMESH_HAVE_PETSC
264  // Make sure that libMesh isn't going to override our preconditioner
265  _transient_sys.nonlinear_solver->jacobian = NULL;
266 
267  PetscNonlinearSolver<Number> & petsc_nonlinear_solver =
268  dynamic_cast<PetscNonlinearSolver<Number> &>(*_transient_sys.nonlinear_solver);
269 
270  // Pointer to underlying PetscMatrix type
271  PetscMatrix<Number> * petsc_mat = dynamic_cast<PetscMatrix<Number> *>(_transient_sys.matrix);
272 
273 #if PETSC_VERSION_LESS_THAN(3, 2, 0)
274  // This variable is only needed for PETSC < 3.2.0
275  PetscVector<Number> * petsc_vec =
276  dynamic_cast<PetscVector<Number> *>(_transient_sys.solution.get());
277 #endif
278 
279  Moose::compute_jacobian(*_transient_sys.current_local_solution, *petsc_mat, _transient_sys);
280 
281  if (!petsc_mat)
282  mooseError("Could not convert to Petsc matrix.");
283 
284  petsc_mat->close();
285 
286  PetscErrorCode ierr = 0;
287  ISColoring iscoloring;
288 
289 #if PETSC_VERSION_LESS_THAN(3, 2, 0)
290  // PETSc 3.2.x
291  ierr = MatGetColoring(petsc_mat->mat(), MATCOLORING_LF, &iscoloring);
292  CHKERRABORT(libMesh::COMM_WORLD, ierr);
293 #elif PETSC_VERSION_LESS_THAN(3, 5, 0)
294  // PETSc 3.3.x, 3.4.x
295  ierr = MatGetColoring(petsc_mat->mat(), MATCOLORINGLF, &iscoloring);
296  CHKERRABORT(_communicator.get(), ierr);
297 #else
298  // PETSc 3.5.x
299  MatColoring matcoloring;
300  ierr = MatColoringCreate(petsc_mat->mat(), &matcoloring);
301  CHKERRABORT(_communicator.get(), ierr);
302  ierr = MatColoringSetType(matcoloring, MATCOLORINGLF);
303  CHKERRABORT(_communicator.get(), ierr);
304  ierr = MatColoringSetFromOptions(matcoloring);
305  CHKERRABORT(_communicator.get(), ierr);
306  ierr = MatColoringApply(matcoloring, &iscoloring);
307  CHKERRABORT(_communicator.get(), ierr);
308  ierr = MatColoringDestroy(&matcoloring);
309  CHKERRABORT(_communicator.get(), ierr);
310 #endif
311 
312  MatFDColoringCreate(petsc_mat->mat(), iscoloring, &_fdcoloring);
313  MatFDColoringSetFromOptions(_fdcoloring);
314  MatFDColoringSetFunction(_fdcoloring,
315  (PetscErrorCode(*)(void)) & libMesh::__libmesh_petsc_snes_residual,
316  &petsc_nonlinear_solver);
317 #if !PETSC_RELEASE_LESS_THAN(3, 5, 0)
318  MatFDColoringSetUp(petsc_mat->mat(), iscoloring, _fdcoloring);
319 #endif
320 #if PETSC_VERSION_LESS_THAN(3, 4, 0)
321  SNESSetJacobian(petsc_nonlinear_solver.snes(),
322  petsc_mat->mat(),
323  petsc_mat->mat(),
324  SNESDefaultComputeJacobianColor,
325  _fdcoloring);
326 #else
327  SNESSetJacobian(petsc_nonlinear_solver.snes(),
328  petsc_mat->mat(),
329  petsc_mat->mat(),
330  SNESComputeJacobianDefaultColor,
331  _fdcoloring);
332 #endif
333 #if PETSC_VERSION_LESS_THAN(3, 2, 0)
334  Mat my_mat = petsc_mat->mat();
335  MatStructure my_struct;
336 
337  SNESComputeJacobian(
338  petsc_nonlinear_solver.snes(), petsc_vec->vec(), &my_mat, &my_mat, &my_struct);
339 #endif
340 
341 #if PETSC_VERSION_LESS_THAN(3, 2, 0)
342  ISColoringDestroy(iscoloring);
343 #else
344  // PETSc 3.3.0
345  ISColoringDestroy(&iscoloring);
346 #endif
347 
348 #endif
349 }
350 
351 bool
353 {
355  return false;
356 
357  return _transient_sys.nonlinear_solver->converged;
358 }
void compute_jacobian(const NumericVector< Number > &soln, SparseMatrix< Number > &jacobian, NonlinearImplicitSystem &sys)
NumericVector< Number > * _Re_time
residual vector for time contributions
virtual void computeNullSpace(NonlinearImplicitSystem &sys, std::vector< NumericVector< Number > * > &sp)
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
if(nl->nonlinearSolver() ->matvec &&nl->nonlinearSolver() ->residual_and_jacobian_object)
void mooseError(Args &&...args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Definition: MooseError.h:182
NonlinearSystem(FEProblemBase &problem, const std::string &name)
void compute_transpose_nullspace(std::vector< NumericVector< Number > * > &sp, NonlinearImplicitSystem &sys)
virtual void setupFiniteDifferencedPreconditioner() override
virtual void computeBounds(NonlinearImplicitSystem &sys, NumericVector< Number > &lower, NumericVector< Number > &upper)
Solving a linear problem.
Definition: MooseTypes.h:251
void compute_bounds(NumericVector< Number > &lower, NumericVector< Number > &upper, NonlinearImplicitSystem &sys)
virtual void stopSolve() override
Quit the current solve as soon as possible.
std::shared_ptr< TimeIntegrator > _time_integrator
Time integrator.
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 void computeTransposeNullSpace(NonlinearImplicitSystem &sys, std::vector< NumericVector< Number > * > &sp)
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)
FEProblemBase & _fe_problem
std::vector< unsigned int > _current_l_its
virtual bool shouldUpdateSolution()
Check to see whether the problem should update the solution.
std::shared_ptr< MoosePreconditioner > _preconditioner
Preconditioner.
void compute_residual(const NumericVector< Number > &soln, NumericVector< Number > &residual, NonlinearImplicitSystem &sys)
virtual void computeNearNullSpace(NonlinearImplicitSystem &sys, std::vector< NumericVector< Number > * > &sp)
bool _use_finite_differenced_preconditioner
Whether or not to use a finite differenced preconditioner.
void compute_nullspace(std::vector< NumericVector< Number > * > &sp, NonlinearImplicitSystem &sys)
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:25
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)
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.
void compute_nearnullspace(std::vector< NumericVector< Number > * > &sp, NonlinearImplicitSystem &sys)
virtual System & system() override
Get the reference to the libMesh system.
virtual TransientNonlinearImplicitSystem & sys()
const NumericVector< Number > * _current_solution
solution vector from nonlinear solver
const InputParameters & parameters() const
Get the parameters of the object.
Definition: MooseObject.h:53
void setupStandardFiniteDifferencedPreconditioner()
Form preconditioning matrix via a standard finite difference method column-by-column.
ierr
Definition: Moose.h:84
const ConsoleStream _console
An instance of helper class to write streams to the Console objects.
virtual void computeResidual(NonlinearImplicitSystem &sys, const NumericVector< Number > &soln, NumericVector< Number > &residual)
virtual void computeJacobian(NonlinearImplicitSystem &sys, const NumericVector< Number > &soln, SparseMatrix< Number > &jacobian)
virtual bool hasException()
Whether or not an exception has occurred.
virtual ~NonlinearSystem()