nonlinear_implicit_system.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2014 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
24 #include "libmesh/diff_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 (NULL),
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-specifiied 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 toleranaces
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  // Solve the nonlinear system.
175  const std::pair<unsigned int, Real> rval =
176  nonlinear_solver->solve (*matrix, *solution, *rhs,
177  nonlinear_solver->relative_residual_tolerance,
178  nonlinear_solver->max_linear_iterations);
179 
180  // Store the number of nonlinear iterations required to
181  // solve and the final residual.
182  _n_nonlinear_iterations = rval.first;
183  _final_nonlinear_residual = rval.second;
184  }
185 
186  // Stop logging the nonlinear solve
187  STOP_LOG("solve()", "System");
188 
189  // Update the system after the solve
190  this->update();
191 }
192 
193 
194 
195 std::pair<unsigned int, Real> NonlinearImplicitSystem::get_linear_solve_parameters() const
196 {
197  if (diff_solver.get())
198  return std::make_pair(this->diff_solver->max_linear_iterations,
199  this->diff_solver->relative_residual_tolerance);
200  return std::make_pair(this->nonlinear_solver->max_linear_iterations,
201  this->nonlinear_solver->relative_residual_tolerance);
202 }
203 
204 
205 
206 void NonlinearImplicitSystem::assembly(bool get_residual,
207  bool get_jacobian)
208 {
209  // Get current_local_solution in sync
210  this->update();
211 
212  //-----------------------------------------------------------------------------
213  // if the user has provided both function pointers and objects only the pointer
214  // will be used, so catch that as an error
215  if (nonlinear_solver->jacobian && nonlinear_solver->jacobian_object)
216  {
217  libMesh::err << "ERROR: cannot specifiy both a function and object to compute the Jacobian!" << std::endl;
218  libmesh_error();
219  }
220 
221  if (nonlinear_solver->residual && nonlinear_solver->residual_object)
222  {
223  libMesh::err << "ERROR: cannot specifiy both a function and object to compute the Residual!" << std::endl;
224  libmesh_error();
225  }
226 
227  if (nonlinear_solver->matvec && nonlinear_solver->residual_and_jacobian_object)
228  {
229  libMesh::err << "ERROR: cannot specifiy both a function and object to compute the combined Residual & Jacobian!" << std::endl;
230  libmesh_error();
231  }
232  //-----------------------------------------------------------------------------
233 
234 
235  if (get_jacobian)
236  {
237  if (nonlinear_solver->jacobian != NULL)
238  nonlinear_solver->jacobian (*current_local_solution.get(), *matrix, *this);
239 
240  else if (nonlinear_solver->jacobian_object != NULL)
241  nonlinear_solver->jacobian_object->jacobian (*current_local_solution.get(), *matrix, *this);
242 
243  else if (nonlinear_solver->matvec != NULL)
244  nonlinear_solver->matvec (*current_local_solution.get(), get_residual?rhs:NULL, matrix, *this);
245 
246  else if (nonlinear_solver->residual_and_jacobian_object != NULL)
247  nonlinear_solver->residual_and_jacobian_object->residual_and_jacobian (*current_local_solution.get(), get_residual?rhs:NULL, matrix, *this);
248 
249  else
250  libmesh_error();
251  }
252 
253  if (get_residual)
254  {
255  if (nonlinear_solver->residual != NULL)
256  nonlinear_solver->residual (*current_local_solution.get(), *rhs, *this);
257 
258  else if (nonlinear_solver->residual_object != NULL)
259  nonlinear_solver->residual_object->residual (*current_local_solution.get(), *rhs, *this);
260 
261  else if (nonlinear_solver->matvec != NULL)
262  {
263  // we might have already grabbed the residual and jacobian together
264  if (!get_jacobian)
265  nonlinear_solver->matvec (*current_local_solution.get(), rhs, NULL, *this);
266  }
267 
268  else if (nonlinear_solver->residual_and_jacobian_object != NULL)
269  {
270  // we might have already grabbed the residual and jacobian together
271  if (!get_jacobian)
272  nonlinear_solver->residual_and_jacobian_object->residual_and_jacobian (*current_local_solution.get(), rhs, NULL, *this);
273  }
274 
275  else
276  libmesh_error();
277  }
278  else
279  libmesh_assert(get_jacobian); // I can't believe you really wanted to assemble *nothing*
280 }
281 
282 
283 
284 
286 {
287  return nonlinear_solver->get_current_nonlinear_iteration_number();
288 }
289 
290 
291 
292 } // namespace libMesh

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:06 UTC

Hosted By:
SourceForge.net Logo