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 688 : LAROMANCEStressUpdateBaseTempl<is_ad>::validParams()
27 : {
28 688 : InputParameters params = RadialReturnCreepStressUpdateBaseTempl<is_ad>::validParams();
29 688 : 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 1376 : params.addRequiredCoupledVar("temperature", "The coupled temperature (K)");
34 1376 : params.addParam<MaterialPropertyName>("environmental_factor",
35 : "Optional coupled environmental factor");
36 :
37 1376 : MooseEnum error_lower_limit_behavior("ERROR EXCEPTION WARN IGNORE DONTHING USELIMIT",
38 : "EXCEPTION");
39 : // Only allow ERROR and EXCEPTION on upper bounds
40 1376 : MooseEnum error_upper_limit_behavior("ERROR EXCEPTION", "EXCEPTION");
41 1376 : 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 1376 : 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 1376 : 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 1376 : 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 1376 : 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 1376 : 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 1376 : MooseEnum extrapolated_lower_limit_behavior(
69 : "ERROR EXCEPTION WARN IGNORE DONOTHING USELIMIT EXTRAPOLATE", "EXTRAPOLATE");
70 1376 : 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 1376 : 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 1376 : 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 1376 : 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 1376 : 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 1376 : 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 1376 : 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 2064 : params.addRangeCheckedParam<Real>(
100 : "cell_dislocations_normal_distribution_width",
101 1376 : 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 2064 : params.addRangeCheckedParam<Real>(
106 : "max_relative_cell_dislocation_increment",
107 1376 : 0.5,
108 : "max_relative_cell_dislocation_increment > 0.0",
109 : "Maximum increment of density of cell (glissile) dislocations.");
110 :
111 1376 : 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 2064 : params.addRangeCheckedParam<Real>(
116 : "wall_dislocations_normal_distribution_width",
117 1376 : 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 2064 : params.addRangeCheckedParam<Real>(
122 : "max_relative_wall_dislocation_increment",
123 1376 : 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 1376 : params.addParam<bool>("verbose", false, "Flag to output verbose information.");
128 :
129 1376 : 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 1376 : 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 1376 : 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 1376 : 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 1376 : params.addParam<unsigned int>("seed", 0, "Random number generator seed");
151 1376 : params.addParam<std::string>("stress_unit", "Pa", "unit of stress");
152 :
153 : // use std::string here to avoid automatic absolute path expansion
154 1376 : params.addParam<DataFileName>("model", "LaRomance model JSON datafile");
155 1376 : params.addParam<FileName>("export_model", "Write LaRomance model to JSON datafile");
156 :
157 1376 : 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 688 : return params;
163 688 : }
164 :
165 : template <bool is_ad>
166 516 : LAROMANCEStressUpdateBaseTempl<is_ad>::LAROMANCEStressUpdateBaseTempl(
167 : const InputParameters & parameters)
168 : : RadialReturnCreepStressUpdateBaseTempl<is_ad>(parameters),
169 516 : _temperature(this->template coupledGenericValue<is_ad>("temperature")),
170 516 : _environmental(
171 516 : this->isParamValid("environmental_factor")
172 516 : ? &this->template getGenericMaterialProperty<Real, is_ad>("environmental_factor")
173 : : nullptr),
174 1032 : _window_failure(_environmental ? 6 : 5),
175 :
176 1032 : _verbose(this->template getParam<bool>("verbose")),
177 516 : _cell_dislocations(
178 516 : this->template declareGenericProperty<Real, is_ad>(this->_base_name + "cell_dislocations")),
179 516 : _cell_dislocations_old(
180 516 : this->template getMaterialPropertyOld<Real>(this->_base_name + "cell_dislocations")),
181 1032 : _max_cell_increment(this->template getParam<Real>("max_relative_cell_dislocation_increment")),
182 516 : _cell_function(this->isParamValid("cell_dislocation_density_forcing_function")
183 630 : ? &this->getFunction("cell_dislocation_density_forcing_function")
184 : : NULL),
185 516 : _cell_dislocation_increment(0.0),
186 516 : _wall_dislocations(
187 516 : this->template declareGenericProperty<Real, is_ad>(this->_base_name + "wall_dislocations")),
188 516 : _wall_dislocations_old(
189 516 : this->template getMaterialPropertyOld<Real>(this->_base_name + "wall_dislocations")),
190 1032 : _max_wall_increment(this->template getParam<Real>("max_relative_wall_dislocation_increment")),
191 516 : _wall_function(this->isParamValid("wall_dislocation_density_forcing_function")
192 630 : ? &this->getFunction("wall_dislocation_density_forcing_function")
193 : : NULL),
194 516 : _stress_function(this->isParamValid("effective_stress_forcing_function")
195 630 : ? &this->getFunction("effective_stress_forcing_function")
196 : : NULL),
197 516 : _wall_dislocation_increment(0.0),
198 516 : _cell_input_index(0),
199 516 : _wall_input_index(1),
200 516 : _stress_input_index(2),
201 516 : _old_strain_input_index(3),
202 516 : _temperature_input_index(4),
203 516 : _environmental_input_index(5),
204 516 : _cell_output_index(0),
205 516 : _wall_output_index(1),
206 516 : _strain_output_index(2),
207 :
208 1032 : _creep_strain_old_forcing_function(this->isParamValid("old_creep_strain_forcing_function")
209 630 : ? &this->getFunction("old_creep_strain_forcing_function")
210 : : NULL),
211 :
212 516 : _creep_rate(
213 516 : this->template declareGenericProperty<Real, is_ad>(this->_base_name + "creep_rate")),
214 1032 : _cell_rate(this->template declareGenericProperty<Real, is_ad>(this->_base_name +
215 : "cell_dislocation_rate")),
216 1032 : _wall_rate(this->template declareGenericProperty<Real, is_ad>(this->_base_name +
217 : "wall_dislocation_rate")),
218 516 : _extrapolation(this->template declareProperty<Real>("ROM_extrapolation")),
219 516 : _second_partition_weight(
220 516 : this->template declareGenericProperty<Real, is_ad>("partition_weight")),
221 :
222 516 : _derivative(0.0),
223 516 : _old_input_values(3),
224 1032 : _wall_dislocations_step(this->template declareGenericProperty<Real, is_ad>(
225 : this->_base_name + "wall_dislocations_step")),
226 1032 : _cell_dislocations_step(this->template declareGenericProperty<Real, is_ad>(
227 : this->_base_name + "cell_dislocations_step")),
228 516 : _plastic_strain_increment(),
229 516 : _number_of_substeps(
230 516 : this->template declareProperty<Real>(this->_base_name + "number_of_substeps")),
231 1032 : _index_name(_window_failure.size())
232 : {
233 516 : this->_check_range = true; // this may not be necessary?
234 :
235 516 : _index_name[_cell_input_index] = "cell";
236 516 : _index_name[_wall_input_index] = "wall";
237 516 : _index_name[_stress_input_index] = "stress";
238 516 : _index_name[_old_strain_input_index] = "old strain";
239 516 : _index_name[_temperature_input_index] = "temperature";
240 :
241 516 : _window_failure[_cell_input_index].first =
242 516 : parameters.get<MooseEnum>("cell_input_window_low_failure").getEnum<WindowFailure>();
243 516 : _window_failure[_cell_input_index].second =
244 516 : parameters.get<MooseEnum>("cell_input_window_high_failure").getEnum<WindowFailure>();
245 516 : _window_failure[_wall_input_index].first =
246 516 : parameters.get<MooseEnum>("wall_input_window_low_failure").getEnum<WindowFailure>();
247 516 : _window_failure[_wall_input_index].second =
248 516 : parameters.get<MooseEnum>("wall_input_window_high_failure").getEnum<WindowFailure>();
249 516 : _window_failure[_stress_input_index].first =
250 516 : parameters.get<MooseEnum>("stress_input_window_low_failure").getEnum<WindowFailure>();
251 516 : _window_failure[_stress_input_index].second =
252 516 : parameters.get<MooseEnum>("stress_input_window_high_failure").getEnum<WindowFailure>();
253 516 : _window_failure[_old_strain_input_index].first =
254 516 : parameters.get<MooseEnum>("old_strain_input_window_low_failure").getEnum<WindowFailure>();
255 516 : _window_failure[_old_strain_input_index].second =
256 516 : parameters.get<MooseEnum>("old_strain_input_window_high_failure").getEnum<WindowFailure>();
257 516 : _window_failure[_temperature_input_index].first =
258 516 : parameters.get<MooseEnum>("temperature_input_window_low_failure").getEnum<WindowFailure>();
259 516 : _window_failure[_temperature_input_index].second =
260 516 : parameters.get<MooseEnum>("temperature_input_window_high_failure").getEnum<WindowFailure>();
261 516 : 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 1032 : if (this->isParamValid("model"))
272 : {
273 126 : const auto model_file_name = this->template getParam<DataFileName>("model");
274 63 : std::ifstream model_file(model_file_name.c_str());
275 63 : model_file >> _json;
276 63 : }
277 :
278 516 : setupUnitConversionFactors(parameters);
279 516 : }
280 :
281 : template <bool is_ad>
282 : void
283 21 : LAROMANCEStressUpdateBaseTempl<is_ad>::exportJSON()
284 : {
285 42 : _json["strain_cutoff"] = getStrainCutoff();
286 42 : _json["transform"] = getTransform();
287 42 : _json["transform_coefs"] = getTransformCoefs();
288 42 : _json["input_limits"] = getInputLimits();
289 42 : _json["normalization_limits"] = getNormalizationLimits();
290 42 : _json["coefs"] = getCoefs();
291 42 : _json["tiling"] = getTilings();
292 42 : _json["cutoff"] = getStrainCutoff();
293 21 : }
294 :
295 : template <bool is_ad>
296 : bool
297 193399 : LAROMANCEStressUpdateBaseTempl<is_ad>::substeppingCapabilityEnabled()
298 : {
299 193399 : return this->_use_substepping != RadialReturnStressUpdateTempl<is_ad>::SubsteppingType::NONE;
300 : }
301 :
302 : template <bool is_ad>
303 : void
304 516 : LAROMANCEStressUpdateBaseTempl<is_ad>::setupUnitConversionFactors(
305 : const InputParameters & parameters)
306 : {
307 : // Stress unit conversion factor
308 1032 : const MooseUnits stress_unit_to("MPa");
309 516 : const MooseUnits stress_unit_from(parameters.get<std::string>("stress_unit"));
310 516 : _stress_ucf = stress_unit_to.convert(1, stress_unit_from);
311 516 : }
312 :
313 : template <bool is_ad>
314 : void
315 516 : LAROMANCEStressUpdateBaseTempl<is_ad>::initialSetup()
316 : {
317 : // export models that are compiled in
318 1032 : if (this->isParamValid("export_model"))
319 : {
320 21 : exportJSON();
321 42 : std::ofstream out(this->template getParam<FileName>("export_model").c_str());
322 21 : out << _json;
323 21 : }
324 :
325 : // Pull in relevant ROM information and do sanity checks
326 516 : _transform = getTransform();
327 516 : _transform_coefs = getTransformCoefs();
328 516 : _input_limits = getInputLimits();
329 516 : _normalization_limits = getNormalizationLimits();
330 516 : _coefs = getCoefs();
331 516 : _tiling = getTilings();
332 516 : _cutoff = getStrainCutoff();
333 : // resize containers to be filled later based on partition dimension
334 : // and immediately run some sanity checks
335 516 : _num_partitions = _transform.size();
336 516 : 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 516 : if (_transform[0].size() < 1)
340 0 : mooseError("In ", _name, ": Transform is not the correct shape");
341 516 : _num_outputs = _transform[0][0].size();
342 516 : 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 516 : _num_inputs = _transform[0][0][0].size();
350 516 : 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 516 : 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 516 : 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 516 : _num_tiles.resize(_num_partitions);
365 516 : _num_coefs.resize(_num_partitions);
366 516 : _degree.resize(_num_partitions);
367 516 : _precomputed_vals.resize(_num_partitions);
368 516 : _rom_inputs.resize(_num_partitions);
369 516 : _polynomial_inputs.resize(_num_partitions);
370 516 : _non_stress_weights.resize(_num_partitions);
371 516 : _weights.resize(_num_partitions);
372 516 : _partition_weights.resize(_num_partitions);
373 516 : _dpartition_weight_dstress.resize(_num_partitions);
374 516 : _transformed_normalization_limits.resize(_num_partitions);
375 516 : _makeframe_helper.resize(_num_partitions);
376 516 : _global_limits.resize(_num_inputs);
377 : // temporarily fill global limits with extreme numerical values, to later update
378 3096 : for (unsigned int i = 0; i < _num_inputs; ++i)
379 2580 : _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 516 : if (_transform.size() != _num_partitions || _transform_coefs.size() != _num_partitions ||
384 516 : _input_limits.size() != _num_partitions || _normalization_limits.size() != _num_partitions ||
385 1032 : _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 1158 : for (unsigned int p = 0; p < _num_partitions; ++p)
391 : {
392 642 : _num_tiles[p] = _transform[p].size();
393 642 : if (!_num_tiles[p])
394 0 : mooseError("In ", _name, ": No tiles detected. Double check your ROM input");
395 :
396 : bool correct_shape = true;
397 642 : if (_transform[p].size() != _num_tiles[p] || _transform_coefs[p].size() != _num_tiles[p] ||
398 642 : _input_limits[p].size() != _num_tiles[p] ||
399 1284 : _normalization_limits[p].size() != _num_tiles[p] || _coefs[p].size() != _num_tiles[p])
400 : correct_shape = false;
401 642 : if (_tiling[p].size() != _num_inputs)
402 : correct_shape = false;
403 642 : if (_coefs[p][0].size() == 0)
404 : correct_shape = false;
405 642 : _num_coefs[p] = _coefs[p][0][0].size();
406 :
407 : // start loop over tiles to perform sanity checks.
408 1914 : for (unsigned int t = 0; t < _num_tiles[p]; ++t)
409 : {
410 1272 : if (_transform[p][t].size() != _num_outputs ||
411 2544 : _transform_coefs[p][t].size() != _num_outputs || _coefs[p][t].size() != _num_outputs)
412 : correct_shape = false;
413 5088 : for (unsigned int o = 0; o < _num_outputs; ++o)
414 3816 : if (_transform[p][t][o].size() != _num_inputs ||
415 7632 : _transform_coefs[p][t][o].size() != _num_inputs ||
416 : _coefs[p][t][o].size() != _num_coefs[p])
417 : correct_shape = false;
418 1272 : if (_input_limits[p][t].size() != _num_inputs ||
419 : _normalization_limits[p][t].size() != _num_inputs)
420 : correct_shape = false;
421 7632 : for (unsigned int i = 0; i < _num_inputs; ++i)
422 6360 : if (_input_limits[p][t][i].size() != 2 || _normalization_limits[p][t][i].size() != 2)
423 : correct_shape = false;
424 : }
425 :
426 642 : if (!correct_shape)
427 0 : mooseError("In ", _name, ": ROM data is not the right shape.");
428 :
429 642 : _degree[p] = std::pow(_num_coefs[p], 1.0 / _num_inputs);
430 642 : 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 3852 : for (unsigned int i = 0; i < _num_inputs; ++i)
438 : {
439 9570 : for (unsigned int t = 0; t < _num_tiles[p]; ++t)
440 : {
441 6360 : if (_input_limits[p][t][i][0] >= _input_limits[p][t][i][1])
442 0 : mooseError("In ", _name, ": Input limits are ordered incorrectly");
443 9192 : _global_limits[i].first = std::min(_global_limits[i].first, _input_limits[p][t][i][0]);
444 9192 : _global_limits[i].second = std::max(_global_limits[i].second, _input_limits[p][t][i][1]);
445 : }
446 : }
447 :
448 : // Precompute helper containers
449 1284 : _transformed_normalization_limits[p] = getTransformedLimits(p, _normalization_limits[p]);
450 1284 : _makeframe_helper[p] = getMakeFrameHelper(p);
451 :
452 : // Prepare containers for each partition
453 642 : _precomputed_vals[p].resize(_num_tiles[p], std::vector<GenericReal<is_ad>>(_num_coefs[p]));
454 642 : _rom_inputs[p].resize(_num_tiles[p], std::vector<GenericReal<is_ad>>(_num_inputs));
455 642 : _polynomial_inputs[p].resize(_num_tiles[p],
456 642 : std::vector<std::vector<GenericReal<is_ad>>>(
457 1284 : _num_inputs, std::vector<GenericReal<is_ad>>(_degree[p])));
458 642 : _non_stress_weights[p].resize(_num_tiles[p]);
459 642 : _weights[p].resize(_num_tiles[p], 0);
460 : }
461 : // Prepare containers independent of partition
462 516 : _input_values.resize(_num_inputs);
463 :
464 516 : if (_verbose)
465 : {
466 84 : Moose::err << "ROM model info: " << _name << "\n";
467 : Moose::err << " number of tiles, partition 1: " << _num_tiles[0] << "\n";
468 84 : 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 84 : 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 84 : 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 84 : << _global_limits[_cell_input_index].first << " - "
480 84 : << _global_limits[_cell_input_index].second << ")\n";
481 84 : Moose::err << " wall dislocations (" << _global_limits[_wall_input_index].first << " - "
482 84 : << _global_limits[_wall_input_index].second << ")\n";
483 84 : Moose::err << " Stress (" << _global_limits[_stress_input_index].first << " - "
484 84 : << _global_limits[_stress_input_index].second << ")\n";
485 84 : Moose::err << " Old strain (" << _global_limits[_old_strain_input_index].first << " - "
486 84 : << _global_limits[_old_strain_input_index].second << ")\n";
487 84 : Moose::err << " Temperature (" << _global_limits[_temperature_input_index].first << " - "
488 84 : << _global_limits[_temperature_input_index].second << ")\n";
489 84 : 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 516 : }
495 :
496 : template <bool is_ad>
497 : void
498 4976 : LAROMANCEStressUpdateBaseTempl<is_ad>::initQpStatefulProperties()
499 : {
500 : MooseRandom rng;
501 9952 : rng.seed(0, this->template getParam<unsigned int>("seed"));
502 :
503 4976 : _cell_dislocations[_qp] = rng.randNormal(
504 9952 : this->template getParam<Real>("initial_cell_dislocation_density"),
505 4976 : this->template getParam<Real>("initial_cell_dislocation_density") *
506 9952 : this->template getParam<Real>("cell_dislocations_normal_distribution_width"));
507 4976 : _wall_dislocations[_qp] = rng.randNormal(
508 9952 : this->template getParam<Real>("initial_wall_dislocation_density"),
509 4976 : this->template getParam<Real>("initial_wall_dislocation_density") *
510 9952 : this->template getParam<Real>("wall_dislocations_normal_distribution_width"));
511 :
512 4976 : RadialReturnCreepStressUpdateBaseTempl<is_ad>::initQpStatefulProperties();
513 4976 : }
514 :
515 : template <bool is_ad>
516 : GenericReal<is_ad>
517 191088 : 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 191088 : return effective_trial_stress / this->_three_shear_modulus * 0.999999;
523 : }
524 :
525 : template <bool is_ad>
526 : void
527 193360 : LAROMANCEStressUpdateBaseTempl<is_ad>::resetIncrementalMaterialProperties()
528 : {
529 193360 : _cell_dislocation_increment = 0.0;
530 193360 : _wall_dislocation_increment = 0.0;
531 :
532 : _plastic_strain_increment.zero();
533 :
534 193360 : _wall_dislocations_step[_qp] = 0.0;
535 193360 : _cell_dislocations_step[_qp] = 0.0;
536 193360 : }
537 :
538 : template <bool is_ad>
539 : void
540 37344 : LAROMANCEStressUpdateBaseTempl<is_ad>::storeIncrementalMaterialProperties(
541 : const unsigned int total_number_of_substeps)
542 : {
543 37344 : _wall_dislocations_step[_qp] += _wall_dislocation_increment;
544 37344 : _cell_dislocations_step[_qp] += _cell_dislocation_increment;
545 37344 : _number_of_substeps[_qp] = total_number_of_substeps;
546 37344 : }
547 :
548 : template <bool is_ad>
549 : void
550 216140 : LAROMANCEStressUpdateBaseTempl<is_ad>::computeStressInitialize(
551 : const GenericReal<is_ad> & effective_trial_stress,
552 : const GenericRankFourTensor<is_ad> & elasticity_tensor)
553 : {
554 84216 : RadialReturnCreepStressUpdateBaseTempl<is_ad>::computeStressInitialize(effective_trial_stress,
555 : elasticity_tensor);
556 : // Previous substep creep strain
557 216140 : RankTwoTensor creep_strain_substep = this->_creep_strain_old[_qp] + _plastic_strain_increment;
558 :
559 : // Prepare old values
560 216140 : _old_input_values[_cell_output_index] =
561 216140 : _cell_function
562 216140 : ? _cell_function->value(_t, _q_point[_qp])
563 203516 : : (_cell_dislocations_old[_qp] + MetaPhysicL::raw_value(_cell_dislocations_step[_qp]));
564 216140 : _old_input_values[_wall_output_index] =
565 216140 : _wall_function
566 216140 : ? _wall_function->value(_t, _q_point[_qp])
567 203516 : : (_wall_dislocations_old[_qp] + MetaPhysicL::raw_value(_wall_dislocations_step[_qp]));
568 216140 : _old_input_values[_strain_output_index] =
569 216140 : _creep_strain_old_forcing_function
570 216140 : ? _creep_strain_old_forcing_function->value(_t, _q_point[_qp])
571 : : MetaPhysicL::raw_value(
572 203516 : std::sqrt((creep_strain_substep).doubleContraction(creep_strain_substep) / 1.5));
573 :
574 : // Prepare input
575 216140 : _input_values[_cell_input_index] = _old_input_values[_cell_output_index];
576 216140 : _input_values[_wall_input_index] = _old_input_values[_wall_output_index];
577 419656 : _input_values[_stress_input_index] = _stress_function ? _stress_function->value(_t, _q_point[_qp])
578 203516 : : effective_trial_stress * _stress_ucf;
579 216140 : _input_values[_old_strain_input_index] = _old_input_values[_strain_output_index];
580 216140 : _input_values[_temperature_input_index] = _temperature[_qp];
581 216140 : if (_environmental)
582 0 : _input_values[_environmental_input_index] = (*_environmental)[_qp];
583 :
584 : // Determine tile weight and check if input is in range
585 445992 : for (unsigned int p = 0; p < _num_partitions; ++p)
586 229852 : std::fill(_non_stress_weights[p].begin(), _non_stress_weights[p].end(), 1.0);
587 1296836 : for (unsigned int i = 0; i < _num_inputs; i++)
588 1080700 : if (i != _stress_input_index)
589 864560 : computeTileWeight(_non_stress_weights, _input_values[i], i);
590 :
591 : // Precompute transformed input and prebuild polynomials for inputs other than strain
592 216136 : precomputeROM(_strain_output_index);
593 216136 : }
594 :
595 : template <bool is_ad>
596 : void
597 1472824 : 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 3035428 : for (unsigned int p = 0; p < _num_partitions; ++p)
607 : {
608 3574132 : for (unsigned int t = 0; t < _num_tiles[p]; ++t)
609 : {
610 : // If input is within a specfic tile's window of applicability
611 2011528 : 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 1790948 : if (_tiling[p][in_index] == 1)
616 : {
617 1579476 : if (derivative)
618 286664 : weights[p][t] = 0.0;
619 : }
620 : else
621 : {
622 : // Flag to ensure weights are applied only once
623 : bool overlap = false;
624 968384 : for (unsigned int tt = 0; tt < _num_tiles[p]; ++tt)
625 : {
626 756912 : 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 461620 : 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 84360 : if (_input_limits[p][t][in_index][0] < _input_limits[p][tt][in_index][0] &&
637 37800 : _input_limits[p][t][in_index][1] > _input_limits[p][tt][in_index][0])
638 : {
639 37800 : 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 46560 : else if (_input_limits[p][t][in_index][0] > _input_limits[p][tt][in_index][0] &&
646 46560 : _input_limits[p][t][in_index][0] < _input_limits[p][tt][in_index][1])
647 : {
648 46560 : if (derivative)
649 17280 : 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 29280 : 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 211472 : if (!overlap && derivative)
665 43380 : weights[p][t] = 0.0;
666 : }
667 : }
668 : // If input is below the lower tile limit
669 220580 : else if (input < _input_limits[p][t][in_index][0])
670 : {
671 : // If the lower tile limit equals the lower global limit
672 175194 : if (_input_limits[p][t][in_index][0] == _global_limits[in_index].first)
673 : {
674 91514 : if (_window_failure[in_index].first == WindowFailure::EXTRAPOLATE)
675 : {
676 90560 : if (derivative)
677 0 : weights[p][t] *= -sigmoid(0.0, _input_limits[p][t][in_index][0], input, derivative);
678 : else
679 159080 : weights[p][t] *= (1.0 - sigmoid(0.0, _input_limits[p][t][in_index][0], input));
680 90560 : input = _input_limits[p][t][in_index][0];
681 : }
682 954 : else if (_window_failure[in_index].first == WindowFailure::USELIMIT)
683 0 : input = _input_limits[p][t][in_index][0];
684 : else
685 : {
686 954 : weights[p][t] = 0.0;
687 954 : std::stringstream msg;
688 954 : msg << "In " << _name << ": " << _index_name[in_index]
689 2862 : << " input parameter with value (" << MetaPhysicL::raw_value(input)
690 954 : << ") 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 954 : input = _input_limits[p][t][in_index][0];
694 :
695 954 : if (_window_failure[in_index].first == WindowFailure::WARN)
696 952 : 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 952 : }
703 : }
704 : // if input below tile limit, update weight of tile to be zero
705 : else
706 83680 : weights[p][t] = 0.0;
707 : }
708 : // If input is above the upper tile limit
709 45386 : else if (input > _input_limits[p][t][in_index][1])
710 : {
711 45386 : 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 2 : 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 45384 : 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 1472820 : }
747 :
748 : template <bool is_ad>
749 : bool
750 758745 : LAROMANCEStressUpdateBaseTempl<is_ad>::checkInTile(const unsigned int p, const unsigned int t) const
751 : {
752 3775020 : for (unsigned int i = 0; i < _num_inputs; ++i)
753 3333225 : if (_input_values[i] < _input_limits[p][t][i][0] ||
754 2792515 : _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 461620 : 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 461620 : if (_input_limits[p][t][in_index][0] != _input_limits[p][tt][in_index][0] &&
767 328188 : _input_limits[p][t][in_index][1] != _input_limits[p][tt][in_index][1])
768 : return true;
769 : else
770 133432 : return false;
771 : }
772 :
773 : template <bool is_ad>
774 : GenericReal<is_ad>
775 304132 : 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 305396 : GenericReal<is_ad> trial_stress_mpa = _stress_function
784 304132 : ? _stress_function->value(_t, _q_point[_qp])
785 291748 : : effective_trial_stress * _stress_ucf;
786 71536 : GenericReal<is_ad> dtrial_stress_dscalar = 0.0;
787 :
788 : // Update stress if strain is being applied, i.e. non-testing simulation
789 304132 : if (this->_apply_strain)
790 : {
791 360580 : trial_stress_mpa -= this->_three_shear_modulus * scalar * _stress_ucf;
792 290788 : dtrial_stress_dscalar -= this->_three_shear_modulus * _stress_ucf;
793 : }
794 304132 : _input_values[_stress_input_index] = trial_stress_mpa;
795 :
796 : // Update weights for each partition with new stress
797 625732 : for (unsigned int p = 0; p < _num_partitions; ++p)
798 321600 : _weights[p] = _non_stress_weights[p];
799 304132 : computeTileWeight(_weights, _input_values[_stress_input_index], _stress_input_index);
800 304132 : auto dweights_dstress = _non_stress_weights;
801 304132 : computeTileWeight(
802 304132 : dweights_dstress, _input_values[_stress_input_index], _stress_input_index, true);
803 :
804 304132 : computePartitionWeights(_partition_weights, _dpartition_weight_dstress);
805 :
806 : // Save extrapolation as a material property in order quantify adequate tiling range
807 304132 : _extrapolation[_qp] = 0.0;
808 625732 : for (unsigned int p = 0; p < _num_partitions; ++p)
809 730540 : for (unsigned int t = 0; t < _num_tiles[p]; ++t)
810 480476 : _extrapolation[_qp] += MetaPhysicL::raw_value(_weights[p][t] * _partition_weights[p]);
811 :
812 71536 : GenericReal<is_ad> total_rom_effective_strain_inc = 0.0;
813 71536 : GenericReal<is_ad> dtotal_rom_effective_strain_inc_dstress = 0.0;
814 :
815 : // Run ROM if all values are within windows.
816 625732 : for (unsigned int p = 0; p < _num_partitions; p++)
817 : {
818 321600 : if (_partition_weights[p])
819 : {
820 : // compute weight normalizing factor
821 71536 : GenericReal<is_ad> weight_normalizer = 0;
822 : unsigned int number_of_active_tiles = 0;
823 683240 : for (unsigned int t = 0; t < _num_tiles[p]; ++t)
824 : {
825 : // count number of active tiles
826 371421 : number_of_active_tiles += checkInTile(p, t);
827 378077 : if (_weights[p][t])
828 : {
829 : // tile normalization factor (sum of tile weights)
830 308531 : weight_normalizer += _weights[p][t];
831 : }
832 : }
833 :
834 : // normalize weights only when 3 tiles overlap
835 311819 : if (number_of_active_tiles == 3)
836 7920 : for (unsigned int t = 0; t < _num_tiles[p]; ++t)
837 : {
838 5940 : _weights[p][t] /= weight_normalizer;
839 5940 : dweights_dstress[p][t] /= weight_normalizer;
840 : }
841 :
842 683240 : for (unsigned int t = 0; t < _num_tiles[p]; ++t)
843 378077 : if (_weights[p][t])
844 : {
845 308531 : const GenericReal<is_ad> rom = computeROM(t, p, _strain_output_index);
846 308531 : if (rom == std::numeric_limits<float>::infinity())
847 0 : mooseError("In ", _name, ": Output for strain increment reaches infinity: ", rom);
848 :
849 308531 : total_rom_effective_strain_inc += _partition_weights[p] * _weights[p][t] * rom;
850 :
851 308531 : dtotal_rom_effective_strain_inc_dstress +=
852 308531 : _partition_weights[p] * _weights[p][t] * computeROM(t, p, _strain_output_index, true);
853 308531 : if (_dpartition_weight_dstress[p])
854 13306 : dtotal_rom_effective_strain_inc_dstress +=
855 13306 : _dpartition_weight_dstress[p] * _weights[p][t] * rom;
856 308531 : if (dweights_dstress[p][t])
857 5760 : dtotal_rom_effective_strain_inc_dstress +=
858 5760 : _partition_weights[p] * dweights_dstress[p][t] * rom;
859 : }
860 : }
861 : }
862 :
863 304132 : if (_verbose)
864 : {
865 : Moose::err << std::setprecision(9);
866 0 : GenericReal<is_ad> environmental = 0.0;
867 9856 : if (_environmental)
868 0 : environmental = (*_environmental)[_qp];
869 9856 : Moose::err << "Verbose information from " << _name << ": \n";
870 9856 : Moose::err << " dt: " << _dt << "\n";
871 9856 : Moose::err << " old cell disl: " << _old_input_values[_cell_output_index] << "\n";
872 9856 : Moose::err << " old wall disl: " << _old_input_values[_wall_output_index] << "\n";
873 : Moose::err << " initial stress (MPa): "
874 9856 : << MetaPhysicL::raw_value(effective_trial_stress) * _stress_ucf << "\n";
875 9856 : 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 9856 : Moose::err << " old effective strain: " << _old_input_values[_strain_output_index] << "\n";
881 9856 : Moose::err << " extrapolation: " << MetaPhysicL::raw_value(_extrapolation[_qp]) << "\n";
882 9856 : Moose::err << " partition 2 weight: " << MetaPhysicL::raw_value(_second_partition_weight[_qp])
883 : << "\n";
884 : Moose::err << " weights by tile, partition 1: ";
885 49280 : 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 9856 : if (_num_partitions > 1)
889 : {
890 : Moose::err << " weights by tile, partition 2: ";
891 39424 : 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 49280 : 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 9856 : if (_num_partitions > 1)
900 : {
901 : Moose::err << " nonstress weights by tile, partition 2: ";
902 39424 : 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 375668 : _creep_rate[_qp] = total_rom_effective_strain_inc / _dt;
912 304132 : _derivative = dtotal_rom_effective_strain_inc_dstress * dtrial_stress_dscalar - 1.0;
913 :
914 304132 : if (!this->_apply_strain)
915 : {
916 13344 : if (_verbose)
917 : Moose::err << " Strain not applied due to apply_strain input parameter!" << std::endl;
918 13344 : _derivative = 1.0;
919 13344 : return 0.0;
920 : }
921 220996 : return total_rom_effective_strain_inc - scalar;
922 304132 : }
923 :
924 : template <bool is_ad>
925 : void
926 598792 : LAROMANCEStressUpdateBaseTempl<is_ad>::precomputeROM(const unsigned out_index)
927 : {
928 1238560 : for (unsigned int p = 0; p < _num_partitions; ++p)
929 : {
930 1484416 : for (unsigned int t = 0; t < _num_tiles[p]; ++t)
931 : {
932 : // Only precompute for tiles that don't have zero weight
933 844648 : if (_non_stress_weights[p][t])
934 : {
935 4371792 : for (unsigned int i = 0; i < _num_inputs; ++i)
936 : {
937 3643160 : if (i != _stress_input_index)
938 : {
939 2914528 : _rom_inputs[p][t][i] =
940 2914528 : 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 2914528 : buildPolynomials(p, _rom_inputs[p][t][i], _polynomial_inputs[p][t][i]);
945 : }
946 : }
947 728632 : precomputeValues(
948 : p, _coefs[p][t][out_index], _polynomial_inputs[p][t], _precomputed_vals[p][t]);
949 : }
950 : }
951 : }
952 598792 : }
953 :
954 : template <bool is_ad>
955 : GenericReal<is_ad>
956 1001974 : 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 1001974 : _rom_inputs[p][t][_stress_input_index] =
963 1001974 : 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 1001974 : _transformed_normalization_limits[p][t][out_index][_stress_input_index]);
967 :
968 1001974 : buildPolynomials(
969 : p, _rom_inputs[p][t][_stress_input_index], _polynomial_inputs[p][t][_stress_input_index]);
970 :
971 : // Compute ROM values
972 1001974 : const GenericReal<is_ad> rom_output =
973 757974 : computeValues(p, _precomputed_vals[p][t], _polynomial_inputs[p][t]);
974 :
975 : // Return converted output if not derivative
976 1001974 : if (!derivative)
977 693443 : return convertOutput(p, _old_input_values, rom_output, out_index);
978 :
979 308531 : const GenericReal<is_ad> drom_input =
980 308531 : 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 308531 : _transformed_normalization_limits[p][t][out_index][_stress_input_index],
984 : derivative);
985 :
986 308531 : std::vector<GenericReal<is_ad>> dpolynomial_inputs(_degree[p], 0.0);
987 308531 : buildPolynomials(
988 308531 : p, _rom_inputs[p][t][_stress_input_index], dpolynomial_inputs, drom_input, derivative);
989 :
990 308531 : const GenericReal<is_ad> drom_output = computeValues(
991 : p, _precomputed_vals[p][t], _polynomial_inputs[p][t], dpolynomial_inputs, derivative);
992 308531 : return convertOutput(p, _old_input_values, rom_output, out_index, drom_output, derivative);
993 308531 : }
994 :
995 : template <bool is_ad>
996 : GenericReal<is_ad>
997 4225033 : 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 4225033 : GenericReal<is_ad> x = input;
1004 4225033 : convertValue(x, transform, transform_coef, derivative);
1005 :
1006 : // transformed_limits[2] = transformed_limits[1] - transformed_limits[0]
1007 4225033 : if (derivative)
1008 243651 : return x / transformed_limits[2];
1009 : else
1010 3916502 : return (x - transformed_limits[0]) / transformed_limits[2] - 1.0;
1011 : }
1012 :
1013 : template <bool is_ad>
1014 : void
1015 4225033 : 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 18675410 : for (unsigned int d = 0; d < _degree[p]; ++d)
1023 : {
1024 14450377 : polynomial_inputs[d] = computePolynomial(rom_input, d);
1025 :
1026 14450377 : if (derivative)
1027 1405219 : polynomial_inputs[d] = drom_input * computePolynomial(rom_input, d, derivative);
1028 : }
1029 4225033 : }
1030 :
1031 : template <bool is_ad>
1032 : void
1033 728632 : 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 572103032 : for (unsigned int c = 0; c < _num_coefs[p]; ++c)
1040 : {
1041 571374400 : precomputed[c] = coefs[c];
1042 3428246400 : for (unsigned int i = 0; i < _num_inputs; ++i)
1043 2856872000 : if (i != _stress_input_index)
1044 2285497600 : precomputed[c] *= polynomial_inputs[i][_makeframe_helper[p][c + _num_coefs[p] * i]];
1045 : }
1046 728632 : }
1047 :
1048 : template <bool is_ad>
1049 : GenericReal<is_ad>
1050 1310505 : 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 308880 : GenericReal<is_ad> rom_output = 0.0;
1058 1206880242 : for (unsigned int c = 0; c < _num_coefs[p]; ++c)
1059 : {
1060 1205569737 : if (!derivative)
1061 919786918 : rom_output +=
1062 919786918 : precomputed[c] *
1063 919786918 : polynomial_inputs[_stress_input_index]
1064 919786918 : [_makeframe_helper[p][c + _num_coefs[p] * _stress_input_index]];
1065 :
1066 : else
1067 285782819 : rom_output +=
1068 285782819 : precomputed[c] *
1069 285782819 : dpolynomial_inputs[_makeframe_helper[p][c + _num_coefs[p] * _stress_input_index]];
1070 : }
1071 1310505 : return rom_output;
1072 : }
1073 :
1074 : template <bool is_ad>
1075 : GenericReal<is_ad>
1076 1001974 : 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 : using std::exp;
1084 1001974 : if (out_index == _strain_output_index)
1085 : {
1086 617062 : if (derivative)
1087 373411 : return exp(rom_output) * _dt * drom_output;
1088 : else
1089 373411 : return exp(rom_output) * _dt;
1090 : }
1091 :
1092 384912 : if (derivative)
1093 0 : return 0.0;
1094 :
1095 384912 : GenericReal<is_ad> expout = exp(rom_output);
1096 : mooseAssert(expout > 0.0, "ROM calculated strain increment must be positive");
1097 :
1098 384912 : if (expout > _cutoff[p])
1099 270672 : expout -= _cutoff[p];
1100 : else
1101 0 : expout = -_cutoff[p] * _cutoff[p] / expout + _cutoff[p];
1102 :
1103 499152 : return -expout * old_input_values[out_index] * _dt;
1104 : }
1105 :
1106 : template <bool is_ad>
1107 : GenericReal<is_ad>
1108 15596076 : LAROMANCEStressUpdateBaseTempl<is_ad>::computePolynomial(const GenericReal<is_ad> & value,
1109 : const unsigned int degree,
1110 : const bool derivative)
1111 : {
1112 15596076 : if (degree == 0)
1113 : {
1114 4533564 : if (derivative)
1115 64880 : return 0.0;
1116 4225033 : return 1.0;
1117 : }
1118 11062512 : else if (degree == 1)
1119 : {
1120 3687504 : if (derivative)
1121 64880 : return 1.0;
1122 3408448 : return value;
1123 : }
1124 7375008 : else if (degree == 2)
1125 : {
1126 3687504 : if (derivative)
1127 343936 : return 3.0 * value;
1128 4564720 : return 1.5 * Utility::pow<2>(value) - 0.5;
1129 : }
1130 : else
1131 : {
1132 3687504 : if (derivative)
1133 343936 : return 7.5 * Utility::pow<2>(value) - 1.5;
1134 5720992 : return 2.5 * Utility::pow<3>(value) - 1.5 * value;
1135 : }
1136 : }
1137 :
1138 : template <bool is_ad>
1139 : std::vector<std::vector<std::vector<std::vector<Real>>>>
1140 642 : LAROMANCEStressUpdateBaseTempl<is_ad>::getTransformedLimits(
1141 : const unsigned int p, const std::vector<std::vector<std::vector<Real>>> limits)
1142 : {
1143 642 : std::vector<std::vector<std::vector<std::vector<Real>>>> transformed_limits(
1144 : _num_tiles[p],
1145 642 : std::vector<std::vector<std::vector<Real>>>(
1146 1284 : _num_outputs, std::vector<std::vector<Real>>(_num_inputs, std::vector<Real>(3))));
1147 1914 : for (unsigned int t = 0; t < _num_tiles[p]; ++t)
1148 : {
1149 5088 : for (unsigned int o = 0; o < _num_outputs; ++o)
1150 : {
1151 22896 : for (unsigned int i = 0; i < _num_inputs; ++i)
1152 : {
1153 57240 : for (unsigned int k = 0; k < 2; ++k)
1154 : {
1155 38160 : transformed_limits[t][o][i][k] = limits[t][i][k];
1156 38160 : convertValue(
1157 : transformed_limits[t][o][i][k], _transform[p][t][o][i], _transform_coefs[p][t][o][i]);
1158 : }
1159 19080 : transformed_limits[t][o][i][2] =
1160 19080 : (transformed_limits[t][o][i][1] - transformed_limits[t][o][i][0]) / 2.0;
1161 : }
1162 : }
1163 : }
1164 642 : return transformed_limits;
1165 : }
1166 :
1167 : template <bool is_ad>
1168 : std::vector<unsigned int>
1169 642 : LAROMANCEStressUpdateBaseTempl<is_ad>::getMakeFrameHelper(const unsigned int p) const
1170 : {
1171 642 : std::vector<unsigned int> v(_num_coefs[p] * _num_inputs);
1172 :
1173 400254 : for (unsigned int numcoeffs = 0; numcoeffs < _num_coefs[p]; ++numcoeffs)
1174 2397672 : for (unsigned int invar = 0; invar < _num_inputs; ++invar)
1175 1998060 : v[numcoeffs + _num_coefs[p] * invar] =
1176 1998060 : numcoeffs / MathUtils::pow(_degree[p], invar) % _degree[p];
1177 :
1178 642 : return v;
1179 : }
1180 :
1181 : template <bool is_ad>
1182 : void
1183 216136 : LAROMANCEStressUpdateBaseTempl<is_ad>::computeStressFinalize(
1184 : const GenericRankTwoTensor<is_ad> & plastic_strain_increment)
1185 : {
1186 216136 : _cell_dislocation_increment = 0.0;
1187 216136 : _wall_dislocation_increment = 0.0;
1188 236536 : if (_input_values[_stress_input_index])
1189 : {
1190 191328 : precomputeROM(_cell_output_index);
1191 396288 : for (unsigned int p = 0; p < _num_partitions; ++p)
1192 205000 : if (_partition_weights[p])
1193 434000 : for (unsigned int t = 0; t < _num_tiles[p]; ++t)
1194 244912 : if (_weights[p][t])
1195 192456 : _cell_dislocation_increment +=
1196 192456 : _partition_weights[p] * _weights[p][t] * computeROM(t, p, _cell_output_index);
1197 191328 : precomputeROM(_wall_output_index);
1198 396288 : for (unsigned int p = 0; p < _num_partitions; ++p)
1199 205000 : if (_partition_weights[p])
1200 434000 : for (unsigned int t = 0; t < _num_tiles[p]; ++t)
1201 244912 : if (_weights[p][t])
1202 192456 : _wall_dislocation_increment +=
1203 192456 : _partition_weights[p] * _weights[p][t] * computeROM(t, p, _wall_output_index);
1204 : }
1205 :
1206 300352 : _cell_rate[_qp] = _cell_dislocation_increment / _dt;
1207 216136 : _wall_rate[_qp] = _wall_dislocation_increment / _dt;
1208 300352 : _cell_dislocations[_qp] = _old_input_values[_cell_output_index] + _cell_dislocation_increment;
1209 300352 : _wall_dislocations[_qp] = _old_input_values[_wall_output_index] + _wall_dislocation_increment;
1210 :
1211 : // For (possibly) substepping.
1212 216136 : _plastic_strain_increment += MetaPhysicL::raw_value(plastic_strain_increment);
1213 :
1214 : // Prevent the ROM from calculating and proceeding with negative dislocations
1215 216136 : if (_apply_strain && (_cell_dislocations[_qp] < 0.0 || _wall_dislocations[_qp] < 0.0))
1216 0 : mooseException("In ",
1217 : _name,
1218 : ": Negative disclocation density calculated for cell (old : ",
1219 : MetaPhysicL::raw_value(_old_input_values[_cell_output_index]),
1220 : " increment: ",
1221 : MetaPhysicL::raw_value(_cell_dislocation_increment),
1222 : " value: ",
1223 : MetaPhysicL::raw_value(_cell_dislocations[_qp]),
1224 : ") or wall (old : ",
1225 : MetaPhysicL::raw_value(_old_input_values[_wall_output_index]),
1226 : " increment: ",
1227 : MetaPhysicL::raw_value(_wall_dislocation_increment),
1228 : " value: ",
1229 : MetaPhysicL::raw_value(_wall_dislocations[_qp]),
1230 : ").");
1231 :
1232 216136 : if (_verbose)
1233 : {
1234 : Moose::err << " Finalized ROM output\n";
1235 : Moose::err << " effective creep strain increment: "
1236 10016 : << std::sqrt(2.0 / 3.0 *
1237 10016 : MetaPhysicL::raw_value(_plastic_strain_increment.doubleContraction(
1238 : _plastic_strain_increment)))
1239 : << "\n";
1240 : Moose::err << " total effective creep strain: "
1241 10016 : << std::sqrt(2.0 / 3.0 *
1242 10016 : MetaPhysicL::raw_value(this->_creep_strain[_qp].doubleContraction(
1243 : this->_creep_strain[_qp])))
1244 : << "\n";
1245 10016 : Moose::err << " creep rate: " << MetaPhysicL::raw_value(_creep_rate[_qp]) << "\n";
1246 10016 : Moose::err << " cell dislocation rate: " << MetaPhysicL::raw_value(_cell_rate[_qp]) << "\n";
1247 10016 : Moose::err << " wall dislocation rate: " << MetaPhysicL::raw_value(_wall_rate[_qp]) << "\n";
1248 10016 : Moose::err << " new cell dislocations: " << MetaPhysicL::raw_value(_cell_dislocations[_qp])
1249 : << "\n";
1250 10016 : Moose::err << " new wall dislocations: " << MetaPhysicL::raw_value(_wall_dislocations[_qp])
1251 : << "\n"
1252 : << std::endl;
1253 : }
1254 :
1255 300352 : RadialReturnCreepStressUpdateBaseTempl<is_ad>::computeStressFinalize(
1256 : MetaPhysicL::raw_value(_plastic_strain_increment));
1257 216136 : }
1258 :
1259 : template <bool is_ad>
1260 : Real
1261 193352 : LAROMANCEStressUpdateBaseTempl<is_ad>::computeTimeStepLimit()
1262 : {
1263 193352 : Real limited_dt = RadialReturnCreepStressUpdateBaseTempl<is_ad>::computeTimeStepLimit();
1264 :
1265 : Real cell_strain_inc = std::abs(MetaPhysicL::raw_value(_cell_dislocation_increment));
1266 193352 : if (cell_strain_inc && _old_input_values[_cell_output_index])
1267 161648 : limited_dt = std::min(limited_dt,
1268 306456 : _dt * _max_cell_increment * _old_input_values[_cell_output_index] /
1269 : cell_strain_inc);
1270 : Real wall_strain_inc = std::abs(MetaPhysicL::raw_value(_wall_dislocation_increment));
1271 193352 : if (wall_strain_inc && _old_input_values[_wall_output_index])
1272 161648 : limited_dt = std::min(limited_dt,
1273 162068 : _dt * _max_wall_increment * _old_input_values[_wall_output_index] /
1274 : wall_strain_inc);
1275 :
1276 193352 : return limited_dt;
1277 : }
1278 :
1279 : template <bool is_ad>
1280 : GenericReal<is_ad>
1281 244792 : LAROMANCEStressUpdateBaseTempl<is_ad>::sigmoid(const Real lower,
1282 : const Real upper,
1283 : const GenericReal<is_ad> & val,
1284 : const bool derivative)
1285 : {
1286 : using std::exp;
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 244792 : const GenericReal<is_ad> x = 2.0 * (val - lower) / (upper - lower) - 1.0;
1294 :
1295 244792 : if (x == -1.0)
1296 : {
1297 36682 : if (derivative)
1298 0 : return 0.0;
1299 : else
1300 17781 : return 1.0;
1301 : }
1302 208110 : else if (x == 1.0)
1303 : {
1304 0 : if (derivative)
1305 0 : return 0.0;
1306 : else
1307 0 : return 0.0;
1308 : }
1309 175388 : else if (x < 1.0 && x > -1.0)
1310 : {
1311 312428 : const GenericReal<is_ad> plus = exp(-2.0 / (1.0 + x));
1312 312428 : const GenericReal<is_ad> minus = exp(-2.0 / (1.0 - x));
1313 :
1314 175388 : if (!derivative)
1315 212694 : return 1.0 - plus / (plus + minus);
1316 :
1317 31214 : const GenericReal<is_ad> dplus_dx = plus * 2.0 / Utility::pow<2>(1.0 + x);
1318 31214 : const GenericReal<is_ad> dminus_dx = -minus * 2.0 / Utility::pow<2>(1.0 - x);
1319 :
1320 31214 : return (plus * dminus_dx - dplus_dx * minus) / Utility::pow<2>(plus + minus) * 2.0 /
1321 31214 : (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 202535 : 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 202535 : SingleVariableReturnMappingSolutionTempl<is_ad>::outputIterationStep(
1387 : iter_output, effective_trial_stress, scalar, reference_residual);
1388 202535 : if (iter_output)
1389 0 : *iter_output << " derivative: "
1390 0 : << MetaPhysicL::raw_value(computeDerivative(effective_trial_stress, scalar))
1391 : << std::endl;
1392 202535 : }
1393 :
1394 : template <bool is_ad>
1395 : void
1396 609 : LAROMANCEStressUpdateBaseTempl<is_ad>::checkJSONKey(const std::string & key)
1397 : {
1398 1218 : if (!this->isParamValid("model"))
1399 0 : this->paramError("model", "Specify a JSON data filename.");
1400 :
1401 609 : const auto model_file_name = this->template getParam<DataFileName>("model");
1402 609 : if (_json.empty())
1403 0 : this->paramError("model", "The specified JSON data file '", model_file_name, "' is empty.");
1404 609 : if (!_json.contains(key))
1405 0 : this->paramError(
1406 : "model", "The key '", key, "' is missing from the JSON data file '", model_file_name, "'.");
1407 609 : }
1408 :
1409 : template class LAROMANCEStressUpdateBaseTempl<false>;
1410 : template class LAROMANCEStressUpdateBaseTempl<true>;
|
