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