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 5834 : SingleVariableReturnMappingSolutionTempl<is_ad>::validParams()
27 : {
28 5834 : InputParameters params = emptyInputParameters();
29 11668 : params.addParam<Real>(
30 11668 : "relative_tolerance", 1e-8, "Relative convergence tolerance for Newton iteration");
31 11668 : params.addParam<Real>(
32 11668 : "absolute_tolerance", 1e-11, "Absolute convergence tolerance for Newton iteration");
33 11668 : params.addParam<Real>("acceptable_multiplier",
34 11668 : 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 11668 : MooseEnum internal_solve_output_on_enum("never on_error always", "on_error");
41 11668 : params.addParam<MooseEnum>("internal_solve_output_on",
42 : internal_solve_output_on_enum,
43 : "When to output internal Newton solve information");
44 11668 : params.addParam<bool>("internal_solve_full_iteration_history",
45 11668 : 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 11668 : params.addParamNamesToGroup("internal_solve_output_on internal_solve_full_iteration_history",
50 : "Debug");
51 :
52 11668 : params.addParam<bool>("automatic_differentiation_return_mapping",
53 11668 : false,
54 : "Whether to use automatic differentiation to compute the derivative.");
55 5834 : return params;
56 5834 : }
57 :
58 : template <bool is_ad>
59 4378 : SingleVariableReturnMappingSolutionTempl<is_ad>::SingleVariableReturnMappingSolutionTempl(
60 : const InputParameters & parameters)
61 4378 : : _check_range(false),
62 4378 : _line_search(true),
63 4378 : _bracket_solution(true),
64 4378 : _internal_solve_output_on(
65 4378 : parameters.get<MooseEnum>("internal_solve_output_on").getEnum<InternalSolveOutput>()),
66 4378 : _max_its(1000), // Far larger than ever expected to be needed
67 4378 : _internal_solve_full_iteration_history(
68 4378 : parameters.get<bool>("internal_solve_full_iteration_history")),
69 4378 : _relative_tolerance(parameters.get<Real>("relative_tolerance")),
70 4378 : _absolute_tolerance(parameters.get<Real>("absolute_tolerance")),
71 4378 : _acceptable_multiplier(parameters.get<Real>("acceptable_multiplier")),
72 4378 : _ad_derivative(parameters.get<bool>("automatic_differentiation_return_mapping")),
73 4378 : _num_resids(30),
74 4378 : _residual_history(_num_resids, std::numeric_limits<Real>::max()),
75 4378 : _iteration(0),
76 4378 : _initial_residual(0.0),
77 4378 : _residual(0.0),
78 8756 : _svrms_name(parameters.get<std::string>("_object_name"))
79 : {
80 4378 : }
81 :
82 : template <bool is_ad>
83 : GenericReal<is_ad>
84 189266 : SingleVariableReturnMappingSolutionTempl<is_ad>::minimumPermissibleValue(
85 : const GenericReal<is_ad> & /*effective_trial_stress*/) const
86 : {
87 189266 : return std::numeric_limits<Real>::lowest();
88 : }
89 :
90 : template <bool is_ad>
91 : GenericReal<is_ad>
92 189266 : SingleVariableReturnMappingSolutionTempl<is_ad>::maximumPermissibleValue(
93 : const GenericReal<is_ad> & /*effective_trial_stress*/) const
94 : {
95 189266 : return std::numeric_limits<Real>::max();
96 : }
97 :
98 : template <bool is_ad>
99 : void
100 56988974 : 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 56988974 : (_internal_solve_output_on == InternalSolveOutput::ALWAYS)
108 56988974 : ? 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 56988974 : internalSolve(effective_trial_stress,
115 : scalar,
116 56988974 : _internal_solve_full_iteration_history ? iter_output.get() : nullptr);
117 56988974 : 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 56988946 : 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 56988974 : }
160 :
161 : template <bool is_ad>
162 : typename SingleVariableReturnMappingSolutionTempl<is_ad>::SolveState
163 56988974 : SingleVariableReturnMappingSolutionTempl<is_ad>::internalSolve(
164 : const GenericReal<is_ad> effective_trial_stress,
165 : GenericReal<is_ad> & scalar,
166 : std::stringstream * iter_output)
167 : {
168 56988974 : scalar = initialGuess(effective_trial_stress);
169 56988974 : GenericReal<is_ad> scalar_old = scalar;
170 56988974 : GenericReal<is_ad> scalar_increment = 0.0;
171 56988974 : const GenericReal<is_ad> min_permissible_scalar = minimumPermissibleValue(effective_trial_stress);
172 56988974 : const GenericReal<is_ad> max_permissible_scalar = maximumPermissibleValue(effective_trial_stress);
173 56988974 : GenericReal<is_ad> scalar_upper_bound = max_permissible_scalar;
174 56988974 : GenericReal<is_ad> scalar_lower_bound = min_permissible_scalar;
175 56988974 : _iteration = 0;
176 :
177 56988974 : computeResidualAndDerivativeHelper(effective_trial_stress, scalar);
178 56988974 : _initial_residual = _residual;
179 :
180 26949530 : GenericReal<is_ad> residual_old = _residual;
181 : Real init_resid_sign = MathUtils::sign(MetaPhysicL::raw_value(_residual));
182 56988974 : Real reference_residual = computeReferenceResidual(effective_trial_stress, scalar);
183 :
184 : if (converged(_residual, reference_residual))
185 : {
186 2937770 : iterationFinalize(scalar);
187 2937770 : outputIterationStep(iter_output, effective_trial_stress, scalar, reference_residual);
188 2201880 : return SolveState::SUCCESS;
189 : }
190 :
191 54051204 : _residual_history.assign(_num_resids, std::numeric_limits<Real>::max());
192 54051204 : _residual_history[0] = MetaPhysicL::raw_value(_residual);
193 :
194 375393424 : while (_iteration < _max_its && !converged(_residual, reference_residual) &&
195 237199118 : !convergedAcceptable(_iteration, reference_residual))
196 : {
197 237198860 : preStep(scalar_old, _residual, _derivative);
198 :
199 374198344 : scalar_increment = -_residual / _derivative;
200 237198860 : scalar = scalar_old + scalar_increment;
201 :
202 237198860 : if (_check_range)
203 3152056 : checkPermissibleRange(scalar,
204 : scalar_increment,
205 : scalar_old,
206 : min_permissible_scalar,
207 : max_permissible_scalar,
208 : iter_output);
209 :
210 237198860 : computeResidualAndDerivativeHelper(effective_trial_stress, scalar);
211 237198860 : reference_residual = computeReferenceResidual(effective_trial_stress, scalar);
212 237198860 : iterationFinalize(scalar);
213 :
214 237198860 : if (_bracket_solution)
215 237198860 : updateBounds(
216 : scalar, _residual, init_resid_sign, scalar_upper_bound, scalar_lower_bound, iter_output);
217 :
218 : if (converged(_residual, reference_residual))
219 : {
220 52856254 : 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 184342606 : if (_line_search)
229 : {
230 184342606 : if (residual_old - _residual != 0.0)
231 : {
232 184319534 : GenericReal<is_ad> alpha = residual_old / (residual_old - _residual);
233 110770522 : alpha = MathUtils::clamp(alpha, 1.0e-2, 1.0);
234 :
235 184319534 : if (alpha != 1.0)
236 : {
237 : modified_increment = true;
238 1205956 : scalar_increment *= alpha;
239 1205956 : 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 184342606 : 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 405899698 : if (scalar_old + scalar_increment >= scalar_upper_bound ||
252 184338862 : scalar_old + scalar_increment <= scalar_lower_bound)
253 : {
254 183147936 : 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 184323852 : if (modified_increment)
272 : {
273 1213436 : scalar = scalar_old + scalar_increment;
274 1213436 : computeResidualAndDerivativeHelper(effective_trial_stress, scalar);
275 1213436 : reference_residual = computeReferenceResidual(effective_trial_stress, scalar);
276 1213436 : iterationFinalize(scalar);
277 :
278 1213436 : if (_bracket_solution)
279 1213436 : updateBounds(scalar,
280 : _residual,
281 : init_resid_sign,
282 : scalar_upper_bound,
283 : scalar_lower_bound,
284 : iter_output);
285 : }
286 : }
287 :
288 184342606 : outputIterationStep(iter_output, effective_trial_stress, scalar, reference_residual);
289 :
290 184342606 : ++_iteration;
291 184342606 : residual_old = _residual;
292 184342606 : scalar_old = scalar;
293 184342606 : _residual_history[_iteration % _num_resids] = MetaPhysicL::raw_value(_residual);
294 : }
295 :
296 54051204 : if (std::isnan(_residual) || std::isinf(MetaPhysicL::raw_value(_residual)))
297 : return SolveState::NAN_INF;
298 :
299 54051204 : if (_iteration == _max_its)
300 : return SolveState::EXCEEDED_ITERATIONS;
301 :
302 : return SolveState::SUCCESS;
303 : }
304 :
305 : template <bool is_ad>
306 : void
307 295401270 : SingleVariableReturnMappingSolutionTempl<is_ad>::computeResidualAndDerivativeHelper(
308 : const GenericReal<is_ad> & effective_trial_stress, const GenericReal<is_ad> & scalar)
309 : {
310 295401270 : if (_ad_derivative)
311 : {
312 6460456 : GenericChainedReal<is_ad> residual_and_derivative =
313 6460456 : computeResidualAndDerivative(effective_trial_stress, GenericChainedReal<is_ad>(scalar, 1));
314 6460456 : _residual = residual_and_derivative.value();
315 6460456 : _derivative = residual_and_derivative.derivatives();
316 : }
317 : else
318 : {
319 288940814 : _residual = computeResidual(effective_trial_stress, scalar);
320 288940814 : _derivative = computeDerivative(effective_trial_stress, scalar);
321 : }
322 295401270 : }
323 :
324 : template <bool is_ad>
325 : bool
326 0 : SingleVariableReturnMappingSolutionTempl<is_ad>::converged(const GenericReal<is_ad> & ad_residual,
327 : const Real reference)
328 : {
329 : const Real residual = MetaPhysicL::raw_value(ad_residual);
330 532607774 : return (std::abs(residual) <= _absolute_tolerance ||
331 476574360 : std::abs(residual / reference) <= _relative_tolerance);
332 : }
333 :
334 : template <bool is_ad>
335 : bool
336 237199118 : SingleVariableReturnMappingSolutionTempl<is_ad>::convergedAcceptable(const unsigned int it,
337 : const Real reference)
338 : {
339 : // Require that we have at least done _num_resids evaluations before we allow for
340 : // acceptable convergence
341 237199118 : if (it < _num_resids)
342 : return false;
343 :
344 : // Check to see whether the residual has dropped by convergence_history_factor over
345 : // the last _num_resids iterations. If it has (which means it's still making progress),
346 : // don't consider it to be converged within the acceptable limits.
347 4721312 : const Real convergence_history_factor = 10.0;
348 11507404 : if (std::abs(_residual * convergence_history_factor) <
349 6786092 : std::abs(_residual_history[(it + 1) % _num_resids]))
350 : return false;
351 :
352 : // Now that it's determined that progress is not being made, treat it as converged if
353 : // we're within the acceptable convergence limits
354 43118 : return converged(_residual / _acceptable_multiplier, reference);
355 : }
356 :
357 : template <bool is_ad>
358 : void
359 3152056 : SingleVariableReturnMappingSolutionTempl<is_ad>::checkPermissibleRange(
360 : GenericReal<is_ad> & scalar,
361 : GenericReal<is_ad> & scalar_increment,
362 : const GenericReal<is_ad> & scalar_old,
363 : const GenericReal<is_ad> min_permissible_scalar,
364 : const GenericReal<is_ad> max_permissible_scalar,
365 : std::stringstream * iter_output)
366 : {
367 3152056 : if (scalar > max_permissible_scalar)
368 : {
369 3520 : scalar_increment = (max_permissible_scalar - scalar_old) / 2.0;
370 1760 : scalar = scalar_old + scalar_increment;
371 1760 : if (iter_output)
372 1760 : *iter_output << "Scalar greater than maximum ("
373 : << MetaPhysicL::raw_value(max_permissible_scalar)
374 1760 : << ") adjusted scalar=" << MetaPhysicL::raw_value(scalar)
375 1760 : << " scalar_increment=" << MetaPhysicL::raw_value(scalar_increment) << std::endl;
376 : }
377 3150296 : else if (scalar < min_permissible_scalar)
378 : {
379 21696 : scalar_increment = (min_permissible_scalar - scalar_old) / 2.0;
380 10912 : scalar = scalar_old + scalar_increment;
381 10912 : if (iter_output)
382 2784 : *iter_output << "Scalar less than minimum (" << MetaPhysicL::raw_value(min_permissible_scalar)
383 2784 : << ") adjusted scalar=" << MetaPhysicL::raw_value(scalar)
384 2784 : << " scalar_increment=" << MetaPhysicL::raw_value(scalar_increment) << std::endl;
385 : }
386 3152056 : }
387 :
388 : template <bool is_ad>
389 : void
390 238412296 : SingleVariableReturnMappingSolutionTempl<is_ad>::updateBounds(
391 : const GenericReal<is_ad> & scalar,
392 : const GenericReal<is_ad> & residual,
393 : const Real init_resid_sign,
394 : GenericReal<is_ad> & scalar_upper_bound,
395 : GenericReal<is_ad> & scalar_lower_bound,
396 : std::stringstream * iter_output)
397 : {
398 : // Update upper/lower bounds as applicable
399 238412296 : if (residual * init_resid_sign < 0.0 && scalar < scalar_upper_bound)
400 : {
401 12183408 : scalar_upper_bound = scalar;
402 12183408 : if (scalar_upper_bound < scalar_lower_bound)
403 : {
404 1408 : scalar_upper_bound = scalar_lower_bound;
405 1408 : scalar_lower_bound = 0.0;
406 1408 : if (iter_output)
407 0 : *iter_output << " Corrected for scalar_upper_bound < scalar_lower_bound" << std::endl;
408 : }
409 : }
410 : // Don't permit setting scalar_lower_bound > scalar_upper_bound (but do permit the reverse).
411 : // This ensures that if we encounter multiple roots, we pick the lowest one.
412 226228888 : else if (residual * init_resid_sign > 0.0 && scalar > scalar_lower_bound &&
413 90998980 : scalar < scalar_upper_bound)
414 224708468 : scalar_lower_bound = scalar;
415 238412296 : }
416 :
417 : template <bool is_ad>
418 : void
419 240136630 : SingleVariableReturnMappingSolutionTempl<is_ad>::outputIterationStep(
420 : std::stringstream * iter_output,
421 : const GenericReal<is_ad> & effective_trial_stress,
422 : const GenericReal<is_ad> & scalar,
423 : const Real reference_residual)
424 : {
425 240136630 : if (iter_output)
426 : {
427 997128 : const unsigned int it = _iteration;
428 : const Real residual = MetaPhysicL::raw_value(_residual);
429 :
430 997128 : *iter_output << " iteration=" << it
431 997128 : << " trial_stress=" << MetaPhysicL::raw_value(effective_trial_stress)
432 1994256 : << " scalar=" << MetaPhysicL::raw_value(scalar) << " residual=" << residual
433 997128 : << " ref_res=" << reference_residual
434 997128 : << " rel_res=" << std::abs(residual) / reference_residual
435 997128 : << " rel_tol=" << _relative_tolerance << " abs_res=" << std::abs(residual)
436 997128 : << " abs_tol=" << _absolute_tolerance << '\n';
437 : }
438 240136630 : }
439 :
440 : template <bool is_ad>
441 : void
442 25364 : SingleVariableReturnMappingSolutionTempl<is_ad>::outputIterationSummary(
443 : std::stringstream * iter_output, const unsigned int total_it)
444 : {
445 25364 : if (iter_output)
446 50728 : *iter_output << "In " << total_it << " iterations the residual went from "
447 25364 : << MetaPhysicL::raw_value(_initial_residual) << " to "
448 50728 : << MetaPhysicL::raw_value(_residual) << " in '" << _svrms_name << "'."
449 : << std::endl;
450 25364 : }
451 :
452 : template class SingleVariableReturnMappingSolutionTempl<false>;
453 : template class SingleVariableReturnMappingSolutionTempl<true>;
|
