libMesh
nonlinear_implicit_system.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 
20 // C++ includes
21 
22 // Local includes
23 #include "libmesh/nonlinear_implicit_system.h"
24 #include "libmesh/diff_solver.h"
25 #include "libmesh/equation_systems.h"
26 #include "libmesh/libmesh_logging.h"
27 #include "libmesh/nonlinear_solver.h"
28 #include "libmesh/sparse_matrix.h"
29 
30 namespace libMesh
31 {
32 
33 // ------------------------------------------------------------
34 // NonlinearImplicitSystem implementation
36  const std::string & name_in,
37  const unsigned int number_in) :
38 
39  Parent (es, name_in, number_in),
40  nonlinear_solver (NonlinearSolver<Number>::build(*this)),
41  diff_solver (),
42  _n_nonlinear_iterations (0),
43  _final_nonlinear_residual (1.e20)
44 {
45  // Set default parameters
46  // These were chosen to match the Petsc defaults
47  es.parameters.set<Real> ("linear solver tolerance") = 1e-5;
48  es.parameters.set<Real> ("linear solver minimum tolerance") = 1e-5;
49  es.parameters.set<unsigned int>("linear solver maximum iterations") = 10000;
50 
51  es.parameters.set<unsigned int>("nonlinear solver maximum iterations") = 50;
52  es.parameters.set<unsigned int>("nonlinear solver maximum function evaluations") = 10000;
53 
54  es.parameters.set<Real>("nonlinear solver absolute residual tolerance") = 1e-35;
55  es.parameters.set<Real>("nonlinear solver relative residual tolerance") = 1e-8;
56  es.parameters.set<Real>("nonlinear solver absolute step tolerance") = 1e-8;
57  es.parameters.set<Real>("nonlinear solver relative step tolerance") = 1e-8;
58 }
59 
60 
61 
63 {
64  // Clear data
65  this->clear();
66 }
67 
68 
69 
71 {
72  // clear the nonlinear solver
73  nonlinear_solver->clear();
74 
75  // clear the parent data
76  Parent::clear();
77 }
78 
79 
80 
82 {
83  // re-initialize the nonlinear solver interface
84  nonlinear_solver->clear();
85 
86  if (diff_solver.get())
87  diff_solver->reinit();
88 
89  // initialize parent data
91 }
92 
93 
94 
96 {
97  // Get a reference to the EquationSystems
98  const EquationSystems & es =
99  this->get_equation_systems();
100 
101  // Get the user-specified nonlinear solver tolerances
102  const unsigned int maxits =
103  es.parameters.get<unsigned int>("nonlinear solver maximum iterations");
104 
105  const unsigned int maxfuncs =
106  es.parameters.get<unsigned int>("nonlinear solver maximum function evaluations");
107 
108  const Real abs_resid_tol =
109  es.parameters.get<Real>("nonlinear solver absolute residual tolerance");
110 
111  const Real rel_resid_tol =
112  es.parameters.get<Real>("nonlinear solver relative residual tolerance");
113 
114  const Real abs_step_tol =
115  es.parameters.get<Real>("nonlinear solver absolute step tolerance");
116 
117  const Real rel_step_tol =
118  es.parameters.get<Real>("nonlinear solver relative step tolerance");
119 
120  // Get the user-specified linear solver tolerances
121  const unsigned int maxlinearits =
122  es.parameters.get<unsigned int>("linear solver maximum iterations");
123 
124  const Real linear_tol =
125  es.parameters.get<Real>("linear solver tolerance");
126 
127  const Real linear_min_tol =
128  es.parameters.get<Real>("linear solver minimum tolerance");
129 
130  // Set all the parameters on the NonlinearSolver
131  nonlinear_solver->max_nonlinear_iterations = maxits;
132  nonlinear_solver->max_function_evaluations = maxfuncs;
133  nonlinear_solver->absolute_residual_tolerance = abs_resid_tol;
134  nonlinear_solver->relative_residual_tolerance = rel_resid_tol;
135  nonlinear_solver->absolute_step_tolerance = abs_step_tol;
136  nonlinear_solver->relative_step_tolerance = rel_step_tol;
137  nonlinear_solver->max_linear_iterations = maxlinearits;
138  nonlinear_solver->initial_linear_tolerance = linear_tol;
139  nonlinear_solver->minimum_linear_tolerance = linear_min_tol;
140 
141  if (diff_solver.get())
142  {
143  diff_solver->max_nonlinear_iterations = maxits;
144  diff_solver->absolute_residual_tolerance = abs_resid_tol;
145  diff_solver->relative_residual_tolerance = rel_resid_tol;
146  diff_solver->absolute_step_tolerance = abs_step_tol;
147  diff_solver->relative_step_tolerance = rel_step_tol;
148  diff_solver->max_linear_iterations = maxlinearits;
149  diff_solver->initial_linear_tolerance = linear_tol;
150  diff_solver->minimum_linear_tolerance = linear_min_tol;
151  }
152 }
153 
154 
155 
157 {
158  // Log how long the nonlinear solve takes.
159  START_LOG("solve()", "System");
160 
161  this->set_solver_parameters();
162 
163  if (diff_solver.get())
164  {
165  diff_solver->solve();
166 
167  // Store the number of nonlinear iterations required to
168  // solve and the final residual.
169  _n_nonlinear_iterations = diff_solver->total_outer_iterations();
170  _final_nonlinear_residual = 0.; // FIXME - support this!
171  }
172  else
173  {
174  if (libMesh::on_command_line("--solver_system_names"))
175  nonlinear_solver->init((this->name()+"_").c_str());
176  else
177  nonlinear_solver->init();
178 
179  // Solve the nonlinear system.
180  const std::pair<unsigned int, Real> rval =
181  nonlinear_solver->solve (*matrix, *solution, *rhs,
182  nonlinear_solver->relative_residual_tolerance,
183  nonlinear_solver->max_linear_iterations);
184 
185  // Store the number of nonlinear iterations required to
186  // solve and the final residual.
187  _n_nonlinear_iterations = rval.first;
188  _final_nonlinear_residual = rval.second;
189  }
190 
191  // Stop logging the nonlinear solve
192  STOP_LOG("solve()", "System");
193 
194  // Update the system after the solve
195  this->update();
196 }
197 
198 
199 
200 std::pair<unsigned int, Real> NonlinearImplicitSystem::get_linear_solve_parameters() const
201 {
202  if (diff_solver.get())
203  return std::make_pair(this->diff_solver->max_linear_iterations,
204  this->diff_solver->relative_residual_tolerance);
205  return std::make_pair(this->nonlinear_solver->max_linear_iterations,
206  this->nonlinear_solver->relative_residual_tolerance);
207 }
208 
209 
210 
211 void NonlinearImplicitSystem::assembly(bool get_residual,
212  bool get_jacobian,
213  bool /*apply_heterogeneous_constraints*/,
214  bool /*apply_no_constraints*/)
215 {
216  // Get current_local_solution in sync
217  this->update();
218 
219  //-----------------------------------------------------------------------------
220  // if the user has provided both function pointers and objects only the pointer
221  // will be used, so catch that as an error
222  if (nonlinear_solver->jacobian && nonlinear_solver->jacobian_object)
223  libmesh_error_msg("ERROR: cannot specify both a function and object to compute the Jacobian!");
224 
225  if (nonlinear_solver->residual && nonlinear_solver->residual_object)
226  libmesh_error_msg("ERROR: cannot specify both a function and object to compute the Residual!");
227 
228  if (nonlinear_solver->matvec && nonlinear_solver->residual_and_jacobian_object)
229  libmesh_error_msg("ERROR: cannot specify both a function and object to compute the combined Residual & Jacobian!");
230 
231 
232  if (get_jacobian)
233  {
234  if (nonlinear_solver->jacobian != libmesh_nullptr)
235  nonlinear_solver->jacobian (*current_local_solution.get(), *matrix, *this);
236 
237  else if (nonlinear_solver->jacobian_object != libmesh_nullptr)
238  nonlinear_solver->jacobian_object->jacobian (*current_local_solution.get(), *matrix, *this);
239 
240  else if (nonlinear_solver->matvec != libmesh_nullptr)
241  nonlinear_solver->matvec (*current_local_solution.get(), get_residual ? rhs : libmesh_nullptr, matrix, *this);
242 
243  else if (nonlinear_solver->residual_and_jacobian_object != libmesh_nullptr)
244  nonlinear_solver->residual_and_jacobian_object->residual_and_jacobian (*current_local_solution.get(), get_residual ? rhs : libmesh_nullptr, matrix, *this);
245 
246  else
247  libmesh_error_msg("Error! Unable to compute residual and/or Jacobian!");
248  }
249 
250  if (get_residual)
251  {
252  if (nonlinear_solver->residual != libmesh_nullptr)
253  nonlinear_solver->residual (*current_local_solution.get(), *rhs, *this);
254 
255  else if (nonlinear_solver->residual_object != libmesh_nullptr)
256  nonlinear_solver->residual_object->residual (*current_local_solution.get(), *rhs, *this);
257 
258  else if (nonlinear_solver->matvec != libmesh_nullptr)
259  {
260  // we might have already grabbed the residual and jacobian together
261  if (!get_jacobian)
263  }
264 
265  else if (nonlinear_solver->residual_and_jacobian_object != libmesh_nullptr)
266  {
267  // we might have already grabbed the residual and jacobian together
268  if (!get_jacobian)
269  nonlinear_solver->residual_and_jacobian_object->residual_and_jacobian (*current_local_solution.get(), rhs, libmesh_nullptr, *this);
270  }
271 
272  else
273  libmesh_error_msg("Error! Unable to compute residual and/or Jacobian!");
274  }
275  else
276  libmesh_assert(get_jacobian); // I can't believe you really wanted to assemble *nothing*
277 }
278 
279 
280 
281 
283 {
284  return nonlinear_solver->get_current_nonlinear_iteration_number();
285 }
286 
287 
288 
289 } // namespace libMesh
void set_solver_parameters()
Copies system parameters into nonlinear solver parameters.
This is the EquationSystems class.
UniquePtr< NumericVector< Number > > current_local_solution
All the values I need to compute my contribution to the simulation at hand.
Definition: system.h:1535
const class libmesh_nullptr_t libmesh_nullptr
NumericVector< Number > * rhs
The system matrix.
This base class can be inherited from to provide interfaces to nonlinear solvers from different packa...
const T & get(const std::string &) const
Definition: parameters.h:430
unsigned get_current_nonlinear_iteration_number() const
If called during the solve(), for example by the user-specified residual or Jacobian function...
NonlinearImplicitSystem(EquationSystems &es, const std::string &name, const unsigned int number)
Constructor.
The libMesh namespace provides an interface to certain functionality in the library.
const std::string & name() const
Definition: system.h:1998
libmesh_assert(j)
UniquePtr< NonlinearSolver< Number > > nonlinear_solver
The NonlinearSolver defines the default interface used to solve the nonlinear_implicit system...
unsigned int _n_nonlinear_iterations
The number of nonlinear iterations required to solve the nonlinear system R(x)=0. ...
This is the base class for classes which contain information related to any physical process that mig...
Definition: system.h:76
Real _final_nonlinear_residual
The final residual for the nonlinear system R(x)
UniquePtr< NumericVector< Number > > solution
Data structure to hold solution values.
Definition: system.h:1523
virtual void solve() libmesh_override
Assembles & solves the nonlinear system R(x) = 0.
virtual std::pair< unsigned int, Real > get_linear_solve_parameters() const libmesh_override
virtual void reinit() libmesh_override
Reinitializes the member data fields associated with the system, so that, e.g., assemble() may be use...
virtual void update()
Update the local values to reflect the solution on neighboring processors.
Definition: system.C:425
virtual void reinit() libmesh_override
Reinitializes the member data fields associated with the system, so that, e.g., assemble() may be use...
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
bool on_command_line(const std::string &arg)
Definition: libmesh.C:921
T & set(const std::string &)
Definition: parameters.h:469
virtual void assembly(bool get_residual, bool get_jacobian, bool apply_heterogeneous_constraints=false, bool apply_no_constraints=false) libmesh_override
Assembles a residual in rhs and/or a jacobian in matrix, as requested.
SparseMatrix< Number > * matrix
The system matrix.
const EquationSystems & get_equation_systems() const
Definition: system.h:712
Parameters parameters
Data structure holding arbitrary parameters.
virtual void clear() libmesh_override
Clear all the data structures associated with the system.
UniquePtr< DiffSolver > diff_solver
The DiffSolver defines an optional interface used to solve the nonlinear_implicit system...
virtual void clear() libmesh_override
Clear all the data structures associated with the system.