Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
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 : #include "SingleVariableReturnMappingSolution.h"
11 :
12 : #include "Moose.h"
13 : #include "MooseEnum.h"
14 : #include "MooseObject.h"
15 : #include "ConsoleStreamInterface.h"
16 : #include "Conversion.h"
17 : #include "MathUtils.h"
18 :
19 : #include <limits>
20 : #include <string>
21 : #include <cmath>
22 : #include <memory>
23 :
24 : template <bool is_ad>
25 : InputParameters
26 6598 : SingleVariableReturnMappingSolutionTempl<is_ad>::validParams()
27 : {
28 6598 : InputParameters params = emptyInputParameters();
29 13196 : params.addParam<Real>(
30 13196 : "relative_tolerance", 1e-8, "Relative convergence tolerance for Newton iteration");
31 13196 : params.addParam<Real>(
32 13196 : "absolute_tolerance", 1e-11, "Absolute convergence tolerance for Newton iteration");
33 13196 : params.addParam<Real>("acceptable_multiplier",
34 13196 : 10,
35 : "Factor applied to relative and absolute "
36 : "tolerance for acceptable convergence if "
37 : "iterations are no longer making progress");
38 :
39 : // diagnostic output parameters
40 13196 : MooseEnum internal_solve_output_on_enum("never on_error always", "on_error");
41 13196 : params.addParam<MooseEnum>("internal_solve_output_on",
42 : internal_solve_output_on_enum,
43 : "When to output internal Newton solve information");
44 13196 : params.addParam<bool>("internal_solve_full_iteration_history",
45 13196 : false,
46 : "Set true to output full internal Newton iteration history at times "
47 : "determined by `internal_solve_output_on`. If false, only a summary is "
48 : "output.");
49 13196 : params.addParamNamesToGroup("internal_solve_output_on internal_solve_full_iteration_history",
50 : "Debug");
51 :
52 13196 : params.addParam<bool>("automatic_differentiation_return_mapping",
53 13196 : false,
54 : "Whether to use automatic differentiation to compute the derivative.");
55 6598 : return params;
56 6598 : }
57 :
58 : template <bool is_ad>
59 4951 : SingleVariableReturnMappingSolutionTempl<is_ad>::SingleVariableReturnMappingSolutionTempl(
60 : const InputParameters & parameters)
61 4951 : : _check_range(false),
62 4951 : _line_search(true),
63 4951 : _bracket_solution(true),
64 4951 : _internal_solve_output_on(
65 4951 : parameters.get<MooseEnum>("internal_solve_output_on").getEnum<InternalSolveOutput>()),
66 4951 : _max_its(1000), // Far larger than ever expected to be needed
67 4951 : _internal_solve_full_iteration_history(
68 4951 : parameters.get<bool>("internal_solve_full_iteration_history")),
69 4951 : _relative_tolerance(parameters.get<Real>("relative_tolerance")),
70 4951 : _absolute_tolerance(parameters.get<Real>("absolute_tolerance")),
71 4951 : _acceptable_multiplier(parameters.get<Real>("acceptable_multiplier")),
72 4951 : _ad_derivative(parameters.get<bool>("automatic_differentiation_return_mapping")),
73 4951 : _num_resids(30),
74 4951 : _residual_history(_num_resids, std::numeric_limits<Real>::max()),
75 4951 : _iteration(0),
76 4951 : _initial_residual(0.0),
77 4951 : _residual(0.0),
78 9902 : _svrms_name(parameters.getObjectName())
79 : {
80 4951 : }
81 :
82 : template <bool is_ad>
83 : GenericReal<is_ad>
84 226787 : SingleVariableReturnMappingSolutionTempl<is_ad>::minimumPermissibleValue(
85 : const GenericReal<is_ad> & /*effective_trial_stress*/) const
86 : {
87 226787 : return std::numeric_limits<Real>::lowest();
88 : }
89 :
90 : template <bool is_ad>
91 : GenericReal<is_ad>
92 226787 : SingleVariableReturnMappingSolutionTempl<is_ad>::maximumPermissibleValue(
93 : const GenericReal<is_ad> & /*effective_trial_stress*/) const
94 : {
95 226787 : return std::numeric_limits<Real>::max();
96 : }
97 :
98 : template <bool is_ad>
99 : void
100 57354569 : SingleVariableReturnMappingSolutionTempl<is_ad>::returnMappingSolve(
101 : const GenericReal<is_ad> & effective_trial_stress,
102 : GenericReal<is_ad> & scalar,
103 : const ConsoleStream & console)
104 : {
105 : // construct the stringstream here only if the debug level is set to ALL
106 25336 : std::unique_ptr<std::stringstream> iter_output =
107 57354569 : (_internal_solve_output_on == InternalSolveOutput::ALWAYS)
108 57354569 : ? std::make_unique<std::stringstream>()
109 : : nullptr;
110 :
111 : // do the internal solve and capture iteration info during the first round
112 : // iff full history output is requested regardless of whether the solve failed or succeeded
113 : auto solve_state =
114 57354569 : internalSolve(effective_trial_stress,
115 : scalar,
116 57354569 : _internal_solve_full_iteration_history ? iter_output.get() : nullptr);
117 57354569 : if (solve_state != SolveState::SUCCESS &&
118 28 : _internal_solve_output_on != InternalSolveOutput::ALWAYS)
119 : {
120 : // output suppressed by user, throw immediately
121 28 : if (_internal_solve_output_on == InternalSolveOutput::NEVER)
122 0 : mooseException("");
123 :
124 : // user expects some kind of output, if necessary setup output stream now
125 28 : if (!iter_output)
126 56 : iter_output = std::make_unique<std::stringstream>();
127 :
128 : // add the appropriate error message to the output
129 28 : switch (solve_state)
130 : {
131 : case SolveState::NAN_INF:
132 0 : *iter_output << "Encountered inf or nan in material return mapping iterations.\n";
133 : break;
134 :
135 : case SolveState::EXCEEDED_ITERATIONS:
136 28 : *iter_output << "Exceeded maximum iterations in material return mapping iterations.\n";
137 : break;
138 :
139 0 : default:
140 0 : mooseError("Unhandled solver state");
141 : }
142 :
143 : // if full history output is only requested for failed solves we have to repeat
144 : // the solve a second time
145 28 : if (_internal_solve_full_iteration_history)
146 0 : internalSolve(effective_trial_stress, scalar, iter_output.get());
147 :
148 : // Append summary and throw exception
149 28 : outputIterationSummary(iter_output.get(), _iteration);
150 84 : mooseException(iter_output->str());
151 : }
152 :
153 57354541 : if (_internal_solve_output_on == InternalSolveOutput::ALWAYS)
154 : {
155 : // the solve did not fail but the user requested debug output anyways
156 25336 : outputIterationSummary(iter_output.get(), _iteration);
157 25336 : console << iter_output->str() << std::flush;
158 : }
159 57354569 : }
160 :
161 : template <bool is_ad>
162 : typename SingleVariableReturnMappingSolutionTempl<is_ad>::SolveState
163 57354569 : SingleVariableReturnMappingSolutionTempl<is_ad>::internalSolve(
164 : const GenericReal<is_ad> effective_trial_stress,
165 : GenericReal<is_ad> & scalar,
166 : std::stringstream * iter_output)
167 : {
168 57354569 : scalar = initialGuess(effective_trial_stress);
169 57354569 : GenericReal<is_ad> scalar_old = scalar;
170 57354569 : GenericReal<is_ad> scalar_increment = 0.0;
171 57354569 : const GenericReal<is_ad> min_permissible_scalar = minimumPermissibleValue(effective_trial_stress);
172 57354569 : const GenericReal<is_ad> max_permissible_scalar = maximumPermissibleValue(effective_trial_stress);
173 57354569 : GenericReal<is_ad> scalar_upper_bound = max_permissible_scalar;
174 57354569 : GenericReal<is_ad> scalar_lower_bound = min_permissible_scalar;
175 57354569 : _iteration = 0;
176 :
177 57354569 : computeResidualAndDerivativeHelper(effective_trial_stress, scalar);
178 57354569 : _initial_residual = _residual;
179 :
180 29231216 : GenericReal<is_ad> residual_old = _residual;
181 : Real init_resid_sign = MathUtils::sign(MetaPhysicL::raw_value(_residual));
182 57354569 : Real reference_residual = computeReferenceResidual(effective_trial_stress, scalar);
183 :
184 : if (converged(_residual, reference_residual))
185 : {
186 2724972 : iterationFinalize(scalar);
187 2724972 : outputIterationStep(iter_output, effective_trial_stress, scalar, reference_residual);
188 1976865 : return SolveState::SUCCESS;
189 : }
190 :
191 54629597 : _residual_history.assign(_num_resids, std::numeric_limits<Real>::max());
192 54629597 : _residual_history[0] = MetaPhysicL::raw_value(_residual);
193 :
194 389614191 : while (_iteration < _max_its && !converged(_residual, reference_residual) &&
195 241287676 : !convergedAcceptable(_iteration, reference_residual))
196 : {
197 241287418 : preStep(scalar_old, _residual, _derivative);
198 :
199 388199721 : scalar_increment = -_residual / _derivative;
200 241287418 : scalar = scalar_old + scalar_increment;
201 :
202 241287418 : if (_check_range)
203 1793261 : checkPermissibleRange(scalar,
204 : scalar_increment,
205 : scalar_old,
206 : min_permissible_scalar,
207 : max_permissible_scalar,
208 : iter_output);
209 :
210 241287418 : computeResidualAndDerivativeHelper(effective_trial_stress, scalar);
211 241287418 : reference_residual = computeReferenceResidual(effective_trial_stress, scalar);
212 241287418 : iterationFinalize(scalar);
213 :
214 241287418 : if (_bracket_solution)
215 241287418 : updateBounds(
216 : scalar, _residual, init_resid_sign, scalar_upper_bound, scalar_lower_bound, iter_output);
217 :
218 : if (converged(_residual, reference_residual))
219 : {
220 53215257 : outputIterationStep(iter_output, effective_trial_stress, scalar, reference_residual);
221 : break;
222 : }
223 : else
224 : {
225 : bool modified_increment = false;
226 :
227 : // Line Search
228 188072161 : if (_line_search)
229 : {
230 188072161 : if (residual_old - _residual != 0.0)
231 : {
232 188049057 : GenericReal<is_ad> alpha = residual_old / (residual_old - _residual);
233 118413872 : alpha = MathUtils::clamp(alpha, 1.0e-2, 1.0);
234 :
235 119833536 : if (alpha != 1.0)
236 : {
237 : modified_increment = true;
238 1425349 : scalar_increment *= alpha;
239 1425349 : if (iter_output)
240 0 : *iter_output << " Line search alpha = " << MetaPhysicL::raw_value(alpha)
241 0 : << " increment = " << MetaPhysicL::raw_value(scalar_increment)
242 : << std::endl;
243 : }
244 : }
245 : }
246 :
247 188072161 : if (_bracket_solution)
248 : {
249 : // Check to see whether trial scalar_increment is outside the bounds, and set it to a point
250 : // within the bounds if it is
251 424915953 : if (scalar_old + scalar_increment >= scalar_upper_bound ||
252 188068417 : scalar_old + scalar_increment <= scalar_lower_bound)
253 : {
254 186658101 : if (scalar_upper_bound != max_permissible_scalar &&
255 9348 : scalar_lower_bound != min_permissible_scalar)
256 : {
257 9406 : const Real frac = 0.5;
258 18754 : scalar_increment =
259 18754 : (1.0 - frac) * scalar_lower_bound + frac * scalar_upper_bound - scalar_old;
260 : modified_increment = true;
261 18754 : if (iter_output)
262 0 : *iter_output << " Trial scalar_increment exceeded bounds. Setting between "
263 : "lower/upper bounds. frac: "
264 : << frac << std::endl;
265 : }
266 : }
267 : }
268 :
269 : // Update the trial scalar and recompute residual if the line search or bounds checking
270 : // modified the increment
271 188053407 : if (modified_increment)
272 : {
273 1432829 : scalar = scalar_old + scalar_increment;
274 1432829 : computeResidualAndDerivativeHelper(effective_trial_stress, scalar);
275 1432829 : reference_residual = computeReferenceResidual(effective_trial_stress, scalar);
276 1432829 : iterationFinalize(scalar);
277 :
278 1432829 : if (_bracket_solution)
279 1432829 : updateBounds(scalar,
280 : _residual,
281 : init_resid_sign,
282 : scalar_upper_bound,
283 : scalar_lower_bound,
284 : iter_output);
285 : }
286 : }
287 :
288 188072161 : outputIterationStep(iter_output, effective_trial_stress, scalar, reference_residual);
289 :
290 188072161 : ++_iteration;
291 188072161 : residual_old = _residual;
292 188072161 : scalar_old = scalar;
293 188072161 : _residual_history[_iteration % _num_resids] = MetaPhysicL::raw_value(_residual);
294 : }
295 :
296 : using std::isnan, std::isinf;
297 54629597 : if (isnan(_residual) || isinf(MetaPhysicL::raw_value(_residual)))
298 : return SolveState::NAN_INF;
299 :
300 54629597 : if (_iteration == _max_its)
301 : return SolveState::EXCEEDED_ITERATIONS;
302 :
303 : return SolveState::SUCCESS;
304 : }
305 :
306 : template <bool is_ad>
307 : void
308 300074816 : SingleVariableReturnMappingSolutionTempl<is_ad>::computeResidualAndDerivativeHelper(
309 : const GenericReal<is_ad> & effective_trial_stress, const GenericReal<is_ad> & scalar)
310 : {
311 300074816 : if (_ad_derivative)
312 : {
313 0 : GenericChainedReal<is_ad> residual_and_derivative =
314 0 : computeResidualAndDerivative(effective_trial_stress, GenericChainedReal<is_ad>(scalar, 1));
315 0 : _residual = residual_and_derivative.value();
316 0 : _derivative = residual_and_derivative.derivatives();
317 : }
318 : else
319 : {
320 300074816 : _residual = computeResidual(effective_trial_stress, scalar);
321 300074816 : _derivative = computeDerivative(effective_trial_stress, scalar);
322 : }
323 300074816 : }
324 :
325 : template <bool is_ad>
326 : bool
327 0 : SingleVariableReturnMappingSolutionTempl<is_ad>::converged(const GenericReal<is_ad> & ad_residual,
328 : const Real reference)
329 : {
330 : const Real residual = MetaPhysicL::raw_value(ad_residual);
331 541369875 : return (std::abs(residual) <= _absolute_tolerance ||
332 484874118 : std::abs(residual / reference) <= _relative_tolerance);
333 : }
334 :
335 : template <bool is_ad>
336 : bool
337 241287676 : SingleVariableReturnMappingSolutionTempl<is_ad>::convergedAcceptable(const unsigned int it,
338 : const Real reference)
339 : {
340 : using std::abs;
341 :
342 : // Require that we have at least done _num_resids evaluations before we allow for
343 : // acceptable convergence
344 241287676 : if (it < _num_resids)
345 : return false;
346 :
347 : // Check to see whether the residual has dropped by convergence_history_factor over
348 : // the last _num_resids iterations. If it has (which means it's still making progress),
349 : // don't consider it to be converged within the acceptable limits.
350 4323496 : const Real convergence_history_factor = 10.0;
351 11612736 : if (abs(_residual * convergence_history_factor) < abs(_residual_history[(it + 1) % _num_resids]))
352 : return false;
353 :
354 : // Now that it's determined that progress is not being made, treat it as converged if
355 : // we're within the acceptable convergence limits
356 43118 : return converged(_residual / _acceptable_multiplier, reference);
357 : }
358 :
359 : template <bool is_ad>
360 : void
361 1793261 : SingleVariableReturnMappingSolutionTempl<is_ad>::checkPermissibleRange(
362 : GenericReal<is_ad> & scalar,
363 : GenericReal<is_ad> & scalar_increment,
364 : const GenericReal<is_ad> & scalar_old,
365 : const GenericReal<is_ad> min_permissible_scalar,
366 : const GenericReal<is_ad> max_permissible_scalar,
367 : std::stringstream * iter_output)
368 : {
369 1793261 : if (scalar > max_permissible_scalar)
370 : {
371 3520 : scalar_increment = (max_permissible_scalar - scalar_old) / 2.0;
372 1760 : scalar = scalar_old + scalar_increment;
373 1760 : if (iter_output)
374 1760 : *iter_output << "Scalar greater than maximum ("
375 : << MetaPhysicL::raw_value(max_permissible_scalar)
376 1760 : << ") adjusted scalar=" << MetaPhysicL::raw_value(scalar)
377 1760 : << " scalar_increment=" << MetaPhysicL::raw_value(scalar_increment) << std::endl;
378 : }
379 1791501 : else if (scalar < min_permissible_scalar)
380 : {
381 21728 : scalar_increment = (min_permissible_scalar - scalar_old) / 2.0;
382 10944 : scalar = scalar_old + scalar_increment;
383 10944 : if (iter_output)
384 2784 : *iter_output << "Scalar less than minimum (" << MetaPhysicL::raw_value(min_permissible_scalar)
385 2784 : << ") adjusted scalar=" << MetaPhysicL::raw_value(scalar)
386 2784 : << " scalar_increment=" << MetaPhysicL::raw_value(scalar_increment) << std::endl;
387 : }
388 1793261 : }
389 :
390 : template <bool is_ad>
391 : void
392 242720247 : SingleVariableReturnMappingSolutionTempl<is_ad>::updateBounds(
393 : const GenericReal<is_ad> & scalar,
394 : const GenericReal<is_ad> & residual,
395 : const Real init_resid_sign,
396 : GenericReal<is_ad> & scalar_upper_bound,
397 : GenericReal<is_ad> & scalar_lower_bound,
398 : std::stringstream * iter_output)
399 : {
400 : // Update upper/lower bounds as applicable
401 242720247 : if (residual * init_resid_sign < 0.0 && scalar < scalar_upper_bound)
402 : {
403 11873469 : scalar_upper_bound = scalar;
404 11873469 : if (scalar_upper_bound < scalar_lower_bound)
405 : {
406 1760 : scalar_upper_bound = scalar_lower_bound;
407 1760 : scalar_lower_bound = 0.0;
408 1760 : if (iter_output)
409 0 : *iter_output << " Corrected for scalar_upper_bound < scalar_lower_bound" << std::endl;
410 : }
411 : }
412 : // Don't permit setting scalar_lower_bound > scalar_upper_bound (but do permit the reverse).
413 : // This ensures that if we encounter multiple roots, we pick the lowest one.
414 230846778 : else if (residual * init_resid_sign > 0.0 && scalar > scalar_lower_bound &&
415 85538622 : scalar < scalar_upper_bound)
416 229326201 : scalar_lower_bound = scalar;
417 242720247 : }
418 :
419 : template <bool is_ad>
420 : void
421 244012390 : SingleVariableReturnMappingSolutionTempl<is_ad>::outputIterationStep(
422 : std::stringstream * iter_output,
423 : const GenericReal<is_ad> & effective_trial_stress,
424 : const GenericReal<is_ad> & scalar,
425 : const Real reference_residual)
426 : {
427 244012390 : if (iter_output)
428 : {
429 997128 : const unsigned int it = _iteration;
430 : const Real residual = MetaPhysicL::raw_value(_residual);
431 :
432 997128 : *iter_output << " iteration=" << it
433 997128 : << " trial_stress=" << MetaPhysicL::raw_value(effective_trial_stress)
434 1994256 : << " scalar=" << MetaPhysicL::raw_value(scalar) << " residual=" << residual
435 997128 : << " ref_res=" << reference_residual
436 997128 : << " rel_res=" << std::abs(residual) / reference_residual
437 997128 : << " rel_tol=" << _relative_tolerance << " abs_res=" << std::abs(residual)
438 997128 : << " abs_tol=" << _absolute_tolerance << '\n';
439 : }
440 244012390 : }
441 :
442 : template <bool is_ad>
443 : void
444 25364 : SingleVariableReturnMappingSolutionTempl<is_ad>::outputIterationSummary(
445 : std::stringstream * iter_output, const unsigned int total_it)
446 : {
447 25364 : if (iter_output)
448 50728 : *iter_output << "In " << total_it << " iterations the residual went from "
449 25364 : << MetaPhysicL::raw_value(_initial_residual) << " to "
450 50728 : << MetaPhysicL::raw_value(_residual) << " in '" << _svrms_name << "'."
451 : << std::endl;
452 25364 : }
453 :
454 : template class SingleVariableReturnMappingSolutionTempl<false>;
455 : template class SingleVariableReturnMappingSolutionTempl<true>;
|
