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