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 290803 : LAROMANCEStressUpdateBaseTempl<is_ad>::substeppingCapabilityEnabled()
298 : {
299 290803 : 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 8568 : LAROMANCEStressUpdateBaseTempl<is_ad>::initQpStatefulProperties()
499 : {
500 : MooseRandom rng;
501 17136 : rng.seed(0, this->template getParam<unsigned int>("seed"));
502 :
503 8568 : _cell_dislocations[_qp] = rng.randNormal(
504 17136 : this->template getParam<Real>("initial_cell_dislocation_density"),
505 8568 : this->template getParam<Real>("initial_cell_dislocation_density") *
506 17136 : this->template getParam<Real>("cell_dislocations_normal_distribution_width"));
507 8568 : _wall_dislocations[_qp] = rng.randNormal(
508 17136 : this->template getParam<Real>("initial_wall_dislocation_density"),
509 8568 : this->template getParam<Real>("initial_wall_dislocation_density") *
510 17136 : this->template getParam<Real>("wall_dislocations_normal_distribution_width"));
511 :
512 8568 : RadialReturnCreepStressUpdateBaseTempl<is_ad>::initQpStatefulProperties();
513 8568 : }
514 :
515 : template <bool is_ad>
516 : GenericReal<is_ad>
517 289424 : 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 289424 : return effective_trial_stress / this->_three_shear_modulus * 0.999999;
523 : }
524 :
525 : template <bool is_ad>
526 : void
527 290764 : LAROMANCEStressUpdateBaseTempl<is_ad>::resetIncrementalMaterialProperties()
528 : {
529 290764 : _cell_dislocation_increment = 0.0;
530 290764 : _wall_dislocation_increment = 0.0;
531 :
532 : _plastic_strain_increment.zero();
533 :
534 290764 : _wall_dislocations_step[_qp] = 0.0;
535 290764 : _cell_dislocations_step[_qp] = 0.0;
536 290764 : }
537 :
538 : template <bool is_ad>
539 : void
540 62652 : LAROMANCEStressUpdateBaseTempl<is_ad>::storeIncrementalMaterialProperties(
541 : const unsigned int total_number_of_substeps)
542 : {
543 62652 : _wall_dislocations_step[_qp] += _wall_dislocation_increment;
544 62652 : _cell_dislocations_step[_qp] += _cell_dislocation_increment;
545 62652 : _number_of_substeps[_qp] = total_number_of_substeps;
546 62652 : }
547 :
548 : template <bool is_ad>
549 : void
550 329828 : LAROMANCEStressUpdateBaseTempl<is_ad>::computeStressInitialize(
551 : const GenericReal<is_ad> & effective_trial_stress,
552 : const GenericRankFourTensor<is_ad> & elasticity_tensor)
553 : {
554 133880 : RadialReturnCreepStressUpdateBaseTempl<is_ad>::computeStressInitialize(effective_trial_stress,
555 : elasticity_tensor);
556 : // Previous substep creep strain
557 329828 : RankTwoTensor creep_strain_substep = this->_creep_strain_old[_qp] + _plastic_strain_increment;
558 :
559 : // Prepare old values
560 329828 : _old_input_values[_cell_output_index] =
561 329828 : _cell_function
562 329828 : ? _cell_function->value(_t, _q_point[_qp])
563 317204 : : (_cell_dislocations_old[_qp] + MetaPhysicL::raw_value(_cell_dislocations_step[_qp]));
564 329828 : _old_input_values[_wall_output_index] =
565 329828 : _wall_function
566 329828 : ? _wall_function->value(_t, _q_point[_qp])
567 317204 : : (_wall_dislocations_old[_qp] + MetaPhysicL::raw_value(_wall_dislocations_step[_qp]));
568 329828 : _old_input_values[_strain_output_index] =
569 329828 : _creep_strain_old_forcing_function
570 329828 : ? _creep_strain_old_forcing_function->value(_t, _q_point[_qp])
571 : : MetaPhysicL::raw_value(
572 317204 : std::sqrt((creep_strain_substep).doubleContraction(creep_strain_substep) / 1.5));
573 :
574 : // Prepare input
575 329828 : _input_values[_cell_input_index] = _old_input_values[_cell_output_index];
576 329828 : _input_values[_wall_input_index] = _old_input_values[_wall_output_index];
577 647032 : _input_values[_stress_input_index] = _stress_function ? _stress_function->value(_t, _q_point[_qp])
578 317204 : : effective_trial_stress * _stress_ucf;
579 329828 : _input_values[_old_strain_input_index] = _old_input_values[_strain_output_index];
580 329828 : _input_values[_temperature_input_index] = _temperature[_qp];
581 329828 : if (_environmental)
582 0 : _input_values[_environmental_input_index] = (*_environmental)[_qp];
583 :
584 : // Determine tile weight and check if input is in range
585 673368 : for (unsigned int p = 0; p < _num_partitions; ++p)
586 343540 : std::fill(_non_stress_weights[p].begin(), _non_stress_weights[p].end(), 1.0);
587 1978964 : for (unsigned int i = 0; i < _num_inputs; i++)
588 1649140 : if (i != _stress_input_index)
589 1319312 : computeTileWeight(_non_stress_weights, _input_values[i], i);
590 :
591 : // Precompute transformed input and prebuild polynomials for inputs other than strain
592 329824 : precomputeROM(_strain_output_index);
593 329824 : }
594 :
595 : template <bool is_ad>
596 : void
597 2252704 : 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 4595188 : for (unsigned int p = 0; p < _num_partitions; ++p)
607 : {
608 5133892 : for (unsigned int t = 0; t < _num_tiles[p]; ++t)
609 : {
610 : // If input is within a specfic tile's window of applicability
611 2791408 : 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 2525312 : if (_tiling[p][in_index] == 1)
616 : {
617 2313840 : if (derivative)
618 449228 : 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 266096 : else if (input < _input_limits[p][t][in_index][0])
670 : {
671 : // If the lower tile limit equals the lower global limit
672 220710 : if (_input_limits[p][t][in_index][0] == _global_limits[in_index].first)
673 : {
674 137030 : if (_window_failure[in_index].first == WindowFailure::EXTRAPOLATE)
675 : {
676 135464 : if (derivative)
677 0 : weights[p][t] *= -sigmoid(0.0, _input_limits[p][t][in_index][0], input, derivative);
678 : else
679 242872 : weights[p][t] *= (1.0 - sigmoid(0.0, _input_limits[p][t][in_index][0], input));
680 135464 : input = _input_limits[p][t][in_index][0];
681 : }
682 1566 : else if (_window_failure[in_index].first == WindowFailure::USELIMIT)
683 0 : input = _input_limits[p][t][in_index][0];
684 : else
685 : {
686 1566 : weights[p][t] = 0.0;
687 1566 : std::stringstream msg;
688 1566 : msg << "In " << _name << ": " << _index_name[in_index]
689 4698 : << " input parameter with value (" << MetaPhysicL::raw_value(input)
690 1566 : << ") 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 1566 : input = _input_limits[p][t][in_index][0];
694 :
695 1566 : if (_window_failure[in_index].first == WindowFailure::WARN)
696 1564 : 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 1564 : }
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 2252700 : }
747 :
748 : template <bool is_ad>
749 : bool
750 921309 : LAROMANCEStressUpdateBaseTempl<is_ad>::checkInTile(const unsigned int p, const unsigned int t) const
751 : {
752 4750404 : for (unsigned int i = 0; i < _num_inputs; ++i)
753 4146045 : if (_input_values[i] < _input_limits[p][t][i][0] ||
754 3408775 : _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 466696 : 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 467960 : GenericReal<is_ad> trial_stress_mpa = _stress_function
784 466696 : ? _stress_function->value(_t, _q_point[_qp])
785 454312 : : effective_trial_stress * _stress_ucf;
786 110848 : GenericReal<is_ad> dtrial_stress_dscalar = 0.0;
787 :
788 : // Update stress if strain is being applied, i.e. non-testing simulation
789 466696 : if (this->_apply_strain)
790 : {
791 561160 : trial_stress_mpa -= this->_three_shear_modulus * scalar * _stress_ucf;
792 452488 : dtrial_stress_dscalar -= this->_three_shear_modulus * _stress_ucf;
793 : }
794 466696 : _input_values[_stress_input_index] = trial_stress_mpa;
795 :
796 : // Update weights for each partition with new stress
797 950860 : for (unsigned int p = 0; p < _num_partitions; ++p)
798 484164 : _weights[p] = _non_stress_weights[p];
799 466696 : computeTileWeight(_weights, _input_values[_stress_input_index], _stress_input_index);
800 466696 : auto dweights_dstress = _non_stress_weights;
801 466696 : computeTileWeight(
802 466696 : dweights_dstress, _input_values[_stress_input_index], _stress_input_index, true);
803 :
804 466696 : computePartitionWeights(_partition_weights, _dpartition_weight_dstress);
805 :
806 : // Save extrapolation as a material property in order quantify adequate tiling range
807 466696 : _extrapolation[_qp] = 0.0;
808 950860 : for (unsigned int p = 0; p < _num_partitions; ++p)
809 1055668 : for (unsigned int t = 0; t < _num_tiles[p]; ++t)
810 682352 : _extrapolation[_qp] += MetaPhysicL::raw_value(_weights[p][t] * _partition_weights[p]);
811 :
812 110848 : GenericReal<is_ad> total_rom_effective_strain_inc = 0.0;
813 110848 : GenericReal<is_ad> dtotal_rom_effective_strain_inc_dstress = 0.0;
814 :
815 : // Run ROM if all values are within windows.
816 950860 : for (unsigned int p = 0; p < _num_partitions; p++)
817 : {
818 484164 : if (_partition_weights[p])
819 : {
820 : // compute weight normalizing factor
821 110848 : GenericReal<is_ad> weight_normalizer = 0;
822 : unsigned int number_of_active_tiles = 0;
823 1008368 : for (unsigned int t = 0; t < _num_tiles[p]; ++t)
824 : {
825 : // count number of active tiles
826 533985 : number_of_active_tiles += checkInTile(p, t);
827 545505 : if (_weights[p][t])
828 : {
829 : // tile normalization factor (sum of tile weights)
830 465619 : weight_normalizer += _weights[p][t];
831 : }
832 : }
833 :
834 : // normalize weights only when 3 tiles overlap
835 474383 : 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 1008368 : for (unsigned int t = 0; t < _num_tiles[p]; ++t)
843 545505 : if (_weights[p][t])
844 : {
845 465619 : const GenericReal<is_ad> rom = computeROM(t, p, _strain_output_index);
846 465619 : if (rom == std::numeric_limits<float>::infinity())
847 0 : mooseError("In ", _name, ": Output for strain increment reaches infinity: ", rom);
848 :
849 465619 : total_rom_effective_strain_inc += _partition_weights[p] * _weights[p][t] * rom;
850 :
851 465619 : dtotal_rom_effective_strain_inc_dstress +=
852 465619 : _partition_weights[p] * _weights[p][t] * computeROM(t, p, _strain_output_index, true);
853 465619 : if (_dpartition_weight_dstress[p])
854 13306 : dtotal_rom_effective_strain_inc_dstress +=
855 13306 : _dpartition_weight_dstress[p] * _weights[p][t] * rom;
856 465619 : 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 466696 : 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 577544 : _creep_rate[_qp] = total_rom_effective_strain_inc / _dt;
912 466696 : _derivative = dtotal_rom_effective_strain_inc_dstress * dtrial_stress_dscalar - 1.0;
913 :
914 466696 : if (!this->_apply_strain)
915 : {
916 14208 : if (_verbose)
917 : Moose::err << " Strain not applied due to apply_strain input parameter!" << std::endl;
918 14208 : _derivative = 1.0;
919 14208 : return 0.0;
920 : }
921 343816 : return total_rom_effective_strain_inc - scalar;
922 466696 : }
923 :
924 : template <bool is_ad>
925 : void
926 909152 : LAROMANCEStressUpdateBaseTempl<is_ad>::precomputeROM(const unsigned out_index)
927 : {
928 1859280 : for (unsigned int p = 0; p < _num_partitions; ++p)
929 : {
930 2105136 : for (unsigned int t = 0; t < _num_tiles[p]; ++t)
931 : {
932 : // Only precompute for tiles that don't have zero weight
933 1155008 : if (_non_stress_weights[p][t])
934 : {
935 6233952 : for (unsigned int i = 0; i < _num_inputs; ++i)
936 : {
937 5194960 : if (i != _stress_input_index)
938 : {
939 4155968 : _rom_inputs[p][t][i] =
940 4155968 : 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 4155968 : buildPolynomials(p, _rom_inputs[p][t][i], _polynomial_inputs[p][t][i]);
945 : }
946 : }
947 1038992 : precomputeValues(
948 : p, _coefs[p][t][out_index], _polynomial_inputs[p][t], _precomputed_vals[p][t]);
949 : }
950 : }
951 : }
952 909152 : }
953 :
954 : template <bool is_ad>
955 : GenericReal<is_ad>
956 1501870 : 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 1501870 : _rom_inputs[p][t][_stress_input_index] =
963 1501870 : 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 1501870 : _transformed_normalization_limits[p][t][out_index][_stress_input_index]);
967 :
968 1501870 : buildPolynomials(
969 : p, _rom_inputs[p][t][_stress_input_index], _polynomial_inputs[p][t][_stress_input_index]);
970 :
971 : // Compute ROM values
972 1501870 : const GenericReal<is_ad> rom_output =
973 1126350 : computeValues(p, _precomputed_vals[p][t], _polynomial_inputs[p][t]);
974 :
975 : // Return converted output if not derivative
976 1501870 : if (!derivative)
977 1036251 : return convertOutput(p, _old_input_values, rom_output, out_index);
978 :
979 465619 : const GenericReal<is_ad> drom_input =
980 465619 : 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 465619 : _transformed_normalization_limits[p][t][out_index][_stress_input_index],
984 : derivative);
985 :
986 465619 : std::vector<GenericReal<is_ad>> dpolynomial_inputs(_degree[p], 0.0);
987 465619 : buildPolynomials(
988 465619 : p, _rom_inputs[p][t][_stress_input_index], dpolynomial_inputs, drom_input, derivative);
989 :
990 465619 : const GenericReal<is_ad> drom_output = computeValues(
991 : p, _precomputed_vals[p][t], _polynomial_inputs[p][t], dpolynomial_inputs, derivative);
992 465619 : return convertOutput(p, _old_input_values, rom_output, out_index, drom_output, derivative);
993 465619 : }
994 :
995 : template <bool is_ad>
996 : GenericReal<is_ad>
997 6123457 : 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 6123457 : GenericReal<is_ad> x = input;
1004 6123457 : convertValue(x, transform, transform_coef, derivative);
1005 :
1006 : // transformed_limits[2] = transformed_limits[1] - transformed_limits[0]
1007 6123457 : if (derivative)
1008 366291 : return x / transformed_limits[2];
1009 : else
1010 5657838 : return (x - transformed_limits[0]) / transformed_limits[2] - 1.0;
1011 : }
1012 :
1013 : template <bool is_ad>
1014 : void
1015 6123457 : 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 28167530 : for (unsigned int d = 0; d < _degree[p]; ++d)
1023 : {
1024 22044073 : polynomial_inputs[d] = computePolynomial(rom_input, d);
1025 :
1026 22044073 : if (derivative)
1027 2171363 : polynomial_inputs[d] = drom_input * computePolynomial(rom_input, d, derivative);
1028 : }
1029 6123457 : }
1030 :
1031 : template <bool is_ad>
1032 : void
1033 1038992 : 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 890222032 : for (unsigned int c = 0; c < _num_coefs[p]; ++c)
1040 : {
1041 889183040 : precomputed[c] = coefs[c];
1042 5335098240 : for (unsigned int i = 0; i < _num_inputs; ++i)
1043 4445915200 : if (i != _stress_input_index)
1044 3556732160 : precomputed[c] *= polynomial_inputs[i][_makeframe_helper[p][c + _num_coefs[p] * i]];
1045 : }
1046 1038992 : }
1047 :
1048 : template <bool is_ad>
1049 : GenericReal<is_ad>
1050 1967489 : 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 474848 : GenericReal<is_ad> rom_output = 0.0;
1058 1880288842 : for (unsigned int c = 0; c < _num_coefs[p]; ++c)
1059 : {
1060 1878321353 : if (!derivative)
1061 1431680422 : rom_output +=
1062 1431680422 : precomputed[c] *
1063 1431680422 : polynomial_inputs[_stress_input_index]
1064 1431680422 : [_makeframe_helper[p][c + _num_coefs[p] * _stress_input_index]];
1065 :
1066 : else
1067 446640931 : rom_output +=
1068 446640931 : precomputed[c] *
1069 446640931 : dpolynomial_inputs[_makeframe_helper[p][c + _num_coefs[p] * _stress_input_index]];
1070 : }
1071 1967489 : return rom_output;
1072 : }
1073 :
1074 : template <bool is_ad>
1075 : GenericReal<is_ad>
1076 1501870 : 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 1501870 : if (out_index == _strain_output_index)
1084 : {
1085 931238 : if (derivative)
1086 564947 : return std::exp(rom_output) * _dt * drom_output;
1087 : else
1088 564947 : return std::exp(rom_output) * _dt;
1089 : }
1090 :
1091 570632 : if (derivative)
1092 0 : return 0.0;
1093 :
1094 570632 : GenericReal<is_ad> expout = std::exp(rom_output);
1095 : mooseAssert(expout > 0.0, "ROM calculated strain increment must be positive");
1096 :
1097 570632 : if (expout > _cutoff[p])
1098 393768 : expout -= _cutoff[p];
1099 : else
1100 0 : expout = -_cutoff[p] * _cutoff[p] / expout + _cutoff[p];
1101 :
1102 747496 : return -expout * old_input_values[out_index] * _dt;
1103 : }
1104 :
1105 : template <bool is_ad>
1106 : GenericReal<is_ad>
1107 23818124 : LAROMANCEStressUpdateBaseTempl<is_ad>::computePolynomial(const GenericReal<is_ad> & value,
1108 : const unsigned int degree,
1109 : const bool derivative)
1110 : {
1111 23818124 : if (degree == 0)
1112 : {
1113 6589076 : if (derivative)
1114 99328 : return 0.0;
1115 6123457 : return 1.0;
1116 : }
1117 17229048 : else if (degree == 1)
1118 : {
1119 5743016 : if (derivative)
1120 99328 : return 1.0;
1121 5306872 : return value;
1122 : }
1123 11486032 : else if (degree == 2)
1124 : {
1125 5743016 : if (derivative)
1126 535472 : return 3.0 * value;
1127 7117176 : return 1.5 * Utility::pow<2>(value) - 0.5;
1128 : }
1129 : else
1130 : {
1131 5743016 : if (derivative)
1132 535472 : return 7.5 * Utility::pow<2>(value) - 1.5;
1133 8927480 : 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 642 : LAROMANCEStressUpdateBaseTempl<is_ad>::getTransformedLimits(
1140 : const unsigned int p, const std::vector<std::vector<std::vector<Real>>> limits)
1141 : {
1142 642 : std::vector<std::vector<std::vector<std::vector<Real>>>> transformed_limits(
1143 : _num_tiles[p],
1144 642 : std::vector<std::vector<std::vector<Real>>>(
1145 1284 : _num_outputs, std::vector<std::vector<Real>>(_num_inputs, std::vector<Real>(3))));
1146 1914 : for (unsigned int t = 0; t < _num_tiles[p]; ++t)
1147 : {
1148 5088 : for (unsigned int o = 0; o < _num_outputs; ++o)
1149 : {
1150 22896 : for (unsigned int i = 0; i < _num_inputs; ++i)
1151 : {
1152 57240 : for (unsigned int k = 0; k < 2; ++k)
1153 : {
1154 38160 : transformed_limits[t][o][i][k] = limits[t][i][k];
1155 38160 : convertValue(
1156 : transformed_limits[t][o][i][k], _transform[p][t][o][i], _transform_coefs[p][t][o][i]);
1157 : }
1158 19080 : transformed_limits[t][o][i][2] =
1159 19080 : (transformed_limits[t][o][i][1] - transformed_limits[t][o][i][0]) / 2.0;
1160 : }
1161 : }
1162 : }
1163 642 : return transformed_limits;
1164 : }
1165 :
1166 : template <bool is_ad>
1167 : std::vector<unsigned int>
1168 642 : LAROMANCEStressUpdateBaseTempl<is_ad>::getMakeFrameHelper(const unsigned int p) const
1169 : {
1170 642 : std::vector<unsigned int> v(_num_coefs[p] * _num_inputs);
1171 :
1172 400254 : for (unsigned int numcoeffs = 0; numcoeffs < _num_coefs[p]; ++numcoeffs)
1173 2397672 : for (unsigned int invar = 0; invar < _num_inputs; ++invar)
1174 1998060 : v[numcoeffs + _num_coefs[p] * invar] =
1175 1998060 : numcoeffs / MathUtils::pow(_degree[p], invar) % _degree[p];
1176 :
1177 642 : return v;
1178 : }
1179 :
1180 : template <bool is_ad>
1181 : void
1182 329824 : LAROMANCEStressUpdateBaseTempl<is_ad>::computeStressFinalize(
1183 : const GenericRankTwoTensor<is_ad> & plastic_strain_increment)
1184 : {
1185 329824 : _cell_dislocation_increment = 0.0;
1186 329824 : _wall_dislocation_increment = 0.0;
1187 363712 : if (_input_values[_stress_input_index])
1188 : {
1189 289664 : precomputeROM(_cell_output_index);
1190 592960 : for (unsigned int p = 0; p < _num_partitions; ++p)
1191 303336 : if (_partition_weights[p])
1192 630672 : for (unsigned int t = 0; t < _num_tiles[p]; ++t)
1193 348112 : if (_weights[p][t])
1194 285316 : _cell_dislocation_increment +=
1195 285316 : _partition_weights[p] * _weights[p][t] * computeROM(t, p, _cell_output_index);
1196 289664 : precomputeROM(_wall_output_index);
1197 592960 : for (unsigned int p = 0; p < _num_partitions; ++p)
1198 303336 : if (_partition_weights[p])
1199 630672 : for (unsigned int t = 0; t < _num_tiles[p]; ++t)
1200 348112 : if (_weights[p][t])
1201 285316 : _wall_dislocation_increment +=
1202 285316 : _partition_weights[p] * _weights[p][t] * computeROM(t, p, _wall_output_index);
1203 : }
1204 :
1205 463704 : _cell_rate[_qp] = _cell_dislocation_increment / _dt;
1206 329824 : _wall_rate[_qp] = _wall_dislocation_increment / _dt;
1207 463704 : _cell_dislocations[_qp] = _old_input_values[_cell_output_index] + _cell_dislocation_increment;
1208 463704 : _wall_dislocations[_qp] = _old_input_values[_wall_output_index] + _wall_dislocation_increment;
1209 :
1210 : // For (possibly) substepping.
1211 329824 : _plastic_strain_increment += MetaPhysicL::raw_value(plastic_strain_increment);
1212 :
1213 : // Prevent the ROM from calculating and proceeding with negative dislocations
1214 329824 : 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 329824 : if (_verbose)
1232 : {
1233 : Moose::err << " Finalized ROM output\n";
1234 : Moose::err << " effective creep strain increment: "
1235 10016 : << std::sqrt(2.0 / 3.0 *
1236 10016 : MetaPhysicL::raw_value(_plastic_strain_increment.doubleContraction(
1237 : _plastic_strain_increment)))
1238 : << "\n";
1239 : Moose::err << " total effective creep strain: "
1240 10016 : << std::sqrt(2.0 / 3.0 *
1241 10016 : MetaPhysicL::raw_value(this->_creep_strain[_qp].doubleContraction(
1242 : this->_creep_strain[_qp])))
1243 : << "\n";
1244 10016 : Moose::err << " creep rate: " << MetaPhysicL::raw_value(_creep_rate[_qp]) << "\n";
1245 10016 : Moose::err << " cell dislocation rate: " << MetaPhysicL::raw_value(_cell_rate[_qp]) << "\n";
1246 10016 : Moose::err << " wall dislocation rate: " << MetaPhysicL::raw_value(_wall_rate[_qp]) << "\n";
1247 10016 : Moose::err << " new cell dislocations: " << MetaPhysicL::raw_value(_cell_dislocations[_qp])
1248 : << "\n";
1249 10016 : Moose::err << " new wall dislocations: " << MetaPhysicL::raw_value(_wall_dislocations[_qp])
1250 : << "\n"
1251 : << std::endl;
1252 : }
1253 :
1254 463704 : RadialReturnCreepStressUpdateBaseTempl<is_ad>::computeStressFinalize(
1255 : MetaPhysicL::raw_value(_plastic_strain_increment));
1256 329824 : }
1257 :
1258 : template <bool is_ad>
1259 : Real
1260 290756 : LAROMANCEStressUpdateBaseTempl<is_ad>::computeTimeStepLimit()
1261 : {
1262 290756 : Real limited_dt = RadialReturnCreepStressUpdateBaseTempl<is_ad>::computeTimeStepLimit();
1263 :
1264 : Real cell_strain_inc = std::abs(MetaPhysicL::raw_value(_cell_dislocation_increment));
1265 290756 : if (cell_strain_inc && _old_input_values[_cell_output_index])
1266 238836 : limited_dt = std::min(limited_dt,
1267 452720 : _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 290756 : if (wall_strain_inc && _old_input_values[_wall_output_index])
1271 238836 : limited_dt = std::min(limited_dt,
1272 239256 : _dt * _max_wall_increment * _old_input_values[_wall_output_index] /
1273 : wall_strain_inc);
1274 :
1275 290756 : return limited_dt;
1276 : }
1277 :
1278 : template <bool is_ad>
1279 : GenericReal<is_ad>
1280 289696 : 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 289696 : const GenericReal<is_ad> x = 2.0 * (val - lower) / (upper - lower) - 1.0;
1291 :
1292 289696 : if (x == -1.0)
1293 : {
1294 36682 : if (derivative)
1295 0 : return 0.0;
1296 : else
1297 17781 : return 1.0;
1298 : }
1299 253014 : else if (x == 1.0)
1300 : {
1301 0 : if (derivative)
1302 0 : return 0.0;
1303 : else
1304 0 : return 0.0;
1305 : }
1306 220292 : else if (x < 1.0 && x > -1.0)
1307 : {
1308 435108 : const GenericReal<is_ad> plus = std::exp(-2.0 / (1.0 + x));
1309 435108 : const GenericReal<is_ad> minus = std::exp(-2.0 / (1.0 - x));
1310 :
1311 220292 : if (!derivative)
1312 296486 : return 1.0 - plus / (plus + minus);
1313 :
1314 31214 : const GenericReal<is_ad> dplus_dx = plus * 2.0 / Utility::pow<2>(1.0 + x);
1315 31214 : const GenericReal<is_ad> dminus_dx = -minus * 2.0 / Utility::pow<2>(1.0 - x);
1316 :
1317 31214 : return (plus * dminus_dx - dplus_dx * minus) / Utility::pow<2>(plus + minus) * 2.0 /
1318 31214 : (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 308731 : 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 308731 : SingleVariableReturnMappingSolutionTempl<is_ad>::outputIterationStep(
1384 : iter_output, effective_trial_stress, scalar, reference_residual);
1385 308731 : if (iter_output)
1386 0 : *iter_output << " derivative: "
1387 0 : << MetaPhysicL::raw_value(computeDerivative(effective_trial_stress, scalar))
1388 : << std::endl;
1389 308731 : }
1390 :
1391 : template <bool is_ad>
1392 : void
1393 609 : LAROMANCEStressUpdateBaseTempl<is_ad>::checkJSONKey(const std::string & key)
1394 : {
1395 1218 : if (!this->isParamValid("model"))
1396 0 : this->paramError("model", "Specify a JSON data filename.");
1397 :
1398 609 : const auto model_file_name = this->template getParam<DataFileName>("model");
1399 609 : if (_json.empty())
1400 0 : this->paramError("model", "The specified JSON data file '", model_file_name, "' is empty.");
1401 609 : if (!_json.contains(key))
1402 0 : this->paramError(
1403 : "model", "The key '", key, "' is missing from the JSON data file '", model_file_name, "'.");
1404 609 : }
1405 :
1406 : template class LAROMANCEStressUpdateBaseTempl<false>;
1407 : template class LAROMANCEStressUpdateBaseTempl<true>;
|
