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