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 "LAROMANCEStressUpdateBase.h"
11 : #include "Function.h"
12 : #include "MathUtils.h"
13 : #include "MooseUtils.h"
14 : #include "MooseRandom.h"
15 : #include "Units.h"
16 :
17 : #include <iostream>
18 :
19 : registerMooseObjectAliased("TensorMechanicsApp",
20 : LAROMANCEStressUpdateBase,
21 : "LAROMANCEStressUpdate");
22 : registerMooseObjectAliased("TensorMechanicsApp",
23 : ADLAROMANCEStressUpdateBase,
24 : "ADLAROMANCEStressUpdate");
25 :
26 : template <bool is_ad>
27 : InputParameters
28 300 : LAROMANCEStressUpdateBaseTempl<is_ad>::validParams()
29 : {
30 300 : InputParameters params = RadialReturnCreepStressUpdateBaseTempl<is_ad>::validParams();
31 300 : params.addClassDescription("Base class to calculate the effective creep strain based on the "
32 : "rates predicted by a material specific Los Alamos Reduced Order "
33 : "Model derived from a Visco-Plastic Self Consistent calculations.");
34 :
35 600 : params.addRequiredCoupledVar("temperature", "The coupled temperature (K)");
36 600 : params.addParam<MaterialPropertyName>("environmental_factor",
37 : "Optional coupled environmental factor");
38 :
39 600 : MooseEnum error_lower_limit_behavior("ERROR EXCEPTION WARN IGNORE DONTHING USELIMIT",
40 : "EXCEPTION");
41 : // Only allow ERROR and EXCEPTION on upper bounds
42 600 : MooseEnum error_upper_limit_behavior("ERROR EXCEPTION", "EXCEPTION");
43 600 : params.addParam<MooseEnum>(
44 : "cell_input_window_low_failure",
45 : error_lower_limit_behavior,
46 : "What to do if cell dislocation concentration is outside the lower global "
47 : "window of applicability.");
48 600 : params.addParam<MooseEnum>(
49 : "cell_input_window_high_failure",
50 : error_upper_limit_behavior,
51 : "What to do if cell dislocation concentration is outside the upper global "
52 : "window of applicability.");
53 600 : params.addParam<MooseEnum>("wall_input_window_low_failure",
54 : error_lower_limit_behavior,
55 : "What to do if wall dislocation concentration is outside the "
56 : "lower global window of applicability.");
57 600 : params.addParam<MooseEnum>("wall_input_window_high_failure",
58 : error_upper_limit_behavior,
59 : "What to do if wall dislocation concentration is outside the "
60 : "upper global window of applicability.");
61 600 : params.addParam<MooseEnum>(
62 : "old_strain_input_window_low_failure",
63 : error_lower_limit_behavior,
64 : "What to do if old strain is outside the lower global window of applicability.");
65 600 : params.addParam<MooseEnum>(
66 : "old_strain_input_window_high_failure",
67 : error_upper_limit_behavior,
68 : "What to do if old strain is outside the upper global window of applicability.");
69 :
70 600 : MooseEnum extrapolated_lower_limit_behavior(
71 : "ERROR EXCEPTION WARN IGNORE DONOTHING USELIMIT EXTRAPOLATE", "EXTRAPOLATE");
72 600 : params.addParam<MooseEnum>(
73 : "stress_input_window_low_failure",
74 : extrapolated_lower_limit_behavior,
75 : "What to do if stress is outside the lower global window of applicability.");
76 600 : params.addParam<MooseEnum>(
77 : "stress_input_window_high_failure",
78 : error_upper_limit_behavior,
79 : "What to do if stress is outside the upper global window of applicability.");
80 600 : params.addParam<MooseEnum>(
81 : "temperature_input_window_low_failure",
82 : extrapolated_lower_limit_behavior,
83 : "What to do if temperature is outside the lower global window of applicability.");
84 600 : params.addParam<MooseEnum>(
85 : "temperature_input_window_high_failure",
86 : error_upper_limit_behavior,
87 : "What to do if temperature is outside the upper global window of applicability.");
88 600 : params.addParam<MooseEnum>(
89 : "environment_input_window_low_failure",
90 : extrapolated_lower_limit_behavior,
91 : "What to do if environmental factor is outside the lower global window of applicability.");
92 600 : params.addParam<MooseEnum>(
93 : "environment_input_window_high_failure",
94 : error_upper_limit_behavior,
95 : "What to do if environmental factor is outside the upper global window of applicability.");
96 :
97 600 : params.addRequiredRangeCheckedParam<Real>(
98 : "initial_cell_dislocation_density",
99 : "initial_cell_dislocation_density >= 0.0",
100 : "Initial density of cell (glissile) dislocations (1/m^2)");
101 900 : params.addRangeCheckedParam<Real>(
102 : "cell_dislocations_normal_distribution_width",
103 600 : 0.0,
104 : "cell_dislocations_normal_distribution_width >= 0.0",
105 : "Width of the normal distribution to assign to the initial cell dislocation value. This is "
106 : "given as a fraction of the initial_cell_dislocation_density.");
107 900 : params.addRangeCheckedParam<Real>(
108 : "max_relative_cell_dislocation_increment",
109 600 : 0.5,
110 : "max_relative_cell_dislocation_increment > 0.0",
111 : "Maximum increment of density of cell (glissile) dislocations.");
112 :
113 600 : params.addRequiredRangeCheckedParam<Real>(
114 : "initial_wall_dislocation_density",
115 : "initial_wall_dislocation_density >= 0.0",
116 : "wall (locked) dislocation density initial value (1/m^2).");
117 900 : params.addRangeCheckedParam<Real>(
118 : "wall_dislocations_normal_distribution_width",
119 600 : 0.0,
120 : "wall_dislocations_normal_distribution_width >= 0.0",
121 : "Width of the normal distribution to assign to the initial wall dislocation value. This is "
122 : "given as a fraction of the initial_wall_dislocation_density.");
123 900 : params.addRangeCheckedParam<Real>(
124 : "max_relative_wall_dislocation_increment",
125 600 : 0.5,
126 : "max_relative_wall_dislocation_increment > 0.0",
127 : "Maximum increment of wall (locked) dislocation density initial value (1/m^2).");
128 :
129 600 : params.addParam<bool>("verbose", false, "Flag to output verbose information.");
130 :
131 600 : params.addParam<FunctionName>(
132 : "cell_dislocation_density_forcing_function",
133 : "Optional forcing function for wall dislocation. If provided, the wall dislocation "
134 : "density will be reset to the function value at the beginning of the timestep. Used for "
135 : "testing purposes only.");
136 600 : params.addParam<FunctionName>(
137 : "wall_dislocation_density_forcing_function",
138 : "Optional forcing function for wall dislocation. If provided, the wall dislocation "
139 : "density will be reset to the function value at the beginning of the timestep. Used for "
140 : "testing purposes only.");
141 600 : params.addParam<FunctionName>(
142 : "old_creep_strain_forcing_function",
143 : "Optional forcing function for the creep strain from the previous timestep. If provided, "
144 : "the old creep strain will be reset to the function value at the beginning of the "
145 : "timestep. Used for testing purposes only.");
146 600 : params.addParam<FunctionName>(
147 : "effective_stress_forcing_function",
148 : "Optional forcing function for the effective stress. If provided, the effective stress will "
149 : "be reset to the function value at the beginning of the timestep. Used for testing purposes "
150 : "only.");
151 :
152 600 : params.addParam<unsigned int>("seed", 0, "Random number generator seed");
153 600 : params.addParam<std::string>("stress_unit", "Pa", "unit of stress");
154 :
155 : // use std::string here to avoid automatic absolute path expansion
156 600 : params.addParam<FileName>("model", "LaRomance model JSON datafile");
157 600 : params.addParam<FileName>("export_model", "Write LaRomance model to JSON datafile");
158 :
159 600 : params.addParamNamesToGroup(
160 : "cell_dislocation_density_forcing_function wall_dislocation_density_forcing_function "
161 : "old_creep_strain_forcing_function effective_stress_forcing_function seed stress_unit",
162 : "Advanced");
163 :
164 300 : return params;
165 300 : }
166 :
167 : template <bool is_ad>
168 225 : LAROMANCEStressUpdateBaseTempl<is_ad>::LAROMANCEStressUpdateBaseTempl(
169 : const InputParameters & parameters)
170 : : RadialReturnCreepStressUpdateBaseTempl<is_ad>(parameters),
171 225 : _temperature(this->template coupledGenericValue<is_ad>("temperature")),
172 225 : _environmental(
173 225 : this->isParamValid("environmental_factor")
174 225 : ? &this->template getGenericMaterialProperty<Real, is_ad>("environmental_factor")
175 : : nullptr),
176 450 : _window_failure(_environmental ? 6 : 5),
177 :
178 450 : _verbose(this->template getParam<bool>("verbose")),
179 225 : _cell_dislocations(
180 225 : this->template declareGenericProperty<Real, is_ad>(this->_base_name + "cell_dislocations")),
181 225 : _cell_dislocations_old(
182 225 : this->template getMaterialPropertyOld<Real>(this->_base_name + "cell_dislocations")),
183 450 : _max_cell_increment(this->template getParam<Real>("max_relative_cell_dislocation_increment")),
184 225 : _cell_function(this->isParamValid("cell_dislocation_density_forcing_function")
185 273 : ? &this->getFunction("cell_dislocation_density_forcing_function")
186 : : NULL),
187 225 : _cell_dislocation_increment(0.0),
188 225 : _wall_dislocations(
189 225 : this->template declareGenericProperty<Real, is_ad>(this->_base_name + "wall_dislocations")),
190 225 : _wall_dislocations_old(
191 225 : this->template getMaterialPropertyOld<Real>(this->_base_name + "wall_dislocations")),
192 450 : _max_wall_increment(this->template getParam<Real>("max_relative_wall_dislocation_increment")),
193 225 : _wall_function(this->isParamValid("wall_dislocation_density_forcing_function")
194 273 : ? &this->getFunction("wall_dislocation_density_forcing_function")
195 : : NULL),
196 225 : _stress_function(this->isParamValid("effective_stress_forcing_function")
197 273 : ? &this->getFunction("effective_stress_forcing_function")
198 : : NULL),
199 225 : _wall_dislocation_increment(0.0),
200 225 : _cell_input_index(0),
201 225 : _wall_input_index(1),
202 225 : _stress_input_index(2),
203 225 : _old_strain_input_index(3),
204 225 : _temperature_input_index(4),
205 225 : _environmental_input_index(5),
206 225 : _cell_output_index(0),
207 225 : _wall_output_index(1),
208 225 : _strain_output_index(2),
209 :
210 450 : _creep_strain_old_forcing_function(this->isParamValid("old_creep_strain_forcing_function")
211 273 : ? &this->getFunction("old_creep_strain_forcing_function")
212 : : NULL),
213 :
214 225 : _creep_rate(
215 450 : this->template declareGenericProperty<Real, is_ad>(this->_base_name + "creep_rate")),
216 450 : _cell_rate(this->template declareGenericProperty<Real, is_ad>(this->_base_name +
217 : "cell_dislocation_rate")),
218 450 : _wall_rate(this->template declareGenericProperty<Real, is_ad>(this->_base_name +
219 : "wall_dislocation_rate")),
220 225 : _extrapolation(this->template declareProperty<Real>("ROM_extrapolation")),
221 225 : _second_partition_weight(
222 225 : this->template declareGenericProperty<Real, is_ad>("partition_weight")),
223 :
224 225 : _derivative(0.0),
225 225 : _old_input_values(3),
226 450 : _wall_dislocations_step(this->template declareGenericProperty<Real, is_ad>(
227 : this->_base_name + "wall_dislocations_step")),
228 450 : _cell_dislocations_step(this->template declareGenericProperty<Real, is_ad>(
229 : this->_base_name + "cell_dislocations_step")),
230 225 : _plastic_strain_increment(),
231 225 : _number_of_substeps(
232 225 : this->template declareProperty<Real>(this->_base_name + "number_of_substeps")),
233 900 : _index_name(_window_failure.size())
234 : {
235 225 : this->_check_range = true; // this may not be necessary?
236 :
237 225 : _index_name[_cell_input_index] = "cell";
238 225 : _index_name[_wall_input_index] = "wall";
239 225 : _index_name[_stress_input_index] = "stress";
240 225 : _index_name[_old_strain_input_index] = "old strain";
241 225 : _index_name[_temperature_input_index] = "temperature";
242 :
243 225 : _window_failure[_cell_input_index].first =
244 225 : parameters.get<MooseEnum>("cell_input_window_low_failure").getEnum<WindowFailure>();
245 225 : _window_failure[_cell_input_index].second =
246 225 : parameters.get<MooseEnum>("cell_input_window_high_failure").getEnum<WindowFailure>();
247 225 : _window_failure[_wall_input_index].first =
248 225 : parameters.get<MooseEnum>("wall_input_window_low_failure").getEnum<WindowFailure>();
249 225 : _window_failure[_wall_input_index].second =
250 225 : parameters.get<MooseEnum>("wall_input_window_high_failure").getEnum<WindowFailure>();
251 225 : _window_failure[_stress_input_index].first =
252 225 : parameters.get<MooseEnum>("stress_input_window_low_failure").getEnum<WindowFailure>();
253 225 : _window_failure[_stress_input_index].second =
254 225 : parameters.get<MooseEnum>("stress_input_window_high_failure").getEnum<WindowFailure>();
255 225 : _window_failure[_old_strain_input_index].first =
256 225 : parameters.get<MooseEnum>("old_strain_input_window_low_failure").getEnum<WindowFailure>();
257 225 : _window_failure[_old_strain_input_index].second =
258 225 : parameters.get<MooseEnum>("old_strain_input_window_high_failure").getEnum<WindowFailure>();
259 225 : _window_failure[_temperature_input_index].first =
260 225 : parameters.get<MooseEnum>("temperature_input_window_low_failure").getEnum<WindowFailure>();
261 225 : _window_failure[_temperature_input_index].second =
262 225 : parameters.get<MooseEnum>("temperature_input_window_high_failure").getEnum<WindowFailure>();
263 225 : if (_environmental)
264 : {
265 0 : _index_name[_environmental_input_index] = "environmental";
266 0 : _window_failure[_environmental_input_index].first =
267 0 : parameters.get<MooseEnum>("environment_input_window_low_failure").getEnum<WindowFailure>();
268 0 : _window_failure[_environmental_input_index].second =
269 0 : parameters.get<MooseEnum>("environment_input_window_high_failure").getEnum<WindowFailure>();
270 : }
271 :
272 : // load JSON datafile
273 450 : if (this->isParamValid("model"))
274 : {
275 54 : const auto model_file_name = this->getDataFileName("model");
276 27 : std::ifstream model_file(model_file_name.c_str());
277 27 : model_file >> _json;
278 27 : }
279 :
280 225 : setupUnitConversionFactors(parameters);
281 225 : }
282 :
283 : template <bool is_ad>
284 : void
285 9 : LAROMANCEStressUpdateBaseTempl<is_ad>::exportJSON()
286 : {
287 9 : _json["strain_cutoff"] = getStrainCutoff();
288 18 : _json["transform"] = getTransform();
289 18 : _json["transform_coefs"] = getTransformCoefs();
290 18 : _json["input_limits"] = getInputLimits();
291 18 : _json["normalization_limits"] = getNormalizationLimits();
292 18 : _json["coefs"] = getCoefs();
293 18 : _json["tiling"] = getTilings();
294 9 : _json["cutoff"] = getStrainCutoff();
295 9 : }
296 :
297 : template <bool is_ad>
298 : bool
299 146734 : LAROMANCEStressUpdateBaseTempl<is_ad>::substeppingCapabilityEnabled()
300 : {
301 146734 : return this->_use_substepping != RadialReturnStressUpdateTempl<is_ad>::SubsteppingType::NONE;
302 : }
303 :
304 : template <bool is_ad>
305 : void
306 225 : LAROMANCEStressUpdateBaseTempl<is_ad>::setupUnitConversionFactors(
307 : const InputParameters & parameters)
308 : {
309 : // Stress unit conversion factor
310 450 : const MooseUnits stress_unit_to("MPa");
311 225 : const MooseUnits stress_unit_from(parameters.get<std::string>("stress_unit"));
312 225 : _stress_ucf = stress_unit_to.convert(1, stress_unit_from);
313 225 : }
314 :
315 : template <bool is_ad>
316 : void
317 225 : LAROMANCEStressUpdateBaseTempl<is_ad>::initialSetup()
318 : {
319 : // export models that are compiled in
320 450 : if (this->isParamValid("export_model"))
321 : {
322 9 : exportJSON();
323 18 : std::ofstream out(this->template getParam<FileName>("export_model").c_str());
324 9 : out << _json;
325 9 : }
326 :
327 : // Pull in relevant ROM information and do sanity checks
328 225 : _transform = getTransform();
329 225 : _transform_coefs = getTransformCoefs();
330 225 : _input_limits = getInputLimits();
331 225 : _normalization_limits = getNormalizationLimits();
332 225 : _coefs = getCoefs();
333 225 : _tiling = getTilings();
334 450 : _cutoff = getStrainCutoff();
335 : // resize containers to be filled later based on partition dimension
336 : // and immediately run some sanity checks
337 225 : _num_partitions = _transform.size();
338 225 : if (_num_partitions < 1 || _num_partitions > 2)
339 0 : mooseError(
340 0 : "In ", _name, ": First dimension of getTransform must be either size one or size two");
341 225 : if (_transform[0].size() < 1)
342 0 : mooseError("In ", _name, ": Transform is not the correct shape");
343 225 : _num_outputs = _transform[0][0].size();
344 225 : if (_num_outputs != 3)
345 0 : mooseError("In ",
346 0 : _name,
347 : ": ",
348 0 : _num_outputs,
349 : " outputs detected. Three and only three outputs are currently supported.");
350 :
351 225 : _num_inputs = _transform[0][0][0].size();
352 225 : if (_num_inputs != 5 && _num_inputs != 6)
353 0 : mooseError("In ",
354 0 : _name,
355 : ": ",
356 0 : _num_inputs,
357 : " inputs detected. Only five or six inputs currently supported.");
358 225 : if (_num_inputs == 6 && !_environmental)
359 0 : this->template paramError(
360 : "environmental_factor",
361 : "Number of ROM inputs indicate environmental factor is required to be coupled.");
362 225 : if (_num_inputs != 6 && _environmental)
363 0 : this->template paramError(
364 : "environmental_factor",
365 : "Number of ROM inputs indicate environmental factor is not implemented, but "
366 : "environmental factor coupled.");
367 225 : _num_tiles.resize(_num_partitions);
368 225 : _num_coefs.resize(_num_partitions);
369 225 : _degree.resize(_num_partitions);
370 225 : _precomputed_vals.resize(_num_partitions);
371 225 : _rom_inputs.resize(_num_partitions);
372 225 : _polynomial_inputs.resize(_num_partitions);
373 225 : _non_stress_weights.resize(_num_partitions);
374 225 : _weights.resize(_num_partitions);
375 225 : _partition_weights.resize(_num_partitions);
376 225 : _dpartition_weight_dstress.resize(_num_partitions);
377 225 : _transformed_normalization_limits.resize(_num_partitions);
378 225 : _makeframe_helper.resize(_num_partitions);
379 225 : _global_limits.resize(_num_inputs);
380 : // temporarily fill global limits with extreme numerical values, to later update
381 1350 : for (unsigned int i = 0; i < _num_inputs; ++i)
382 1125 : _global_limits[i] = {std::numeric_limits<Real>::max(), 0.0};
383 :
384 : // start loop over partitions to perform sanity checks, set global limits,
385 : // and print global configurations
386 225 : if (_transform.size() != _num_partitions || _transform_coefs.size() != _num_partitions ||
387 225 : _input_limits.size() != _num_partitions || _normalization_limits.size() != _num_partitions ||
388 450 : _coefs.size() != _num_partitions || _tiling.size() != _num_partitions ||
389 : _cutoff.size() != _num_partitions)
390 0 : mooseError(
391 0 : "In ", _name, ": one of the ROM inputs does not have the correct number of partitions");
392 :
393 504 : for (unsigned int p = 0; p < _num_partitions; ++p)
394 : {
395 279 : _num_tiles[p] = _transform[p].size();
396 279 : if (!_num_tiles[p])
397 0 : mooseError("In ", _name, ": No tiles detected. Double check your ROM input");
398 :
399 : bool correct_shape = true;
400 279 : if (_transform[p].size() != _num_tiles[p] || _transform_coefs[p].size() != _num_tiles[p] ||
401 279 : _input_limits[p].size() != _num_tiles[p] ||
402 558 : _normalization_limits[p].size() != _num_tiles[p] || _coefs[p].size() != _num_tiles[p])
403 : correct_shape = false;
404 279 : if (_tiling[p].size() != _num_inputs)
405 : correct_shape = false;
406 279 : if (_coefs[p][0].size() == 0)
407 : correct_shape = false;
408 279 : _num_coefs[p] = _coefs[p][0][0].size();
409 :
410 : // start loop over tiles to perform sanity checks.
411 828 : for (unsigned int t = 0; t < _num_tiles[p]; ++t)
412 : {
413 549 : if (_transform[p][t].size() != _num_outputs ||
414 1098 : _transform_coefs[p][t].size() != _num_outputs || _coefs[p][t].size() != _num_outputs)
415 : correct_shape = false;
416 2196 : for (unsigned int o = 0; o < _num_outputs; ++o)
417 1647 : if (_transform[p][t][o].size() != _num_inputs ||
418 3294 : _transform_coefs[p][t][o].size() != _num_inputs ||
419 : _coefs[p][t][o].size() != _num_coefs[p])
420 : correct_shape = false;
421 549 : if (_input_limits[p][t].size() != _num_inputs ||
422 : _normalization_limits[p][t].size() != _num_inputs)
423 : correct_shape = false;
424 3294 : for (unsigned int i = 0; i < _num_inputs; ++i)
425 2745 : if (_input_limits[p][t][i].size() != 2 || _normalization_limits[p][t][i].size() != 2)
426 : correct_shape = false;
427 : }
428 :
429 279 : if (!correct_shape)
430 0 : mooseError("In ", _name, ": ROM data is not the right shape.");
431 :
432 279 : _degree[p] = std::pow(_num_coefs[p], 1.0 / _num_inputs);
433 279 : if (!_degree[p] || _degree[p] > 4)
434 0 : mooseError("In ", _name, ": degree must be 1, 2, 3 or 4.");
435 :
436 : // Check input limits and find global limits of all partitions. Note that this will return the
437 : // extremes of the model! If the model is not flat along one input limit, then global limits
438 : // will not truely capture the mutli-dimensionality of the problem. Consequently, the user may
439 : // find input values that result in errors in the partition weight computation.
440 1674 : for (unsigned int i = 0; i < _num_inputs; ++i)
441 : {
442 4140 : for (unsigned int t = 0; t < _num_tiles[p]; ++t)
443 : {
444 2745 : if (_input_limits[p][t][i][0] >= _input_limits[p][t][i][1])
445 0 : mooseError("In ", _name, ": Input limits are ordered incorrectly");
446 3978 : _global_limits[i].first = std::min(_global_limits[i].first, _input_limits[p][t][i][0]);
447 3978 : _global_limits[i].second = std::max(_global_limits[i].second, _input_limits[p][t][i][1]);
448 : }
449 : }
450 :
451 : // Precompute helper containers
452 558 : _transformed_normalization_limits[p] = getTransformedLimits(p, _normalization_limits[p]);
453 558 : _makeframe_helper[p] = getMakeFrameHelper(p);
454 :
455 : // Prepare containers for each partition
456 558 : _precomputed_vals[p].resize(_num_tiles[p], std::vector<GenericReal<is_ad>>(_num_coefs[p]));
457 558 : _rom_inputs[p].resize(_num_tiles[p], std::vector<GenericReal<is_ad>>(_num_inputs));
458 279 : _polynomial_inputs[p].resize(_num_tiles[p],
459 279 : std::vector<std::vector<GenericReal<is_ad>>>(
460 558 : _num_inputs, std::vector<GenericReal<is_ad>>(_degree[p])));
461 279 : _non_stress_weights[p].resize(_num_tiles[p]);
462 279 : _weights[p].resize(_num_tiles[p], 0);
463 : }
464 : // Prepare containers independent of partition
465 225 : _input_values.resize(_num_inputs);
466 :
467 225 : if (_verbose)
468 : {
469 36 : Moose::err << "ROM model info: " << _name << "\n";
470 : Moose::err << " number of tiles, partition 1: " << _num_tiles[0] << "\n";
471 36 : if (_num_partitions > 1)
472 : Moose::err << " number of tiles, partition 2: " << _num_tiles[1] << "\n";
473 : Moose::err << " number of outputs: " << _num_outputs << "\n";
474 : Moose::err << " number of inputs: " << _num_inputs << "\n";
475 : Moose::err << " degree (max Legendre degree + constant), partition 1: " << _degree[0] << "\n";
476 36 : if (_num_partitions > 1)
477 : Moose::err << " degree (max Legendre degree + constant), partition 2: " << _degree[1] << "\n";
478 : Moose::err << " number of coefficients, partition 1: " << _num_coefs[0] << "\n";
479 36 : if (_num_partitions > 1)
480 : Moose::err << " number of coefficients, partition 2: " << _num_coefs[1] << "\n";
481 : Moose::err << " Global limits:\n cell dislocations ("
482 36 : << _global_limits[_cell_input_index].first << " - "
483 36 : << _global_limits[_cell_input_index].second << ")\n";
484 36 : Moose::err << " wall dislocations (" << _global_limits[_wall_input_index].first << " - "
485 36 : << _global_limits[_wall_input_index].second << ")\n";
486 36 : Moose::err << " Stress (" << _global_limits[_stress_input_index].first << " - "
487 36 : << _global_limits[_stress_input_index].second << ")\n";
488 36 : Moose::err << " Old strain (" << _global_limits[_old_strain_input_index].first << " - "
489 36 : << _global_limits[_old_strain_input_index].second << ")\n";
490 36 : Moose::err << " Temperature (" << _global_limits[_temperature_input_index].first << " - "
491 36 : << _global_limits[_temperature_input_index].second << ")\n";
492 36 : if (_environmental)
493 0 : Moose::err << " Environmental factor (" << _global_limits[_environmental_input_index].first
494 0 : << " - " << _global_limits[_environmental_input_index].second << ")\n";
495 : Moose::err << std::endl;
496 : }
497 225 : }
498 :
499 : template <bool is_ad>
500 : void
501 3560 : LAROMANCEStressUpdateBaseTempl<is_ad>::initQpStatefulProperties()
502 : {
503 : MooseRandom rng;
504 7120 : rng.seed(0, this->template getParam<unsigned int>("seed"));
505 :
506 3560 : _cell_dislocations[_qp] = rng.randNormal(
507 7120 : this->template getParam<Real>("initial_cell_dislocation_density"),
508 3560 : this->template getParam<Real>("initial_cell_dislocation_density") *
509 7120 : this->template getParam<Real>("cell_dislocations_normal_distribution_width"));
510 3560 : _wall_dislocations[_qp] = rng.randNormal(
511 7120 : this->template getParam<Real>("initial_wall_dislocation_density"),
512 3560 : this->template getParam<Real>("initial_wall_dislocation_density") *
513 7120 : this->template getParam<Real>("wall_dislocations_normal_distribution_width"));
514 :
515 3560 : RadialReturnCreepStressUpdateBaseTempl<is_ad>::initQpStatefulProperties();
516 3560 : }
517 :
518 : template <bool is_ad>
519 : GenericReal<is_ad>
520 143872 : LAROMANCEStressUpdateBaseTempl<is_ad>::maximumPermissibleValue(
521 : const GenericReal<is_ad> & effective_trial_stress) const
522 : {
523 : // Make maximum allowed scalar a little bit less than the deformation that would reduce the
524 : // trial stress to zero. This prevents negative trial stresses.
525 143872 : return effective_trial_stress / this->_three_shear_modulus * 0.999999;
526 : }
527 :
528 : template <bool is_ad>
529 : void
530 146716 : LAROMANCEStressUpdateBaseTempl<is_ad>::resetIncrementalMaterialProperties()
531 : {
532 146716 : _cell_dislocation_increment = 0.0;
533 146716 : _wall_dislocation_increment = 0.0;
534 :
535 : _plastic_strain_increment.zero();
536 :
537 146716 : _wall_dislocations_step[_qp] = 0.0;
538 146716 : _cell_dislocations_step[_qp] = 0.0;
539 146716 : }
540 :
541 : template <bool is_ad>
542 : void
543 30960 : LAROMANCEStressUpdateBaseTempl<is_ad>::storeIncrementalMaterialProperties(
544 : const unsigned int total_number_of_substeps)
545 : {
546 30960 : _wall_dislocations_step[_qp] += _wall_dislocation_increment;
547 30960 : _cell_dislocations_step[_qp] += _cell_dislocation_increment;
548 30960 : _number_of_substeps[_qp] = total_number_of_substeps;
549 30960 : }
550 :
551 : template <bool is_ad>
552 : void
553 164738 : LAROMANCEStressUpdateBaseTempl<is_ad>::computeStressInitialize(
554 : const GenericReal<is_ad> & effective_trial_stress,
555 : const GenericRankFourTensor<is_ad> & elasticity_tensor)
556 : {
557 67168 : RadialReturnCreepStressUpdateBaseTempl<is_ad>::computeStressInitialize(effective_trial_stress,
558 : elasticity_tensor);
559 : // Previous substep creep strain
560 164738 : RankTwoTensor creep_strain_substep = this->_creep_strain_old[_qp] + _plastic_strain_increment;
561 :
562 : // Prepare old values
563 164738 : _old_input_values[_cell_output_index] =
564 164738 : _cell_function
565 164738 : ? _cell_function->value(_t, _q_point[_qp])
566 158242 : : (_cell_dislocations_old[_qp] + MetaPhysicL::raw_value(_cell_dislocations_step[_qp]));
567 164738 : _old_input_values[_wall_output_index] =
568 164738 : _wall_function
569 164738 : ? _wall_function->value(_t, _q_point[_qp])
570 158242 : : (_wall_dislocations_old[_qp] + MetaPhysicL::raw_value(_wall_dislocations_step[_qp]));
571 164738 : _old_input_values[_strain_output_index] =
572 164738 : _creep_strain_old_forcing_function
573 164738 : ? _creep_strain_old_forcing_function->value(_t, _q_point[_qp])
574 : : MetaPhysicL::raw_value(
575 158242 : std::sqrt((creep_strain_substep).doubleContraction(creep_strain_substep) / 1.5));
576 :
577 : // Prepare input
578 164738 : _input_values[_cell_input_index] = _old_input_values[_cell_output_index];
579 164738 : _input_values[_wall_input_index] = _old_input_values[_wall_output_index];
580 322980 : _input_values[_stress_input_index] = _stress_function ? _stress_function->value(_t, _q_point[_qp])
581 158242 : : effective_trial_stress * 1.0e-6;
582 164738 : _input_values[_old_strain_input_index] = _old_input_values[_strain_output_index];
583 164738 : _input_values[_temperature_input_index] = _temperature[_qp];
584 164738 : if (_environmental)
585 0 : _input_values[_environmental_input_index] = (*_environmental)[_qp];
586 :
587 : // Determine tile weight and check if input is in range
588 336228 : for (unsigned int p = 0; p < _num_partitions; ++p)
589 171490 : std::fill(_non_stress_weights[p].begin(), _non_stress_weights[p].end(), 1.0);
590 988426 : for (unsigned int i = 0; i < _num_inputs; i++)
591 823690 : if (i != _stress_input_index)
592 658952 : computeTileWeight(_non_stress_weights, _input_values[i], i);
593 :
594 : // Precompute transformed input and prebuild polynomials for inputs other than strain
595 164736 : precomputeROM(_strain_output_index);
596 164736 : }
597 :
598 : template <bool is_ad>
599 : void
600 1122392 : LAROMANCEStressUpdateBaseTempl<is_ad>::computeTileWeight(
601 : std::vector<std::vector<GenericReal<is_ad>>> & weights,
602 : GenericReal<is_ad> & input,
603 : const unsigned int in_index,
604 : const bool derivative)
605 : {
606 : mooseAssert(std::isfinite(MetaPhysicL::raw_value(input)),
607 : "computeTileWeight input must be finite");
608 :
609 2288254 : for (unsigned int p = 0; p < _num_partitions; ++p)
610 : {
611 2549086 : for (unsigned int t = 0; t < _num_tiles[p]; ++t)
612 : {
613 : // If input is within a specfic tile's window of applicability
614 1383224 : if (input >= _input_limits[p][t][in_index][0] && input <= _input_limits[p][t][in_index][1])
615 : {
616 : // If there is not more than one tile for this input, set derivative of
617 : // weight to zero
618 1257358 : if (_tiling[p][in_index] == 1)
619 : {
620 1153374 : if (derivative)
621 223488 : weights[p][t] = 0.0;
622 : }
623 : else
624 : {
625 : // Flag to ensure weights are applied only once
626 : bool overlap = false;
627 475808 : for (unsigned int tt = 0; tt < _num_tiles[p]; ++tt)
628 : {
629 371824 : if (!overlap && t != tt)
630 : {
631 : // If input is within another tile's window of applicability, check to see if inputs
632 : // place us in that tile and ensure the two tiles are different in the dimension of
633 : // interest
634 223136 : if (areTilesNotIdentical(p, t, tt, in_index) && checkInTile(p, tt))
635 : {
636 : overlap = true;
637 :
638 : // If current tile is below the second tile's window of applicability
639 44992 : if (_input_limits[p][t][in_index][0] < _input_limits[p][tt][in_index][0] &&
640 20160 : _input_limits[p][t][in_index][1] > _input_limits[p][tt][in_index][0])
641 : {
642 20160 : weights[p][t] *= sigmoid(_input_limits[p][tt][in_index][0],
643 : _input_limits[p][t][in_index][1],
644 : input,
645 : derivative);
646 : }
647 : // If current tile is above the second tile's window of applicability
648 24832 : else if (_input_limits[p][t][in_index][0] > _input_limits[p][tt][in_index][0] &&
649 24832 : _input_limits[p][t][in_index][0] < _input_limits[p][tt][in_index][1])
650 : {
651 24832 : if (derivative)
652 9216 : weights[p][t] *= -sigmoid(_input_limits[p][t][in_index][0],
653 : _input_limits[p][tt][in_index][1],
654 : input,
655 : derivative);
656 : else
657 15616 : weights[p][t] *= (1.0 - sigmoid(_input_limits[p][t][in_index][0],
658 : _input_limits[p][tt][in_index][1],
659 : input));
660 : }
661 : }
662 : }
663 : }
664 :
665 : // If not overlapping, weight is not updated, and the derivative of tile weight is set to
666 : // zero
667 103984 : if (!overlap && derivative)
668 19832 : weights[p][t] = 0.0;
669 : }
670 : }
671 : // If input is below the lower tile limit
672 125866 : else if (input < _input_limits[p][t][in_index][0])
673 : {
674 : // If the lower tile limit equals the lower global limit
675 103497 : if (_input_limits[p][t][in_index][0] == _global_limits[in_index].first)
676 : {
677 67337 : if (_window_failure[in_index].first == WindowFailure::EXTRAPOLATE)
678 : {
679 66600 : if (derivative)
680 0 : weights[p][t] *= -sigmoid(0.0, _input_limits[p][t][in_index][0], input, derivative);
681 : else
682 120072 : weights[p][t] *= (1.0 - sigmoid(0.0, _input_limits[p][t][in_index][0], input));
683 66600 : input = _input_limits[p][t][in_index][0];
684 : }
685 737 : else if (_window_failure[in_index].first == WindowFailure::USELIMIT)
686 0 : input = _input_limits[p][t][in_index][0];
687 : else
688 : {
689 737 : weights[p][t] = 0.0;
690 737 : std::stringstream msg;
691 737 : msg << "In " << _name << ": " << _index_name[in_index]
692 2211 : << " input parameter with value (" << MetaPhysicL::raw_value(input)
693 737 : << ") is out of lower global range (" << _input_limits[p][t][in_index][0] << ")";
694 :
695 : // Set input to edge of limit so it is found later in computePartitionWeights
696 737 : input = _input_limits[p][t][in_index][0];
697 :
698 737 : if (_window_failure[in_index].first == WindowFailure::WARN)
699 736 : mooseWarning(msg.str());
700 1 : else if (_window_failure[in_index].first == WindowFailure::ERROR)
701 1 : mooseError(msg.str());
702 0 : else if (_window_failure[in_index].first == WindowFailure::EXCEPTION)
703 0 : mooseException(msg.str());
704 : // if (_window_failure[in_index].first == WindowFailure::DONOTHING) <- nothing is done
705 736 : }
706 : }
707 : // if input below tile limit, update weight of tile to be zero
708 : else
709 36160 : weights[p][t] = 0.0;
710 : }
711 : // If input is above the upper tile limit
712 22369 : else if (input > _input_limits[p][t][in_index][1])
713 : {
714 22369 : if (_input_limits[p][t][in_index][1] == _global_limits[in_index].second)
715 : {
716 1 : if (_window_failure[in_index].second == WindowFailure::EXTRAPOLATE)
717 0 : mooseError("Internal error. Extrapolate not appropriate for upper bound");
718 1 : else if (_window_failure[in_index].second == WindowFailure::USELIMIT)
719 0 : input = _input_limits[p][t][in_index][1];
720 : else
721 : {
722 1 : weights[p][t] = 0.0;
723 1 : std::stringstream msg;
724 1 : msg << "In " << _name << ": " << _index_name[in_index]
725 3 : << " input parameter with value (" << MetaPhysicL::raw_value(input)
726 1 : << ") is out of upper global range (" << _input_limits[p][t][in_index][1] << ")";
727 :
728 : // Set input to edge of limit so it is found later in computePartitionWeights
729 1 : input = _input_limits[p][t][in_index][1];
730 :
731 1 : if (_window_failure[in_index].second == WindowFailure::WARN)
732 0 : mooseWarning(msg.str());
733 1 : else if (_window_failure[in_index].second == WindowFailure::ERROR)
734 0 : mooseError(msg.str());
735 1 : else if (_window_failure[in_index].second == WindowFailure::EXCEPTION)
736 3 : mooseException(msg.str());
737 : // if (_window_failure[in_index].second == WindowFailure::DONOTHING) <- nothing is done
738 1 : }
739 : }
740 : // if input above tile limit, update weight of tile to be zero
741 : else
742 22368 : weights[p][t] = 0.0;
743 : }
744 : // If input is outside window of applicability, weight is zero
745 : else
746 0 : mooseError("In ", _name, ": Internal error. Outside input limits, input=", input);
747 : }
748 : }
749 1122390 : }
750 :
751 : template <bool is_ad>
752 : bool
753 449194 : LAROMANCEStressUpdateBaseTempl<is_ad>::checkInTile(const unsigned int p, const unsigned int t) const
754 : {
755 2340408 : for (unsigned int i = 0; i < _num_inputs; ++i)
756 2038170 : if (_input_values[i] < _input_limits[p][t][i][0] ||
757 1683406 : _input_values[i] > _input_limits[p][t][i][1])
758 : return false;
759 : return true;
760 : }
761 :
762 : template <bool is_ad>
763 : bool
764 223136 : LAROMANCEStressUpdateBaseTempl<is_ad>::areTilesNotIdentical(const unsigned int p,
765 : const unsigned int t,
766 : const unsigned int tt,
767 : const unsigned int in_index)
768 : {
769 223136 : if (_input_limits[p][t][in_index][0] != _input_limits[p][tt][in_index][0] &&
770 158656 : _input_limits[p][t][in_index][1] != _input_limits[p][tt][in_index][1])
771 : return true;
772 : else
773 64480 : return false;
774 : }
775 :
776 : template <bool is_ad>
777 : GenericReal<is_ad>
778 231720 : LAROMANCEStressUpdateBaseTempl<is_ad>::computeResidual(
779 : const GenericReal<is_ad> & effective_trial_stress, const GenericReal<is_ad> & scalar)
780 : {
781 : mooseAssert(std::isfinite(MetaPhysicL::raw_value(effective_trial_stress)),
782 : "computeResidual: effective_trial_stress must be finite");
783 : mooseAssert(std::isfinite(MetaPhysicL::raw_value(scalar)),
784 : "computeResidual: scalar must be finite!");
785 : // Update new stress
786 232344 : GenericReal<is_ad> trial_stress_mpa = _stress_function
787 231720 : ? _stress_function->value(_t, _q_point[_qp])
788 225416 : : effective_trial_stress * _stress_ucf;
789 55056 : GenericReal<is_ad> dtrial_stress_dscalar = 0.0;
790 :
791 : // Update stress if strain is being applied, i.e. non-testing simulation
792 231720 : if (this->_apply_strain)
793 : {
794 278936 : trial_stress_mpa -= this->_three_shear_modulus * scalar * _stress_ucf;
795 224808 : dtrial_stress_dscalar -= this->_three_shear_modulus * _stress_ucf;
796 : }
797 231720 : _input_values[_stress_input_index] = trial_stress_mpa;
798 :
799 : // Update weights for each partition with new stress
800 471672 : for (unsigned int p = 0; p < _num_partitions; ++p)
801 239952 : _weights[p] = _non_stress_weights[p];
802 231720 : computeTileWeight(_weights, _input_values[_stress_input_index], _stress_input_index);
803 231720 : auto dweights_dstress = _non_stress_weights;
804 231720 : computeTileWeight(
805 231720 : dweights_dstress, _input_values[_stress_input_index], _stress_input_index, true);
806 :
807 231720 : computePartitionWeights(_partition_weights, _dpartition_weight_dstress);
808 :
809 : // Save extrapolation as a material property in order quantify adequate tiling range
810 231720 : _extrapolation[_qp] = 0.0;
811 471672 : for (unsigned int p = 0; p < _num_partitions; ++p)
812 521064 : for (unsigned int t = 0; t < _num_tiles[p]; ++t)
813 336168 : _extrapolation[_qp] += MetaPhysicL::raw_value(_weights[p][t] * _partition_weights[p]);
814 :
815 55056 : GenericReal<is_ad> total_rom_effective_strain_inc = 0.0;
816 55056 : GenericReal<is_ad> dtotal_rom_effective_strain_inc_dstress = 0.0;
817 :
818 : // Run ROM if all values are within windows.
819 471672 : for (unsigned int p = 0; p < _num_partitions; p++)
820 : {
821 239952 : if (_partition_weights[p])
822 : {
823 : // compute weight normalizing factor
824 55056 : GenericReal<is_ad> weight_normalizer = 0;
825 : unsigned int number_of_active_tiles = 0;
826 497440 : for (unsigned int t = 0; t < _num_tiles[p]; ++t)
827 : {
828 : // count number of active tiles
829 262354 : number_of_active_tiles += checkInTile(p, t);
830 268114 : if (_weights[p][t])
831 : {
832 : // tile normalization factor (sum of tile weights)
833 230894 : weight_normalizer += _weights[p][t];
834 : }
835 : }
836 :
837 : // normalize weights only when 3 tiles overlap
838 235086 : if (number_of_active_tiles == 3)
839 4224 : for (unsigned int t = 0; t < _num_tiles[p]; ++t)
840 : {
841 3168 : _weights[p][t] /= weight_normalizer;
842 3168 : dweights_dstress[p][t] /= weight_normalizer;
843 : }
844 :
845 497440 : for (unsigned int t = 0; t < _num_tiles[p]; ++t)
846 268114 : if (_weights[p][t])
847 : {
848 230894 : const GenericReal<is_ad> rom = computeROM(t, p, _strain_output_index);
849 230894 : if (rom == std::numeric_limits<float>::infinity())
850 0 : mooseError("In ", _name, ": Output for strain increment reaches infinity: ", rom);
851 :
852 230894 : total_rom_effective_strain_inc += _partition_weights[p] * _weights[p][t] * rom;
853 :
854 230894 : dtotal_rom_effective_strain_inc_dstress +=
855 230894 : _partition_weights[p] * _weights[p][t] * computeROM(t, p, _strain_output_index, true);
856 230894 : if (_dpartition_weight_dstress[p])
857 6012 : dtotal_rom_effective_strain_inc_dstress +=
858 6012 : _dpartition_weight_dstress[p] * _weights[p][t] * rom;
859 230894 : if (dweights_dstress[p][t])
860 3072 : dtotal_rom_effective_strain_inc_dstress +=
861 3072 : _partition_weights[p] * dweights_dstress[p][t] * rom;
862 : }
863 : }
864 : }
865 :
866 231720 : if (_verbose)
867 : {
868 : Moose::err << std::setprecision(9);
869 0 : GenericReal<is_ad> environmental = 0.0;
870 5056 : if (_environmental)
871 0 : environmental = (*_environmental)[_qp];
872 5056 : Moose::err << "Verbose information from " << _name << ": \n";
873 5056 : Moose::err << " dt: " << _dt << "\n";
874 5056 : Moose::err << " old cell disl: " << _old_input_values[_cell_output_index] << "\n";
875 5056 : Moose::err << " old wall disl: " << _old_input_values[_wall_output_index] << "\n";
876 : Moose::err << " initial stress (MPa): "
877 5056 : << MetaPhysicL::raw_value(effective_trial_stress) * _stress_ucf << "\n";
878 5056 : Moose::err << " temperature: " << MetaPhysicL::raw_value(_temperature[_qp]) << "\n";
879 : Moose::err << " environmental factor: " << MetaPhysicL::raw_value(environmental) << "\n";
880 : Moose::err << " calculated scalar strain value: " << MetaPhysicL::raw_value(scalar) << "\n";
881 : Moose::err << " trial stress into rom (MPa): " << MetaPhysicL::raw_value(trial_stress_mpa)
882 : << "\n";
883 5056 : Moose::err << " old effective strain: " << _old_input_values[_strain_output_index] << "\n";
884 5056 : Moose::err << " extrapolation: " << MetaPhysicL::raw_value(_extrapolation[_qp]) << "\n";
885 5056 : Moose::err << " partition 2 weight: " << MetaPhysicL::raw_value(_second_partition_weight[_qp])
886 : << "\n";
887 : Moose::err << " weights by tile, partition 1: ";
888 25280 : for (unsigned int t = 0; t < _num_tiles[0]; ++t)
889 : Moose::err << " (" << t << ", " << MetaPhysicL::raw_value(_weights[0][t]) << ") ";
890 : Moose::err << "\n";
891 5056 : if (_num_partitions > 1)
892 : {
893 : Moose::err << " weights by tile, partition 2: ";
894 20224 : for (unsigned int t = 0; t < _num_tiles[1]; ++t)
895 : Moose::err << " (" << t << ", " << MetaPhysicL::raw_value(_weights[1][t]) << ") ";
896 : }
897 : Moose::err << "\n";
898 : Moose::err << " nonstress weights by tile, partition 1: ";
899 25280 : for (unsigned int t = 0; t < _num_tiles[0]; ++t)
900 : Moose::err << " (" << t << ", " << MetaPhysicL::raw_value(_non_stress_weights[0][t]) << ") ";
901 : Moose::err << "\n";
902 5056 : if (_num_partitions > 1)
903 : {
904 : Moose::err << " nonstress weights by tile, partition 2: ";
905 20224 : for (unsigned int t = 0; t < _num_tiles[1]; ++t)
906 : Moose::err << " (" << t << ", " << MetaPhysicL::raw_value(_non_stress_weights[1][t])
907 : << ") ";
908 : }
909 : Moose::err << "\n";
910 : Moose::err << " effective strain increment: "
911 : << MetaPhysicL::raw_value(total_rom_effective_strain_inc) << std::endl;
912 : }
913 :
914 286776 : _creep_rate[_qp] = total_rom_effective_strain_inc / _dt;
915 231720 : _derivative = dtotal_rom_effective_strain_inc_dstress * dtrial_stress_dscalar - 1.0;
916 :
917 231720 : if (!this->_apply_strain)
918 : {
919 6912 : if (_verbose)
920 : Moose::err << " Strain not applied due to apply_strain input parameter!" << std::endl;
921 6912 : _derivative = 1.0;
922 6912 : return 0.0;
923 : }
924 170680 : return total_rom_effective_strain_inc - scalar;
925 231720 : }
926 :
927 : template <bool is_ad>
928 : void
929 452864 : LAROMANCEStressUpdateBaseTempl<is_ad>::precomputeROM(const unsigned out_index)
930 : {
931 925856 : for (unsigned int p = 0; p < _num_partitions; ++p)
932 : {
933 1046624 : for (unsigned int t = 0; t < _num_tiles[p]; ++t)
934 : {
935 : // Only precompute for tiles that don't have zero weight
936 573632 : if (_non_stress_weights[p][t])
937 : {
938 3102528 : for (unsigned int i = 0; i < _num_inputs; ++i)
939 : {
940 2585440 : if (i != _stress_input_index)
941 : {
942 2068352 : _rom_inputs[p][t][i] =
943 2068352 : normalizeInput(_input_values[i],
944 : _transform[p][t][out_index][i],
945 : _transform_coefs[p][t][out_index][i],
946 : _transformed_normalization_limits[p][t][out_index][i]);
947 2068352 : buildPolynomials(p, _rom_inputs[p][t][i], _polynomial_inputs[p][t][i]);
948 : }
949 : }
950 517088 : precomputeValues(
951 : p, _coefs[p][t][out_index], _polynomial_inputs[p][t], _precomputed_vals[p][t]);
952 : }
953 : }
954 : }
955 452864 : }
956 :
957 : template <bool is_ad>
958 : GenericReal<is_ad>
959 745372 : LAROMANCEStressUpdateBaseTempl<is_ad>::computeROM(const unsigned int t,
960 : const unsigned int p,
961 : const unsigned out_index,
962 : const bool derivative)
963 : {
964 : // Update due to new stress
965 745372 : _rom_inputs[p][t][_stress_input_index] =
966 745372 : normalizeInput(_input_values[_stress_input_index],
967 : _transform[p][t][out_index][_stress_input_index],
968 : _transform_coefs[p][t][out_index][_stress_input_index],
969 745372 : _transformed_normalization_limits[p][t][out_index][_stress_input_index]);
970 :
971 745372 : buildPolynomials(
972 : p, _rom_inputs[p][t][_stress_input_index], _polynomial_inputs[p][t][_stress_input_index]);
973 :
974 : // Compute ROM values
975 745372 : const GenericReal<is_ad> rom_output =
976 558620 : computeValues(p, _precomputed_vals[p][t], _polynomial_inputs[p][t]);
977 :
978 : // Return converted output if not derivative
979 745372 : if (!derivative)
980 514478 : return convertOutput(p, _old_input_values, rom_output, out_index);
981 :
982 230894 : const GenericReal<is_ad> drom_input =
983 230894 : normalizeInput(_input_values[_stress_input_index],
984 : _transform[p][t][out_index][_stress_input_index],
985 : _transform_coefs[p][t][out_index][_stress_input_index],
986 230894 : _transformed_normalization_limits[p][t][out_index][_stress_input_index],
987 : derivative);
988 :
989 230894 : std::vector<GenericReal<is_ad>> dpolynomial_inputs(_degree[p], 0.0);
990 230894 : buildPolynomials(
991 230894 : p, _rom_inputs[p][t][_stress_input_index], dpolynomial_inputs, drom_input, derivative);
992 :
993 230894 : const GenericReal<is_ad> drom_output = computeValues(
994 : p, _precomputed_vals[p][t], _polynomial_inputs[p][t], dpolynomial_inputs, derivative);
995 230894 : return convertOutput(p, _old_input_values, rom_output, out_index, drom_output, derivative);
996 : }
997 :
998 : template <bool is_ad>
999 : GenericReal<is_ad>
1000 3044618 : LAROMANCEStressUpdateBaseTempl<is_ad>::normalizeInput(const GenericReal<is_ad> & input,
1001 : const ROMInputTransform transform,
1002 : const Real transform_coef,
1003 : const std::vector<Real> & transformed_limits,
1004 : const bool derivative)
1005 : {
1006 3044618 : GenericReal<is_ad> x = input;
1007 3044618 : convertValue(x, transform, transform_coef, derivative);
1008 :
1009 : // transformed_limits[2] = transformed_limits[1] - transformed_limits[0]
1010 3044618 : if (derivative)
1011 181598 : return x / transformed_limits[2];
1012 : else
1013 2813724 : return (x - transformed_limits[0]) / transformed_limits[2] - 1.0;
1014 : }
1015 :
1016 : template <bool is_ad>
1017 : void
1018 3044618 : LAROMANCEStressUpdateBaseTempl<is_ad>::buildPolynomials(
1019 : const unsigned int p,
1020 : const GenericReal<is_ad> & rom_input,
1021 : std::vector<GenericReal<is_ad>> & polynomial_inputs,
1022 : const GenericReal<is_ad> & drom_input,
1023 : const bool derivative)
1024 : {
1025 14019892 : for (unsigned int d = 0; d < _degree[p]; ++d)
1026 : {
1027 10975274 : polynomial_inputs[d] = computePolynomial(rom_input, d);
1028 :
1029 10975274 : if (derivative)
1030 1079054 : polynomial_inputs[d] = drom_input * computePolynomial(rom_input, d, derivative);
1031 : }
1032 3044618 : }
1033 :
1034 : template <bool is_ad>
1035 : void
1036 517088 : LAROMANCEStressUpdateBaseTempl<is_ad>::precomputeValues(
1037 : const unsigned int p,
1038 : const std::vector<Real> & coefs,
1039 : const std::vector<std::vector<GenericReal<is_ad>>> & polynomial_inputs,
1040 : std::vector<GenericReal<is_ad>> & precomputed)
1041 : {
1042 443723104 : for (unsigned int c = 0; c < _num_coefs[p]; ++c)
1043 : {
1044 443206016 : precomputed[c] = coefs[c];
1045 2659236096 : for (unsigned int i = 0; i < _num_inputs; ++i)
1046 2216030080 : if (i != _stress_input_index)
1047 1772824064 : precomputed[c] *= polynomial_inputs[i][_makeframe_helper[p][c + _num_coefs[p] * i]];
1048 : }
1049 517088 : }
1050 :
1051 : template <bool is_ad>
1052 : GenericReal<is_ad>
1053 976266 : LAROMANCEStressUpdateBaseTempl<is_ad>::computeValues(
1054 : const unsigned int p,
1055 : const std::vector<GenericReal<is_ad>> & precomputed,
1056 : const std::vector<std::vector<GenericReal<is_ad>>> & polynomial_inputs,
1057 : const std::vector<GenericReal<is_ad>> & dpolynomial_inputs,
1058 : const bool derivative)
1059 : {
1060 236048 : GenericReal<is_ad> rom_output = 0.0;
1061 935550516 : for (unsigned int c = 0; c < _num_coefs[p]; ++c)
1062 : {
1063 934574250 : if (!derivative)
1064 712360540 : rom_output +=
1065 712360540 : precomputed[c] *
1066 712360540 : polynomial_inputs[_stress_input_index]
1067 712360540 : [_makeframe_helper[p][c + _num_coefs[p] * _stress_input_index]];
1068 :
1069 : else
1070 222213710 : rom_output +=
1071 222213710 : precomputed[c] *
1072 222213710 : dpolynomial_inputs[_makeframe_helper[p][c + _num_coefs[p] * _stress_input_index]];
1073 : }
1074 976266 : return rom_output;
1075 : }
1076 :
1077 : template <bool is_ad>
1078 : GenericReal<is_ad>
1079 745372 : LAROMANCEStressUpdateBaseTempl<is_ad>::convertOutput(const unsigned int p,
1080 : const std::vector<Real> & old_input_values,
1081 : const GenericReal<is_ad> & rom_output,
1082 : const unsigned out_index,
1083 : const GenericReal<is_ad> & drom_output,
1084 : const bool derivative)
1085 : {
1086 745372 : if (out_index == _strain_output_index)
1087 : {
1088 461788 : if (derivative)
1089 280190 : return std::exp(rom_output) * _dt * drom_output;
1090 : else
1091 280190 : return std::exp(rom_output) * _dt;
1092 : }
1093 :
1094 283584 : if (derivative)
1095 0 : return 0.0;
1096 :
1097 283584 : GenericReal<is_ad> expout = std::exp(rom_output);
1098 : mooseAssert(expout > 0.0, "ROM calculated strain increment must be positive");
1099 :
1100 283584 : if (expout > _cutoff[p])
1101 195424 : expout -= _cutoff[p];
1102 : else
1103 0 : expout = -_cutoff[p] * _cutoff[p] / expout + _cutoff[p];
1104 :
1105 371744 : return -expout * old_input_values[out_index] * _dt;
1106 : }
1107 :
1108 : template <bool is_ad>
1109 : GenericReal<is_ad>
1110 11857144 : LAROMANCEStressUpdateBaseTempl<is_ad>::computePolynomial(const GenericReal<is_ad> & value,
1111 : const unsigned int degree,
1112 : const bool derivative)
1113 : {
1114 11857144 : if (degree == 0)
1115 : {
1116 3275512 : if (derivative)
1117 49296 : return 0.0;
1118 3044618 : return 1.0;
1119 : }
1120 8581632 : else if (degree == 1)
1121 : {
1122 2860544 : if (derivative)
1123 49296 : return 1.0;
1124 2643552 : return value;
1125 : }
1126 5721088 : else if (degree == 2)
1127 : {
1128 2860544 : if (derivative)
1129 266288 : return 3.0 * value;
1130 3547248 : return 1.5 * Utility::pow<2>(value) - 0.5;
1131 : }
1132 : else
1133 : {
1134 2860544 : if (derivative)
1135 266288 : return 7.5 * Utility::pow<2>(value) - 1.5;
1136 4450944 : return 2.5 * Utility::pow<3>(value) - 1.5 * value;
1137 : }
1138 : }
1139 :
1140 : template <bool is_ad>
1141 : std::vector<std::vector<std::vector<std::vector<Real>>>>
1142 279 : LAROMANCEStressUpdateBaseTempl<is_ad>::getTransformedLimits(
1143 : const unsigned int p, const std::vector<std::vector<std::vector<Real>>> limits)
1144 : {
1145 279 : std::vector<std::vector<std::vector<std::vector<Real>>>> transformed_limits(
1146 : _num_tiles[p],
1147 279 : std::vector<std::vector<std::vector<Real>>>(
1148 558 : _num_outputs, std::vector<std::vector<Real>>(_num_inputs, std::vector<Real>(3))));
1149 828 : for (unsigned int t = 0; t < _num_tiles[p]; ++t)
1150 : {
1151 2196 : for (unsigned int o = 0; o < _num_outputs; ++o)
1152 : {
1153 9882 : for (unsigned int i = 0; i < _num_inputs; ++i)
1154 : {
1155 24705 : for (unsigned int k = 0; k < 2; ++k)
1156 : {
1157 16470 : transformed_limits[t][o][i][k] = limits[t][i][k];
1158 16470 : convertValue(
1159 : transformed_limits[t][o][i][k], _transform[p][t][o][i], _transform_coefs[p][t][o][i]);
1160 : }
1161 8235 : transformed_limits[t][o][i][2] =
1162 8235 : (transformed_limits[t][o][i][1] - transformed_limits[t][o][i][0]) / 2.0;
1163 : }
1164 : }
1165 : }
1166 279 : return transformed_limits;
1167 : }
1168 :
1169 : template <bool is_ad>
1170 : std::vector<unsigned int>
1171 279 : LAROMANCEStressUpdateBaseTempl<is_ad>::getMakeFrameHelper(const unsigned int p) const
1172 : {
1173 279 : std::vector<unsigned int> v(_num_coefs[p] * _num_inputs);
1174 :
1175 175491 : for (unsigned int numcoeffs = 0; numcoeffs < _num_coefs[p]; ++numcoeffs)
1176 1051272 : for (unsigned int invar = 0; invar < _num_inputs; ++invar)
1177 876060 : v[numcoeffs + _num_coefs[p] * invar] =
1178 876060 : numcoeffs / MathUtils::pow(_degree[p], invar) % _degree[p];
1179 :
1180 279 : return v;
1181 : }
1182 :
1183 : template <bool is_ad>
1184 : void
1185 164736 : LAROMANCEStressUpdateBaseTempl<is_ad>::computeStressFinalize(
1186 : const GenericRankTwoTensor<is_ad> & plastic_strain_increment)
1187 : {
1188 164736 : _cell_dislocation_increment = 0.0;
1189 164736 : _wall_dislocation_increment = 0.0;
1190 182032 : if (_input_values[_stress_input_index])
1191 : {
1192 144064 : precomputeROM(_cell_output_index);
1193 294816 : for (unsigned int p = 0; p < _num_partitions; ++p)
1194 150784 : if (_partition_weights[p])
1195 312128 : for (unsigned int t = 0; t < _num_tiles[p]; ++t)
1196 171904 : if (_weights[p][t])
1197 141792 : _cell_dislocation_increment +=
1198 141792 : _partition_weights[p] * _weights[p][t] * computeROM(t, p, _cell_output_index);
1199 144064 : precomputeROM(_wall_output_index);
1200 294816 : for (unsigned int p = 0; p < _num_partitions; ++p)
1201 150784 : if (_partition_weights[p])
1202 312128 : for (unsigned int t = 0; t < _num_tiles[p]; ++t)
1203 171904 : if (_weights[p][t])
1204 141792 : _wall_dislocation_increment +=
1205 141792 : _partition_weights[p] * _weights[p][t] * computeROM(t, p, _wall_output_index);
1206 : }
1207 :
1208 231904 : _cell_rate[_qp] = _cell_dislocation_increment / _dt;
1209 164736 : _wall_rate[_qp] = _wall_dislocation_increment / _dt;
1210 231904 : _cell_dislocations[_qp] = _old_input_values[_cell_output_index] + _cell_dislocation_increment;
1211 231904 : _wall_dislocations[_qp] = _old_input_values[_wall_output_index] + _wall_dislocation_increment;
1212 :
1213 : // For (possibly) substepping.
1214 164736 : _plastic_strain_increment += MetaPhysicL::raw_value(plastic_strain_increment);
1215 :
1216 : // Prevent the ROM from calculating and proceeding with negative dislocations
1217 164736 : if (_apply_strain && (_cell_dislocations[_qp] < 0.0 || _wall_dislocations[_qp] < 0.0))
1218 0 : mooseException("In ",
1219 : _name,
1220 : ": Negative disclocation density calculated for cell (old : ",
1221 : MetaPhysicL::raw_value(_old_input_values[_cell_output_index]),
1222 : " increment: ",
1223 : MetaPhysicL::raw_value(_cell_dislocation_increment),
1224 : " value: ",
1225 : MetaPhysicL::raw_value(_cell_dislocations[_qp]),
1226 : ") or wall (old : ",
1227 : MetaPhysicL::raw_value(_old_input_values[_wall_output_index]),
1228 : " increment: ",
1229 : MetaPhysicL::raw_value(_wall_dislocation_increment),
1230 : " value: ",
1231 : MetaPhysicL::raw_value(_wall_dislocations[_qp]),
1232 : ").");
1233 :
1234 164736 : if (_verbose)
1235 : {
1236 : Moose::err << " Finalized ROM output\n";
1237 : Moose::err << " effective creep strain increment: "
1238 5184 : << std::sqrt(2.0 / 3.0 *
1239 5184 : MetaPhysicL::raw_value(_plastic_strain_increment.doubleContraction(
1240 : _plastic_strain_increment)))
1241 : << "\n";
1242 : Moose::err << " total effective creep strain: "
1243 5184 : << std::sqrt(2.0 / 3.0 *
1244 5184 : MetaPhysicL::raw_value(this->_creep_strain[_qp].doubleContraction(
1245 : this->_creep_strain[_qp])))
1246 : << "\n";
1247 5184 : Moose::err << " creep rate: " << MetaPhysicL::raw_value(_creep_rate[_qp]) << "\n";
1248 5184 : Moose::err << " cell dislocation rate: " << MetaPhysicL::raw_value(_cell_rate[_qp]) << "\n";
1249 5184 : Moose::err << " wall dislocation rate: " << MetaPhysicL::raw_value(_wall_rate[_qp]) << "\n";
1250 5184 : Moose::err << " new cell dislocations: " << MetaPhysicL::raw_value(_cell_dislocations[_qp])
1251 : << "\n";
1252 5184 : Moose::err << " new wall dislocations: " << MetaPhysicL::raw_value(_wall_dislocations[_qp])
1253 : << "\n"
1254 : << std::endl;
1255 : }
1256 :
1257 231904 : RadialReturnCreepStressUpdateBaseTempl<is_ad>::computeStressFinalize(
1258 : MetaPhysicL::raw_value(_plastic_strain_increment));
1259 164736 : }
1260 :
1261 : template <bool is_ad>
1262 : Real
1263 146712 : LAROMANCEStressUpdateBaseTempl<is_ad>::computeTimeStepLimit()
1264 : {
1265 146712 : Real limited_dt = RadialReturnCreepStressUpdateBaseTempl<is_ad>::computeTimeStepLimit();
1266 :
1267 : Real cell_strain_inc = std::abs(MetaPhysicL::raw_value(_cell_dislocation_increment));
1268 146712 : if (cell_strain_inc && _old_input_values[_cell_output_index])
1269 120088 : limited_dt = std::min(limited_dt,
1270 228136 : _dt * _max_cell_increment * _old_input_values[_cell_output_index] /
1271 : cell_strain_inc);
1272 : Real wall_strain_inc = std::abs(MetaPhysicL::raw_value(_wall_dislocation_increment));
1273 146712 : if (wall_strain_inc && _old_input_values[_wall_output_index])
1274 120088 : limited_dt = std::min(limited_dt,
1275 120312 : _dt * _max_wall_increment * _old_input_values[_wall_output_index] /
1276 : wall_strain_inc);
1277 :
1278 146712 : return limited_dt;
1279 : }
1280 :
1281 : template <bool is_ad>
1282 : GenericReal<is_ad>
1283 144520 : LAROMANCEStressUpdateBaseTempl<is_ad>::sigmoid(const Real lower,
1284 : const Real upper,
1285 : const GenericReal<is_ad> & val,
1286 : const bool derivative)
1287 : {
1288 : mooseAssert(std::isfinite(MetaPhysicL::raw_value(val)), "Sigmoid value should must be infinite");
1289 : mooseAssert(MetaPhysicL::raw_value(val) >= lower, "Input value must be greater than lower limit");
1290 : mooseAssert(MetaPhysicL::raw_value(val) <= upper, "Input value must be greater than upper limit");
1291 :
1292 : // normalize val between 0 and 1, then shift between -1 and 1
1293 144520 : const GenericReal<is_ad> x = 2.0 * (val - lower) / (upper - lower) - 1.0;
1294 :
1295 144520 : if (x == -1.0)
1296 : {
1297 19012 : if (derivative)
1298 0 : return 0.0;
1299 : else
1300 9282 : return 1.0;
1301 : }
1302 125508 : else if (x == 1.0)
1303 : {
1304 0 : if (derivative)
1305 0 : return 0.0;
1306 : else
1307 0 : return 0.0;
1308 : }
1309 108608 : else if (x < 1.0 && x > -1.0)
1310 : {
1311 215552 : const GenericReal<is_ad> plus = std::exp(-2.0 / (1.0 + x));
1312 215552 : const GenericReal<is_ad> minus = std::exp(-2.0 / (1.0 - x));
1313 :
1314 108608 : if (!derivative)
1315 146900 : return 1.0 - plus / (plus + minus);
1316 :
1317 15180 : const GenericReal<is_ad> dplus_dx = plus * 2.0 / Utility::pow<2>(1.0 + x);
1318 15180 : const GenericReal<is_ad> dminus_dx = -minus * 2.0 / Utility::pow<2>(1.0 - x);
1319 :
1320 15180 : return (plus * dminus_dx - dplus_dx * minus) / Utility::pow<2>(plus + minus) * 2.0 /
1321 15180 : (upper - lower);
1322 : }
1323 : else
1324 0 : mooseError("Internal error: Sigmoid, value: x is out of bounds. val=",
1325 : val,
1326 : " low=",
1327 : lower,
1328 : " high=",
1329 : upper);
1330 : }
1331 :
1332 : template <bool is_ad>
1333 : void
1334 0 : LAROMANCEStressUpdateBaseTempl<is_ad>::outputIterationSummary(std::stringstream * iter_output,
1335 : const unsigned int total_it)
1336 : {
1337 0 : if (iter_output)
1338 : {
1339 0 : *iter_output << "At element " << this->_current_elem->id() << " _qp=" << _qp << " Coordinates "
1340 0 : << _q_point[_qp] << " block=" << this->_current_elem->subdomain_id() << '\n';
1341 0 : *iter_output << " dt " << _dt << " old cell disl: " << _old_input_values[_cell_output_index]
1342 0 : << " old wall disl: " << _old_input_values[_wall_output_index]
1343 0 : << " old effective strain: " << _old_input_values[_strain_output_index] << "\n";
1344 :
1345 0 : *iter_output << " temp: " << MetaPhysicL::raw_value(_temperature[_qp]) << " environmental: "
1346 0 : << (_environmental ? MetaPhysicL::raw_value((*_environmental)[_qp]) : 0.0)
1347 0 : << " trial stress into rom (MPa): "
1348 0 : << MetaPhysicL::raw_value(_input_values[_stress_input_index])
1349 0 : << " cell: " << MetaPhysicL::raw_value(_input_values[_cell_input_index])
1350 0 : << " wall: " << MetaPhysicL::raw_value(_input_values[_wall_input_index])
1351 0 : << " old strain: "
1352 0 : << MetaPhysicL::raw_value(_input_values[_old_strain_input_index]) << "\n";
1353 0 : *iter_output << " partition 2 weight: "
1354 0 : << MetaPhysicL::raw_value(_second_partition_weight[_qp]) << "\n";
1355 0 : *iter_output << " weights by tile, partition 1: ";
1356 0 : for (unsigned int t = 0; t < _num_tiles[0]; ++t)
1357 0 : *iter_output << " (" << t << ", " << MetaPhysicL::raw_value(_weights[0][t]) << ") ";
1358 0 : *iter_output << "\n";
1359 0 : *iter_output << " weights by tile, partition 2: ";
1360 0 : for (unsigned int t = 0; t < _num_tiles[1]; ++t)
1361 0 : *iter_output << " (" << t << ", " << MetaPhysicL::raw_value(_weights[1][t]) << ") ";
1362 0 : *iter_output << "\n";
1363 0 : *iter_output << " nonstress weights by tile, partition 1: ";
1364 0 : for (unsigned int t = 0; t < _num_tiles[0]; ++t)
1365 0 : *iter_output << " (" << t << ", " << MetaPhysicL::raw_value(_non_stress_weights[0][t])
1366 0 : << ") ";
1367 0 : *iter_output << "\n";
1368 0 : *iter_output << " nonstress weights by tile, partition 2: ";
1369 0 : for (unsigned int t = 0; t < _num_tiles[1]; ++t)
1370 0 : *iter_output << " (" << t << ", " << MetaPhysicL::raw_value(_non_stress_weights[1][t])
1371 0 : << ") ";
1372 0 : *iter_output << "\n";
1373 : }
1374 :
1375 0 : SingleVariableReturnMappingSolutionTempl<is_ad>::outputIterationSummary(iter_output, total_it);
1376 0 : }
1377 :
1378 : template <bool is_ad>
1379 : void
1380 153126 : LAROMANCEStressUpdateBaseTempl<is_ad>::outputIterationStep(
1381 : std::stringstream * iter_output,
1382 : const GenericReal<is_ad> & effective_trial_stress,
1383 : const GenericReal<is_ad> & scalar,
1384 : const Real reference_residual)
1385 : {
1386 153126 : SingleVariableReturnMappingSolutionTempl<is_ad>::outputIterationStep(
1387 : iter_output, effective_trial_stress, scalar, reference_residual);
1388 153126 : if (iter_output)
1389 0 : *iter_output << " derivative: "
1390 0 : << MetaPhysicL::raw_value(computeDerivative(effective_trial_stress, scalar))
1391 : << std::endl;
1392 153126 : }
1393 :
1394 : template <bool is_ad>
1395 : void
1396 261 : LAROMANCEStressUpdateBaseTempl<is_ad>::checkJSONKey(const std::string & key)
1397 : {
1398 522 : if (!this->isParamValid("model"))
1399 0 : this->paramError("model", "Specify a JSON data filename.");
1400 :
1401 522 : const auto model_file_name = this->_pars.rawParamVal("model");
1402 261 : if (_json.empty())
1403 0 : this->paramError("model", "The specified JSON data file '", model_file_name, "' is empty.");
1404 261 : if (!_json.contains(key))
1405 0 : this->paramError(
1406 : "model", "The key '", key, "' is missing from the JSON data file '", model_file_name, "'.");
1407 261 : }
1408 :
1409 : template class LAROMANCEStressUpdateBaseTempl<false>;
1410 : template class LAROMANCEStressUpdateBaseTempl<true>;
|
