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 "GeneralizedReturnMappingSolution.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 188 : GeneralizedReturnMappingSolutionTempl<is_ad>::validParams()
29 : {
30 188 : InputParameters params = emptyInputParameters();
31 376 : params.addParam<Real>(
32 376 : "relative_tolerance", 1e-8, "Relative convergence tolerance for Newton iteration");
33 376 : params.addParam<Real>(
34 376 : "absolute_tolerance", 1e-11, "Absolute convergence tolerance for Newton iteration");
35 376 : params.addParam<Real>("acceptable_multiplier",
36 376 : 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 376 : MooseEnum internal_solve_output_on_enum("never on_error always", "on_error");
43 376 : params.addParam<MooseEnum>("internal_solve_output_on",
44 : internal_solve_output_on_enum,
45 : "When to output internal Newton solve information");
46 376 : params.addParam<bool>("internal_solve_full_iteration_history",
47 376 : 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 376 : params.addParamNamesToGroup("internal_solve_output_on internal_solve_full_iteration_history",
52 : "Debug");
53 188 : return params;
54 188 : }
55 : template <bool is_ad>
56 141 : GeneralizedReturnMappingSolutionTempl<is_ad>::GeneralizedReturnMappingSolutionTempl(
57 : const InputParameters & parameters)
58 141 : : _check_range(false),
59 141 : _line_search(true),
60 141 : _bracket_solution(false),
61 141 : _internal_solve_output_on(
62 141 : parameters.get<MooseEnum>("internal_solve_output_on").getEnum<InternalSolveOutput>()),
63 141 : _max_its(50),
64 141 : _internal_solve_full_iteration_history(
65 141 : parameters.get<bool>("internal_solve_full_iteration_history")),
66 141 : _relative_tolerance(parameters.get<Real>("relative_tolerance")),
67 141 : _absolute_tolerance(parameters.get<Real>("absolute_tolerance")),
68 141 : _acceptable_multiplier(parameters.get<Real>("acceptable_multiplier")),
69 141 : _num_resids(30),
70 141 : _residual_history(_num_resids, std::numeric_limits<Real>::max()),
71 141 : _iteration(0),
72 141 : _initial_residual(0.0),
73 141 : _residual(0.0),
74 282 : _svrms_name(parameters.get<std::string>("_object_name"))
75 : {
76 141 : }
77 :
78 : template <bool is_ad>
79 : GenericReal<is_ad>
80 0 : GeneralizedReturnMappingSolutionTempl<is_ad>::minimumPermissibleValue(
81 : const GenericDenseVector<is_ad> & /*effective_trial_stress*/) const
82 : {
83 0 : return std::numeric_limits<Real>::lowest();
84 : }
85 :
86 : template <bool is_ad>
87 : GenericReal<is_ad>
88 0 : GeneralizedReturnMappingSolutionTempl<is_ad>::maximumPermissibleValue(
89 : const GenericDenseVector<is_ad> & /*effective_trial_stress*/) const
90 : {
91 0 : return std::numeric_limits<Real>::max();
92 : }
93 :
94 : template <bool is_ad>
95 : void
96 1018592 : GeneralizedReturnMappingSolutionTempl<is_ad>::returnMappingSolve(
97 : const GenericDenseVector<is_ad> & stress_dev,
98 : const GenericDenseVector<is_ad> & stress_new,
99 : GenericReal<is_ad> & scalar,
100 : const ConsoleStream & console)
101 : {
102 : // construct the stringstream here only if the debug level is set to ALL
103 0 : std::unique_ptr<std::stringstream> iter_output =
104 1018592 : (_internal_solve_output_on == InternalSolveOutput::ALWAYS)
105 1018592 : ? std::make_unique<std::stringstream>()
106 : : nullptr;
107 :
108 : // do the internal solve and capture iteration info during the first round
109 : // iff full history output is requested regardless of whether the solve failed or succeeded
110 : auto solve_state =
111 1018592 : internalSolve(stress_dev,
112 : stress_new,
113 : scalar,
114 1018592 : _internal_solve_full_iteration_history ? iter_output.get() : nullptr);
115 1018592 : if (solve_state != SolveState::SUCCESS &&
116 0 : _internal_solve_output_on != InternalSolveOutput::ALWAYS)
117 : {
118 : // output suppressed by user, throw immediately
119 0 : if (_internal_solve_output_on == InternalSolveOutput::NEVER)
120 0 : mooseException("");
121 :
122 : // user expects some kind of output, if necessary setup output stream now
123 0 : if (!iter_output)
124 0 : iter_output = std::make_unique<std::stringstream>();
125 :
126 : // add the appropriate error message to the output
127 0 : switch (solve_state)
128 : {
129 : case SolveState::NAN_INF:
130 0 : *iter_output << "Encountered inf or nan in material return mapping iterations.\n";
131 : break;
132 :
133 : case SolveState::EXCEEDED_ITERATIONS:
134 0 : *iter_output << "Exceeded maximum iterations in material return mapping iterations.\n";
135 : break;
136 :
137 0 : default:
138 0 : mooseError("Unhandled solver state");
139 : }
140 :
141 : // if full history output is only requested for failed solves we have to repeat
142 : // the solve a second time
143 0 : if (_internal_solve_full_iteration_history)
144 0 : internalSolve(stress_dev, stress_new, scalar, iter_output.get());
145 :
146 : // Append summary and throw exception
147 0 : outputIterationSummary(iter_output.get(), _iteration);
148 0 : mooseException(iter_output->str());
149 : }
150 :
151 1018592 : if (_internal_solve_output_on == InternalSolveOutput::ALWAYS)
152 : {
153 : // the solve did not fail but the user requested debug output anyways
154 0 : outputIterationSummary(iter_output.get(), _iteration);
155 0 : console << iter_output->str();
156 : }
157 1018592 : }
158 :
159 : template <bool is_ad>
160 : typename GeneralizedReturnMappingSolutionTempl<is_ad>::SolveState
161 1018592 : GeneralizedReturnMappingSolutionTempl<is_ad>::internalSolve(
162 : const GenericDenseVector<is_ad> & stress_dev,
163 : const GenericDenseVector<is_ad> & stress_new,
164 : GenericReal<is_ad> & delta_gamma,
165 : std::stringstream * iter_output)
166 : {
167 1018592 : delta_gamma = initialGuess(stress_dev);
168 1018592 : GenericReal<is_ad> scalar_old = delta_gamma;
169 1018592 : GenericReal<is_ad> scalar_increment = 0.0;
170 1018592 : const GenericReal<is_ad> min_permissible_scalar = minimumPermissibleValue(stress_dev);
171 1018592 : const GenericReal<is_ad> max_permissible_scalar = maximumPermissibleValue(stress_dev);
172 1018592 : GenericReal<is_ad> scalar_upper_bound = max_permissible_scalar;
173 1018592 : GenericReal<is_ad> scalar_lower_bound = min_permissible_scalar;
174 1018592 : _iteration = 0;
175 :
176 1154784 : _initial_residual = _residual = computeResidual(stress_dev, stress_new, delta_gamma);
177 :
178 : Real init_resid_sign = MathUtils::sign(MetaPhysicL::raw_value(_residual));
179 1018592 : Real reference_residual =
180 1018592 : computeReferenceResidual(stress_dev, stress_new, _residual, delta_gamma);
181 :
182 : if (converged(_residual, reference_residual))
183 : {
184 45920 : iterationFinalize(delta_gamma);
185 45920 : outputIterationStep(iter_output, stress_dev, delta_gamma, reference_residual);
186 45920 : return SolveState::SUCCESS;
187 : }
188 :
189 972672 : _residual_history.assign(_num_resids, std::numeric_limits<Real>::max());
190 972672 : _residual_history[0] = MetaPhysicL::raw_value(_residual);
191 :
192 2758556 : while (_iteration < _max_its && !converged(_residual, reference_residual) &&
193 1919886 : !convergedAcceptable(_iteration, reference_residual))
194 : {
195 2758104 : scalar_increment = -_residual / computeDerivative(stress_dev, stress_new, delta_gamma);
196 1919660 : delta_gamma = scalar_old + scalar_increment;
197 :
198 1919660 : if (_check_range)
199 0 : checkPermissibleRange(delta_gamma,
200 : scalar_increment,
201 : scalar_old,
202 : min_permissible_scalar,
203 : max_permissible_scalar,
204 : iter_output);
205 :
206 1919660 : _residual = computeResidual(stress_dev, stress_new, delta_gamma);
207 :
208 1919660 : reference_residual = computeReferenceResidual(stress_dev, stress_new, _residual, delta_gamma);
209 1919660 : iterationFinalize(delta_gamma);
210 :
211 1919660 : if (_bracket_solution)
212 0 : updateBounds(delta_gamma,
213 : _residual,
214 : init_resid_sign,
215 : scalar_upper_bound,
216 : scalar_lower_bound,
217 : iter_output);
218 :
219 : if (converged(_residual, reference_residual))
220 : {
221 972446 : outputIterationStep(iter_output, stress_dev, delta_gamma, reference_residual);
222 431838 : break;
223 : }
224 : else
225 : {
226 : bool modified_increment = false;
227 :
228 947214 : if (_bracket_solution)
229 : {
230 : // Check to see whether trial scalar_increment is outside the bounds, and set it to a point
231 : // within the bounds if it is
232 0 : if (scalar_old + scalar_increment >= scalar_upper_bound ||
233 0 : scalar_old + scalar_increment <= scalar_lower_bound)
234 : {
235 0 : if (scalar_upper_bound != max_permissible_scalar &&
236 0 : scalar_lower_bound != min_permissible_scalar)
237 : {
238 0 : const Real frac = 0.5;
239 0 : scalar_increment =
240 0 : (1.0 - frac) * scalar_lower_bound + frac * scalar_upper_bound - scalar_old;
241 : modified_increment = true;
242 0 : if (iter_output)
243 0 : *iter_output << " Trial scalar_increment exceeded bounds. Setting between "
244 : "lower/upper bounds. frac: "
245 : << frac << std::endl;
246 : }
247 : }
248 : }
249 :
250 : // Update the trial scalar and recompute residual if the line search or bounds checking
251 : // modified the increment
252 : if (modified_increment)
253 : {
254 0 : delta_gamma = scalar_old + scalar_increment;
255 0 : _residual = computeResidual(stress_dev, stress_new, delta_gamma);
256 0 : reference_residual =
257 0 : computeReferenceResidual(stress_dev, stress_new, _residual, delta_gamma);
258 0 : iterationFinalize(delta_gamma);
259 :
260 0 : if (_bracket_solution)
261 0 : updateBounds(delta_gamma,
262 : _residual,
263 : init_resid_sign,
264 : scalar_upper_bound,
265 : scalar_lower_bound,
266 : iter_output);
267 : }
268 : }
269 :
270 947214 : outputIterationStep(iter_output, stress_dev, delta_gamma, reference_residual);
271 :
272 947214 : ++_iteration;
273 947214 : scalar_old = delta_gamma;
274 947214 : _residual_history[_iteration % _num_resids] = MetaPhysicL::raw_value(_residual);
275 : }
276 :
277 972672 : if (std::isnan(_residual) || std::isinf(MetaPhysicL::raw_value(_residual)))
278 : return SolveState::NAN_INF;
279 :
280 972672 : if (_iteration == _max_its)
281 : return SolveState::EXCEEDED_ITERATIONS;
282 :
283 : return SolveState::SUCCESS;
284 : }
285 :
286 : template <bool is_ad>
287 : bool
288 0 : GeneralizedReturnMappingSolutionTempl<is_ad>::converged(const GenericReal<is_ad> & ad_residual,
289 : const Real & reference)
290 : {
291 : const Real residual = MetaPhysicL::raw_value(ad_residual);
292 4858364 : return (std::abs(residual) <= _absolute_tolerance ||
293 3853167 : std::abs(residual / reference) <= _relative_tolerance);
294 : }
295 :
296 : template <bool is_ad>
297 : bool
298 1919886 : GeneralizedReturnMappingSolutionTempl<is_ad>::convergedAcceptable(const unsigned int it,
299 : const Real & reference)
300 : {
301 : // Require that we have at least done _num_resids evaluations before we allow for
302 : // acceptable convergence
303 1919886 : if (it < _num_resids)
304 : return false;
305 :
306 : // Check to see whether the residual has dropped by convergence_history_factor over
307 : // the last _num_resids iterations. If it has (which means it's still making progress),
308 : // don't consider it to be converged within the acceptable limits.
309 452 : const Real convergence_history_factor = 10.0;
310 904 : if (std::abs(_residual * convergence_history_factor) <
311 452 : std::abs(_residual_history[(it + 1) % _num_resids]))
312 : return false;
313 :
314 : // Now that it's determined that progress is not being made, treat it as converged if
315 : // we're within the acceptable convergence limits
316 452 : return converged(_residual / _acceptable_multiplier, reference);
317 : }
318 :
319 : template <bool is_ad>
320 : void
321 0 : GeneralizedReturnMappingSolutionTempl<is_ad>::checkPermissibleRange(
322 : GenericReal<is_ad> & scalar,
323 : GenericReal<is_ad> & scalar_increment,
324 : const GenericReal<is_ad> & scalar_old,
325 : const GenericReal<is_ad> min_permissible_scalar,
326 : const GenericReal<is_ad> max_permissible_scalar,
327 : std::stringstream * iter_output)
328 : {
329 0 : if (scalar > max_permissible_scalar)
330 : {
331 0 : scalar_increment = (max_permissible_scalar - scalar_old) / 2.0;
332 0 : scalar = scalar_old + scalar_increment;
333 0 : if (iter_output)
334 0 : *iter_output << "Scalar greater than maximum ("
335 : << MetaPhysicL::raw_value(max_permissible_scalar)
336 0 : << ") adjusted scalar=" << MetaPhysicL::raw_value(scalar)
337 0 : << " scalar_increment=" << MetaPhysicL::raw_value(scalar_increment) << std::endl;
338 : }
339 0 : else if (scalar < min_permissible_scalar)
340 : {
341 0 : scalar_increment = (min_permissible_scalar - scalar_old) / 2.0;
342 0 : scalar = scalar_old + scalar_increment;
343 0 : if (iter_output)
344 0 : *iter_output << "Scalar less than minimum (" << MetaPhysicL::raw_value(min_permissible_scalar)
345 0 : << ") adjusted scalar=" << MetaPhysicL::raw_value(scalar)
346 0 : << " scalar_increment=" << MetaPhysicL::raw_value(scalar_increment) << std::endl;
347 : }
348 0 : }
349 :
350 : template <bool is_ad>
351 : void
352 0 : GeneralizedReturnMappingSolutionTempl<is_ad>::updateBounds(const GenericReal<is_ad> & scalar,
353 : const GenericReal<is_ad> & residual,
354 : const Real init_resid_sign,
355 : GenericReal<is_ad> & scalar_upper_bound,
356 : GenericReal<is_ad> & scalar_lower_bound,
357 : std::stringstream * iter_output)
358 : {
359 : // Update upper/lower bounds as applicable
360 0 : if (residual * init_resid_sign < 0.0 && scalar < scalar_upper_bound)
361 : {
362 0 : scalar_upper_bound = scalar;
363 0 : if (scalar_upper_bound < scalar_lower_bound)
364 : {
365 0 : scalar_upper_bound = scalar_lower_bound;
366 0 : scalar_lower_bound = 0.0;
367 0 : if (iter_output)
368 0 : *iter_output << " Corrected for scalar_upper_bound < scalar_lower_bound" << std::endl;
369 : }
370 : }
371 : // Don't permit setting scalar_lower_bound > scalar_upper_bound (but do permit the reverse).
372 : // This ensures that if we encounter multiple roots, we pick the lowest one.
373 0 : else if (residual * init_resid_sign > 0.0 && scalar > scalar_lower_bound &&
374 0 : scalar < scalar_upper_bound)
375 0 : scalar_lower_bound = scalar;
376 0 : }
377 :
378 : template <bool is_ad>
379 : void
380 1965580 : GeneralizedReturnMappingSolutionTempl<is_ad>::outputIterationStep(
381 : std::stringstream * iter_output,
382 : const GenericDenseVector<is_ad> & effective_trial_stress,
383 : const GenericReal<is_ad> & scalar,
384 : const GenericReal<is_ad> reference_residual)
385 : {
386 1965580 : if (iter_output)
387 : {
388 0 : const unsigned int it = _iteration;
389 : const Real residual = MetaPhysicL::raw_value(_residual);
390 :
391 0 : *iter_output << " iteration=" << it
392 0 : << " trial_stress=" << MetaPhysicL::raw_value(effective_trial_stress)
393 0 : << " scalar=" << MetaPhysicL::raw_value(scalar) << " residual=" << residual
394 0 : << " ref_res=" << reference_residual
395 0 : << " rel_res=" << std::abs(residual) / reference_residual
396 0 : << " rel_tol=" << _relative_tolerance << " abs_res=" << std::abs(residual)
397 0 : << " abs_tol=" << _absolute_tolerance << '\n';
398 : }
399 1965580 : }
400 :
401 : template <bool is_ad>
402 : void
403 0 : GeneralizedReturnMappingSolutionTempl<is_ad>::outputIterationSummary(
404 : std::stringstream * iter_output, const unsigned int total_it)
405 : {
406 0 : if (iter_output)
407 0 : *iter_output << "In " << total_it << " iterations the residual went from "
408 0 : << MetaPhysicL::raw_value(_initial_residual) << " to "
409 0 : << MetaPhysicL::raw_value(_residual) << " in '" << _svrms_name << "'.\n";
410 0 : }
411 :
412 : template class GeneralizedReturnMappingSolutionTempl<false>;
413 : template class GeneralizedReturnMappingSolutionTempl<true>;
|
