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 3525 : SingleVariableReturnMappingSolutionTempl<is_ad>::validParams()
27 : {
28 3525 : InputParameters params = emptyInputParameters();
29 7050 : params.addParam<Real>(
30 7050 : "relative_tolerance", 1e-8, "Relative convergence tolerance for Newton iteration");
31 7050 : params.addParam<Real>(
32 7050 : "absolute_tolerance", 1e-11, "Absolute convergence tolerance for Newton iteration");
33 7050 : params.addParam<Real>("acceptable_multiplier",
34 7050 : 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 7050 : MooseEnum internal_solve_output_on_enum("never on_error always", "on_error");
41 7050 : params.addParam<MooseEnum>("internal_solve_output_on",
42 : internal_solve_output_on_enum,
43 : "When to output internal Newton solve information");
44 7050 : params.addParam<bool>("internal_solve_full_iteration_history",
45 7050 : 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 7050 : params.addParamNamesToGroup("internal_solve_output_on internal_solve_full_iteration_history",
50 : "Debug");
51 :
52 7050 : params.addParam<bool>("automatic_differentiation_return_mapping",
53 7050 : false,
54 : "Whether to use automatic differentiation to compute the derivative.");
55 3525 : return params;
56 3525 : }
57 :
58 : template <bool is_ad>
59 2645 : SingleVariableReturnMappingSolutionTempl<is_ad>::SingleVariableReturnMappingSolutionTempl(
60 : const InputParameters & parameters)
61 2645 : : _check_range(false),
62 2645 : _line_search(true),
63 2645 : _bracket_solution(true),
64 2645 : _internal_solve_output_on(
65 2645 : parameters.get<MooseEnum>("internal_solve_output_on").getEnum<InternalSolveOutput>()),
66 2645 : _max_its(1000), // Far larger than ever expected to be needed
67 2645 : _internal_solve_full_iteration_history(
68 2645 : parameters.get<bool>("internal_solve_full_iteration_history")),
69 2645 : _relative_tolerance(parameters.get<Real>("relative_tolerance")),
70 2645 : _absolute_tolerance(parameters.get<Real>("absolute_tolerance")),
71 2645 : _acceptable_multiplier(parameters.get<Real>("acceptable_multiplier")),
72 2645 : _ad_derivative(parameters.get<bool>("automatic_differentiation_return_mapping")),
73 2645 : _num_resids(30),
74 2645 : _residual_history(_num_resids, std::numeric_limits<Real>::max()),
75 2645 : _iteration(0),
76 2645 : _initial_residual(0.0),
77 2645 : _residual(0.0),
78 5290 : _svrms_name(parameters.getObjectName())
79 : {
80 2645 : }
81 :
82 : template <bool is_ad>
83 : GenericReal<is_ad>
84 131956 : SingleVariableReturnMappingSolutionTempl<is_ad>::minimumPermissibleValue(
85 : const GenericReal<is_ad> & /*effective_trial_stress*/) const
86 : {
87 131956 : return std::numeric_limits<Real>::lowest();
88 : }
89 :
90 : template <bool is_ad>
91 : GenericReal<is_ad>
92 131956 : SingleVariableReturnMappingSolutionTempl<is_ad>::maximumPermissibleValue(
93 : const GenericReal<is_ad> & /*effective_trial_stress*/) const
94 : {
95 131956 : return std::numeric_limits<Real>::max();
96 : }
97 :
98 : template <bool is_ad>
99 : void
100 32177129 : 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 12668 : std::unique_ptr<std::stringstream> iter_output =
107 32177129 : (_internal_solve_output_on == InternalSolveOutput::ALWAYS)
108 32177129 : ? 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 32177129 : internalSolve(effective_trial_stress,
115 : scalar,
116 32177129 : _internal_solve_full_iteration_history ? iter_output.get() : nullptr);
117 32177129 : if (solve_state != SolveState::SUCCESS &&
118 14 : _internal_solve_output_on != InternalSolveOutput::ALWAYS)
119 : {
120 : // output suppressed by user, throw immediately
121 14 : 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 14 : if (!iter_output)
126 28 : iter_output = std::make_unique<std::stringstream>();
127 :
128 : // add the appropriate error message to the output
129 14 : 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 14 : *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 14 : if (_internal_solve_full_iteration_history)
146 0 : internalSolve(effective_trial_stress, scalar, iter_output.get());
147 :
148 : // Append summary and throw exception
149 14 : outputIterationSummary(iter_output.get(), _iteration);
150 42 : mooseException(iter_output->str());
151 : }
152 :
153 32177115 : if (_internal_solve_output_on == InternalSolveOutput::ALWAYS)
154 : {
155 : // the solve did not fail but the user requested debug output anyways
156 12668 : outputIterationSummary(iter_output.get(), _iteration);
157 12668 : console << iter_output->str() << std::flush;
158 : }
159 32177129 : }
160 :
161 : template <bool is_ad>
162 : typename SingleVariableReturnMappingSolutionTempl<is_ad>::SolveState
163 32177129 : SingleVariableReturnMappingSolutionTempl<is_ad>::internalSolve(
164 : const GenericReal<is_ad> effective_trial_stress,
165 : GenericReal<is_ad> & scalar,
166 : std::stringstream * iter_output)
167 : {
168 32177129 : scalar = initialGuess(effective_trial_stress);
169 32177129 : GenericReal<is_ad> scalar_old = scalar;
170 32177129 : GenericReal<is_ad> scalar_increment = 0.0;
171 32177129 : const GenericReal<is_ad> min_permissible_scalar = minimumPermissibleValue(effective_trial_stress);
172 32177129 : const GenericReal<is_ad> max_permissible_scalar = maximumPermissibleValue(effective_trial_stress);
173 32177129 : GenericReal<is_ad> scalar_upper_bound = max_permissible_scalar;
174 32177129 : GenericReal<is_ad> scalar_lower_bound = min_permissible_scalar;
175 32177129 : _iteration = 0;
176 :
177 32177129 : computeResidualAndDerivativeHelper(effective_trial_stress, scalar);
178 32177129 : _initial_residual = _residual;
179 :
180 16018455 : GenericReal<is_ad> residual_old = _residual;
181 : Real init_resid_sign = MathUtils::sign(MetaPhysicL::raw_value(_residual));
182 32177129 : Real reference_residual = computeReferenceResidual(effective_trial_stress, scalar);
183 :
184 : if (converged(_residual, reference_residual))
185 : {
186 1552839 : iterationFinalize(scalar);
187 1552839 : outputIterationStep(iter_output, effective_trial_stress, scalar, reference_residual);
188 1155533 : return SolveState::SUCCESS;
189 : }
190 :
191 30624290 : _residual_history.assign(_num_resids, std::numeric_limits<Real>::max());
192 30624290 : _residual_history[0] = MetaPhysicL::raw_value(_residual);
193 :
194 139140489 : while (_iteration < _max_its && !converged(_residual, reference_residual) &&
195 138293583 : !convergedAcceptable(_iteration, reference_residual))
196 : {
197 138293446 : preStep(scalar_old, _residual, _derivative);
198 :
199 220749891 : scalar_increment = -_residual / _derivative;
200 138293446 : scalar = scalar_old + scalar_increment;
201 :
202 138293446 : if (_check_range)
203 951900 : checkPermissibleRange(scalar,
204 : scalar_increment,
205 : scalar_old,
206 : min_permissible_scalar,
207 : max_permissible_scalar,
208 : iter_output);
209 :
210 138293446 : computeResidualAndDerivativeHelper(effective_trial_stress, scalar);
211 138293446 : reference_residual = computeReferenceResidual(effective_trial_stress, scalar);
212 138293446 : iterationFinalize(scalar);
213 :
214 138293446 : if (_bracket_solution)
215 138293446 : updateBounds(
216 : scalar, _residual, init_resid_sign, scalar_upper_bound, scalar_lower_bound, iter_output);
217 :
218 : if (converged(_residual, reference_residual))
219 : {
220 29777247 : 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 108516199 : if (_line_search)
229 : {
230 108516199 : if (residual_old - _residual != 0.0)
231 : {
232 108504343 : GenericReal<is_ad> alpha = residual_old / (residual_old - _residual);
233 66827551 : alpha = MathUtils::clamp(alpha, 1.0e-2, 1.0);
234 :
235 67677248 : if (alpha != 1.0)
236 : {
237 : modified_increment = true;
238 852571 : scalar_increment *= alpha;
239 852571 : 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 108516199 : 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 242181227 : if (scalar_old + scalar_increment >= scalar_upper_bound ||
252 108513357 : scalar_old + scalar_increment <= scalar_lower_bound)
253 : {
254 107669304 : if (scalar_upper_bound != max_permissible_scalar &&
255 4702 : scalar_lower_bound != min_permissible_scalar)
256 : {
257 4782 : const Real frac = 0.5;
258 9484 : scalar_increment =
259 9484 : (1.0 - frac) * scalar_lower_bound + frac * scalar_upper_bound - scalar_old;
260 : modified_increment = true;
261 9484 : 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 108506715 : if (modified_increment)
272 : {
273 856397 : scalar = scalar_old + scalar_increment;
274 856397 : computeResidualAndDerivativeHelper(effective_trial_stress, scalar);
275 856397 : reference_residual = computeReferenceResidual(effective_trial_stress, scalar);
276 856397 : iterationFinalize(scalar);
277 :
278 856397 : if (_bracket_solution)
279 856397 : updateBounds(scalar,
280 : _residual,
281 : init_resid_sign,
282 : scalar_upper_bound,
283 : scalar_lower_bound,
284 : iter_output);
285 : }
286 : }
287 :
288 108516199 : outputIterationStep(iter_output, effective_trial_stress, scalar, reference_residual);
289 :
290 108516199 : ++_iteration;
291 108516199 : residual_old = _residual;
292 108516199 : scalar_old = scalar;
293 108516199 : _residual_history[_iteration % _num_resids] = MetaPhysicL::raw_value(_residual);
294 : }
295 :
296 : using std::isnan, std::isinf;
297 30624290 : if (isnan(_residual) || isinf(MetaPhysicL::raw_value(_residual)))
298 : return SolveState::NAN_INF;
299 :
300 30624290 : if (_iteration == _max_its)
301 : return SolveState::EXCEEDED_ITERATIONS;
302 :
303 : return SolveState::SUCCESS;
304 : }
305 :
306 : template <bool is_ad>
307 : void
308 171326972 : SingleVariableReturnMappingSolutionTempl<is_ad>::computeResidualAndDerivativeHelper(
309 : const GenericReal<is_ad> & effective_trial_stress, const GenericReal<is_ad> & scalar)
310 : {
311 171326972 : 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 171326972 : _residual = computeResidual(effective_trial_stress, scalar);
321 171326972 : _derivative = computeDerivative(effective_trial_stress, scalar);
322 : }
323 171326972 : }
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 309624135 : return (std::abs(residual) <= _absolute_tolerance ||
332 277961226 : std::abs(residual / reference) <= _relative_tolerance);
333 : }
334 :
335 : template <bool is_ad>
336 : bool
337 138293583 : 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 138293583 : 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 2795919 : const Real convergence_history_factor = 10.0;
351 7525472 : 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 21568 : return converged(_residual / _acceptable_multiplier, reference);
357 : }
358 :
359 : template <bool is_ad>
360 : void
361 951900 : 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 951900 : if (scalar > max_permissible_scalar)
370 : {
371 1763 : scalar_increment = (max_permissible_scalar - scalar_old) / 2.0;
372 883 : scalar = scalar_old + scalar_increment;
373 883 : 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 951017 : else if (scalar < min_permissible_scalar)
380 : {
381 10946 : scalar_increment = (min_permissible_scalar - scalar_old) / 2.0;
382 5554 : scalar = scalar_old + scalar_increment;
383 5554 : 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 951900 : }
389 :
390 : template <bool is_ad>
391 : void
392 139149843 : 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 139149843 : if (residual * init_resid_sign < 0.0 && scalar < scalar_upper_bound)
402 : {
403 6488279 : scalar_upper_bound = scalar;
404 6488279 : if (scalar_upper_bound < scalar_lower_bound)
405 : {
406 1026 : scalar_upper_bound = scalar_lower_bound;
407 1026 : scalar_lower_bound = 0.0;
408 1026 : 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 132661564 : else if (residual * init_resid_sign > 0.0 && scalar > scalar_lower_bound &&
415 50917866 : scalar < scalar_upper_bound)
416 131808442 : scalar_lower_bound = scalar;
417 139149843 : }
418 :
419 : template <bool is_ad>
420 : void
421 139846285 : 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 139846285 : 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 139846285 : }
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>;
|
