www.mooseframework.org
SingleVariableReturnMappingSolution.C
Go to the documentation of this file.
1 /****************************************************************/
2 /* MOOSE - Multiphysics Object Oriented Simulation Environment */
3 /* */
4 /* All contents are licensed under LGPL V2.1 */
5 /* See LICENSE for full restrictions */
6 /****************************************************************/
8 
9 #include "InputParameters.h"
10 #include <cmath>
11 #include "Conversion.h"
12 
13 template <>
14 InputParameters
16 {
17  InputParameters params = emptyInputParameters();
18 
19  // Newton iteration control parameters
20  params.addParam<unsigned int>("max_its", 30, "Maximum number of Newton iterations");
21  params.addParam<bool>(
22  "output_iteration_info", false, "Set true to output Newton iteration information");
23  params.addDeprecatedParam<bool>(
24  "output_iteration_info_on_error",
25  false,
26  "Set true to output Newton iteration information when those iterations fail",
27  "This information is always output when iterations fail");
28  params.addParam<Real>(
29  "relative_tolerance", 1e-8, "Relative convergence tolerance for Newton iteration");
30  params.addParam<Real>(
31  "absolute_tolerance", 1e-11, "Absolute convergence tolerance for Newton iteration");
32  params.addParam<Real>("acceptable_multiplier",
33  10,
34  "Factor applied to relative and absolute "
35  "tolerance for acceptable convergence if "
36  "iterations are no longer making progress");
37  params.addParam<bool>("legacy_return_mapping",
38  false,
39  "Perform iterations and compute residual "
40  "the same way as the previous "
41  "algorithm. Also use same old defaults for relative_tolerance, "
42  "absolute_tolerance, and max_its.");
43 
44  return params;
45 }
46 
48  const InputParameters & parameters)
49  : _legacy_return_mapping(false),
50  _check_range(false),
51  _max_its(parameters.get<unsigned int>("max_its")),
52  _fixed_max_its(1000), // Far larger than ever expected to be needed
53  _output_iteration_info(parameters.get<bool>("output_iteration_info")),
54  _relative_tolerance(parameters.get<Real>("relative_tolerance")),
55  _absolute_tolerance(parameters.get<Real>("absolute_tolerance")),
56  _acceptable_multiplier(parameters.get<Real>("acceptable_multiplier")),
57  _line_search(true),
58  _bracket_solution(true),
59  _num_resids(30),
60  _residual_history(_num_resids, std::numeric_limits<Real>::max())
61 {
62  if (parameters.get<bool>("legacy_return_mapping") == true)
63  {
64  if (!parameters.isParamSetByUser("relative_tolerance"))
65  _relative_tolerance = 1.e-5;
66  if (!parameters.isParamSetByUser("absolute_tolerance"))
67  _absolute_tolerance = 1.e-20;
68  if (!parameters.isParamSetByUser("max_its"))
69  _max_its = 30;
70  _line_search = false;
71  _bracket_solution = false;
72  _check_range = false;
74  }
75  else
76  {
77  if (parameters.isParamSetByUser("max_its"))
78  mooseWarning("Please remove the parameter 'max_its', as it is no longer used in the return "
79  "mapping procedure.");
80  }
81 }
82 
83 Real
85  const Real /*effective_trial_stress*/) const
86 {
87  return std::numeric_limits<Real>::max();
88 }
89 
90 void
92  Real & scalar,
93  const ConsoleStream & console)
94 {
95  std::stringstream iter_output;
96  std::stringstream * iter_output_ptr = (_output_iteration_info ? &iter_output : nullptr);
97 
99  {
100  if (!internalSolve(effective_trial_stress, scalar, iter_output_ptr))
101  {
102  if (iter_output_ptr)
103  mooseError(iter_output_ptr->str());
104  else
105  {
106  internalSolve(effective_trial_stress, scalar, &iter_output);
107  mooseError(iter_output.str());
108  }
109  }
110  else if (iter_output_ptr)
111  console << iter_output_ptr->str();
112  }
113  else
114  {
115  if (!internalSolveLegacy(effective_trial_stress, scalar, iter_output_ptr))
116  {
117  if (iter_output_ptr)
118  mooseError(iter_output_ptr->str());
119  else
120  {
121  internalSolveLegacy(effective_trial_stress, scalar, &iter_output);
122  mooseError(iter_output.str());
123  }
124  }
125  else if (iter_output_ptr)
126  console << iter_output_ptr->str();
127  }
128 }
129 
130 bool
131 SingleVariableReturnMappingSolution::internalSolve(const Real effective_trial_stress,
132  Real & scalar,
133  std::stringstream * iter_output)
134 {
135  scalar = 0.0;
136  Real scalar_old = 0.0;
137  Real scalar_increment = 0.0;
138  const Real min_permissible_scalar = 0.0;
139  const Real max_permissible_scalar = maximumPermissibleValue(effective_trial_stress);
140  Real scalar_upper_bound = max_permissible_scalar;
141  Real scalar_lower_bound = 0.0;
142  unsigned int it = 0;
143 
144  if (effective_trial_stress == 0.0)
145  {
146  outputIterInfo(iter_output, it, effective_trial_stress, scalar, 0.0, 1.0);
147  return true;
148  }
149 
150  Real residual = computeResidual(effective_trial_stress, scalar);
151  Real residual_old = residual;
152 
153  Real init_resid_norm = std::abs(residual);
154  if (init_resid_norm == 0.0)
155  init_resid_norm = 1.0;
156  Real init_resid_sign = (residual < 0.0 ? -1.0 : 1.0);
157 
158  Real reference_residual = computeReferenceResidual(effective_trial_stress, scalar);
159 
160  if (converged(residual, reference_residual))
161  {
162  iterationFinalize(scalar);
163  outputIterInfo(iter_output, it, effective_trial_stress, scalar, residual, reference_residual);
164  return true;
165  }
166 
167  _residual_history.assign(_num_resids, std::numeric_limits<Real>::max());
168  _residual_history[0] = residual;
169 
170  while (it < _fixed_max_its && !converged(residual, reference_residual) &&
171  !convergedAcceptable(it, residual, reference_residual))
172  {
173  scalar_increment = -residual / computeDerivative(effective_trial_stress, scalar);
174  scalar = scalar_old + scalar_increment;
175 
176  if (_check_range)
177  checkPermissibleRange(scalar,
178  scalar_increment,
179  scalar_old,
180  min_permissible_scalar,
181  max_permissible_scalar,
182  iter_output);
183 
184  residual = computeResidual(effective_trial_stress, scalar);
185  reference_residual = computeReferenceResidual(effective_trial_stress, scalar);
186  iterationFinalize(scalar);
187 
188  if (_bracket_solution)
189  updateBounds(
190  scalar, residual, init_resid_sign, scalar_upper_bound, scalar_lower_bound, iter_output);
191 
192  if (converged(residual, reference_residual))
193  {
194  outputIterInfo(iter_output, it, effective_trial_stress, scalar, residual, reference_residual);
195  break;
196  }
197  else
198  {
199  bool modified_increment = false;
200 
201  // Line Search
202  if (_line_search)
203  {
204  if (residual_old - residual != 0.0)
205  {
206  Real alpha = residual_old / (residual_old - residual);
207  if (alpha > 1.0) // upper bound for alpha
208  alpha = 1.0;
209  else if (alpha < 1e-2) // lower bound for alpha
210  alpha = 1e-2;
211  if (alpha != 1.0)
212  {
213  modified_increment = true;
214  scalar_increment *= alpha;
215  if (iter_output)
216  *iter_output << " Line search alpha = " << alpha
217  << " increment = " << scalar_increment << std::endl;
218  }
219  }
220  }
221 
222  if (_bracket_solution)
223  {
224  // Check to see whether trial scalar_increment is outside the bounds, and set it to a point
225  // within the bounds if it is
226  if (scalar_old + scalar_increment >= scalar_upper_bound ||
227  scalar_old + scalar_increment <= scalar_lower_bound)
228  {
229  if (scalar_upper_bound != max_permissible_scalar && scalar_lower_bound != 0.0)
230  {
231  Real frac = 0.5;
232  scalar_increment =
233  (1.0 - frac) * scalar_lower_bound + frac * scalar_upper_bound - scalar_old;
234  modified_increment = true;
235  if (iter_output)
236  *iter_output << " Trial scalar_increment exceeded bounds. Setting between "
237  "lower/upper bounds. frac: "
238  << frac << std::endl;
239  }
240  }
241  }
242 
243  // Update the trial scalar and recompute residual if the line search or bounds checking
244  // modified the increment
245  if (modified_increment)
246  {
247  scalar = scalar_old + scalar_increment;
248  residual = computeResidual(effective_trial_stress, scalar);
249  reference_residual = computeReferenceResidual(effective_trial_stress, scalar);
250  iterationFinalize(scalar);
251 
252  if (_bracket_solution)
253  updateBounds(scalar,
254  residual,
255  init_resid_sign,
256  scalar_upper_bound,
257  scalar_lower_bound,
258  iter_output);
259  }
260  }
261 
262  outputIterInfo(iter_output, it, effective_trial_stress, scalar, residual, reference_residual);
263 
264  ++it;
265  residual_old = residual;
266  scalar_old = scalar;
267  _residual_history[it % _num_resids] = residual;
268  }
269 
270  bool has_converged = true;
271  if (std::isnan(residual) || std::isinf(residual))
272  {
273  has_converged = false;
274  if (iter_output)
275  *iter_output << "Encountered inf or nan in material return mapping iterations." << std::endl;
276  }
277 
278  if (it == _fixed_max_its)
279  {
280  has_converged = false;
281  if (iter_output)
282  *iter_output << "Exceeded maximum iterations in material return mapping iterations."
283  << std::endl;
284  ;
285  }
286 
287  return has_converged;
288 }
289 
290 bool
292  Real & scalar,
293  std::stringstream * iter_output)
294 {
295  scalar = 0.0;
296  unsigned int it = 0;
297  Real residual = 10.0;
298  Real norm_residual = 10.0;
299  Real first_norm_residual = 10.0;
300  std::string iter_str;
301 
302  while (it < _max_its && norm_residual > _absolute_tolerance &&
303  (norm_residual / first_norm_residual) > _relative_tolerance)
304  {
305  residual = computeResidual(effective_trial_stress, scalar);
306  norm_residual = std::abs(residual);
307  if (it == 0)
308  {
309  first_norm_residual = norm_residual;
310  if (first_norm_residual == 0)
311  first_norm_residual = 1;
312  }
313 
314  scalar -= residual / computeDerivative(effective_trial_stress, scalar);
315 
316  if (iter_output)
317  {
318  iter_str = "Return mapping solve:\n iteration = " + Moose::stringify(it) + "\n" +
319  +" effective trial stress = " + Moose::stringify(effective_trial_stress) + "\n" +
320  +" scalar effective inelastic strain = " + Moose::stringify(scalar) + "\n" +
321  +" relative residual = " + Moose::stringify(norm_residual / first_norm_residual) +
322  "\n" + +" relative tolerance = " + Moose::stringify(_relative_tolerance) + "\n" +
323  +" absolute residual = " + Moose::stringify(norm_residual) + "\n" +
324  +" absolute tolerance = " + Moose::stringify(_absolute_tolerance) + "\n";
325  }
326  iterationFinalize(scalar);
327  ++it;
328  }
329 
330  if (iter_output)
331  *iter_output << iter_str;
332 
333  bool has_converged = true;
334 
335  if (it == _max_its && norm_residual > _absolute_tolerance &&
336  (norm_residual / first_norm_residual) > _relative_tolerance)
337  {
338  has_converged = false;
339  if (iter_output)
340  *iter_output << "Exceeded maximum iterations in material return mapping iterations."
341  << std::endl;
342  }
343  return has_converged;
344 }
345 
346 bool
347 SingleVariableReturnMappingSolution::converged(const Real residual, const Real reference)
348 {
349  return (std::abs(residual) <= _absolute_tolerance ||
350  (std::abs(residual) / reference) <= _relative_tolerance);
351 }
352 
353 bool
355  const Real residual,
356  const Real reference)
357 {
358  // Require that we have at least done _num_resids evaluations before we allow for
359  // acceptable convergence
360  if (it < _num_resids)
361  return false;
362 
363  // Check to see whether the residual has dropped by convergence_history_factor over
364  // the last _num_resids iterations. If it has (which means it's still making progress),
365  // don't consider it to be converged within the acceptable limits.
366  const Real convergence_history_factor = 10.0;
367  if (std::abs(residual * convergence_history_factor) <
368  std::abs(_residual_history[(it + 1) % _num_resids]))
369  return false;
370 
371  // Now that it's determined that progress is not being made, treat it as converged if
372  // we're within the acceptable convergence limits
373  return converged(residual / _acceptable_multiplier, reference);
374 }
375 
376 void
378  const unsigned int it,
379  const Real effective_trial_stress,
380  const Real scalar,
381  const Real residual,
382  const Real reference_residual)
383 {
384  if (iter_output)
385  {
386  *iter_output << " iteration=" << it << " trial_stress=" << effective_trial_stress
387  << " scalar=" << scalar << " residual=" << residual
388  << " ref_res=" << reference_residual
389  << " rel_res=" << std::abs(residual) / reference_residual
390  << " rel_tol=" << _relative_tolerance << " abs_res=" << std::abs(residual)
391  << " abs_tol=" << _absolute_tolerance << std::endl;
392  }
393 }
394 
395 void
397  Real & scalar_increment,
398  const Real scalar_old,
399  const Real min_permissible_scalar,
400  const Real max_permissible_scalar,
401  std::stringstream * iter_output)
402 {
403  if (scalar > max_permissible_scalar)
404  {
405  scalar_increment = (max_permissible_scalar - scalar_old) / 2.0;
406  scalar = scalar_old + scalar_increment;
407  if (iter_output)
408  *iter_output << "Scalar greater than maximum (" << max_permissible_scalar
409  << ") adjusted scalar=" << scalar << " scalar_increment=" << scalar_increment
410  << std::endl;
411  }
412  else if (scalar < min_permissible_scalar)
413  {
414  scalar_increment = (min_permissible_scalar - scalar_old) / 2.0;
415  scalar = scalar_old + scalar_increment;
416  if (iter_output)
417  *iter_output << "Scalar less than minimum (" << min_permissible_scalar
418  << ") adjusted scalar=" << scalar << " scalar_increment=" << scalar_increment
419  << std::endl;
420  }
421 }
422 
423 void
425  const Real residual,
426  const Real init_resid_sign,
427  Real & scalar_upper_bound,
428  Real & scalar_lower_bound,
429  std::stringstream * iter_output)
430 {
431  // Update upper/lower bounds as applicable
432  if (residual * init_resid_sign < 0.0 && scalar < scalar_upper_bound)
433  {
434  scalar_upper_bound = scalar;
435  if (scalar_upper_bound < scalar_lower_bound)
436  {
437  scalar_upper_bound = scalar_lower_bound;
438  scalar_lower_bound = 0.0;
439  if (iter_output)
440  *iter_output << " Corrected for scalar_upper_bound < scalar_lower_bound" << std::endl;
441  }
442  }
443  // Don't permit setting scalar_lower_bound > scalar_upper_bound (but do permit the reverse).
444  // This ensures that if we encounter multiple roots, we pick the lowest one.
445  else if (residual * init_resid_sign > 0.0 && scalar > scalar_lower_bound &&
446  scalar < scalar_upper_bound)
447  scalar_lower_bound = scalar;
448 }
SingleVariableReturnMappingSolution(const InputParameters &parameters)
bool internalSolveLegacy(const Real effective_trial_stress, Real &scalar, std::stringstream *iter_output)
Method called from within this class to perform the actual return mappping iterations.
bool internalSolve(const Real effective_trial_stress, Real &scalar, std::stringstream *iter_output)
Method called from within this class to perform the actual return mappping iterations.
InputParameters validParams< SingleVariableReturnMappingSolution >()
Real _acceptable_multiplier
Multiplier applied to relative and absolute tolerances for acceptable convergence.
bool _legacy_return_mapping
Whether to use the legacy return mapping algorithm and compute residuals in the legacy manner...
void outputIterInfo(std::stringstream *iter_output, const unsigned int it, const Real effective_trial_stress, const Real scalar, const Real residual, const Real reference_residual)
Output information about convergence history of the model.
const bool _output_iteration_info
Whether to output iteration information all the time (regardless of whether iterations converge) ...
bool converged(const Real residual, const Real reference)
Check to see whether the residual is within the convergence limits.
virtual Real computeDerivative(const Real effective_trial_stress, const Real scalar)=0
Compute the derivative of the residual as a function of the scalar variable.
void checkPermissibleRange(Real &scalar, Real &scalar_increment, const Real scalar_old, const Real min_permissible_scalar, const Real max_permissible_scalar, std::stringstream *iter_output)
Check to see whether solution is within admissible range, and set it within that range if it is not...
void updateBounds(const Real scalar, const Real residual, const Real init_resid_sign, Real &scalar_upper_bound, Real &scalar_lower_bound, std::stringstream *iter_output)
Update the upper and lower bounds of the root for the effective inelastic strain. ...
virtual void iterationFinalize(Real)
Finalize internal state variables for a model for a given iteration.
void returnMappingSolve(const Real effective_trial_stress, Real &scalar, const ConsoleStream &console)
Perform the return mapping iterations.
virtual Real maximumPermissibleValue(const Real effective_trial_stress) const
Compute the maximum permissible value of the scalar.
bool _line_search
Whether to use line searches to improve convergence.
virtual Real computeResidual(const Real effective_trial_stress, const Real scalar)=0
Compute the residual for a predicted value of the scalar.
std::vector< Real > _residual_history
History of residuals used to check whether progress is still being made on decreasing the residual...
Real _relative_tolerance
Relative convergence tolerance.
unsigned int _max_its
Maximum number of return mapping iterations (used only in legacy return mapping)
const std::size_t _num_resids
Number of residuals to be stored in history.
const unsigned int _fixed_max_its
Maximum number of return mapping iterations used in current procedure. Not settable by user...
bool _check_range
Whether to check to see whether iterative solution is within admissible range, and set within that ra...
Real _absolute_tolerance
Absolute convergence tolerance.
virtual Real computeReferenceResidual(const Real effective_trial_stress, const Real scalar)=0
Compute a reference quantity to be used for checking relative convergence.
bool convergedAcceptable(const unsigned int it, const Real residual, const Real reference)
Check to see whether the residual is within acceptable convergence limits.
bool _bracket_solution
Whether to save upper and lower bounds of root for scalar, and set solution to the midpoint between t...