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 "FiniteStrainCrystalPlasticity.h"
11 : #include "petscblaslapack.h"
12 : #include "libmesh/utility.h"
13 :
14 : #include <fstream>
15 : #include <cmath>
16 :
17 : registerMooseObjectDeprecated("SolidMechanicsApp",
18 : FiniteStrainCrystalPlasticity,
19 : "11/15/2024 12:00");
20 :
21 : InputParameters
22 280 : FiniteStrainCrystalPlasticity::validParams()
23 : {
24 280 : InputParameters params = ComputeStressBase::validParams();
25 280 : params.addClassDescription("Deprecated class: please use CrystalPlasticityKalidindiUpdate and "
26 : "ComputeMultipleCrystalPlasticityStress instead. Crystal Plasticity "
27 : "base class: FCC system with power law flow rule implemented");
28 560 : params.addRequiredParam<int>("nss", "Number of slip systems");
29 560 : params.addParam<std::vector<Real>>("gprops", {}, "Initial values of slip system resistances");
30 560 : params.addParam<std::vector<Real>>("hprops", {}, "Hardening properties");
31 560 : params.addParam<std::vector<Real>>("flowprops", {}, "Parameters used in slip rate equations");
32 560 : params.addRequiredParam<FileName>("slip_sys_file_name",
33 : "Name of the file containing the slip system");
34 560 : params.addParam<FileName>(
35 : "slip_sys_res_prop_file_name",
36 : "",
37 : "Name of the file containing the initial values of slip system resistances");
38 560 : params.addParam<FileName>(
39 : "slip_sys_flow_prop_file_name",
40 : "",
41 : "Name of the file containing the values of slip rate equation parameters");
42 560 : params.addParam<FileName>(
43 : "slip_sys_hard_prop_file_name",
44 : "",
45 : "Name of the file containing the values of hardness evolution parameters");
46 560 : params.addParam<Real>("rtol", 1e-6, "Constitutive stress residue relative tolerance");
47 560 : params.addParam<Real>("abs_tol", 1e-6, "Constitutive stress residue absolute tolerance");
48 560 : params.addParam<Real>("gtol", 1e2, "Constitutive slip system resistance residual tolerance");
49 560 : params.addParam<Real>("slip_incr_tol", 2e-2, "Maximum allowable slip in an increment");
50 560 : params.addParam<unsigned int>("maxiter", 100, "Maximum number of iterations for stress update");
51 560 : params.addParam<unsigned int>(
52 560 : "maxitergss", 100, "Maximum number of iterations for slip system resistance update");
53 560 : params.addParam<unsigned int>(
54 : "num_slip_sys_flowrate_props",
55 560 : 2,
56 : "Number of flow rate properties for a slip system"); // Used for reading flow rate parameters
57 560 : params.addParam<UserObjectName>("read_prop_user_object",
58 : "The ElementReadPropertyFile "
59 : "GeneralUserObject to read element "
60 : "specific property values from file");
61 560 : MooseEnum tan_mod_options("exact none", "none"); // Type of read
62 560 : params.addParam<MooseEnum>("tan_mod_type",
63 : tan_mod_options,
64 : "Type of tangent moduli for preconditioner: default elastic");
65 560 : MooseEnum intvar_read_options("slip_sys_file slip_sys_res_file none", "none");
66 560 : params.addParam<MooseEnum>(
67 : "intvar_read_type",
68 : intvar_read_options,
69 : "Read from options for initial value of internal variables: Default from .i file");
70 560 : params.addParam<unsigned int>("num_slip_sys_props",
71 560 : 0,
72 : "Number of slip system specific properties provided in the file "
73 : "containing slip system normals and directions");
74 560 : params.addParam<bool>(
75 : "gen_random_stress_flag",
76 560 : false,
77 : "Flag to generate random stress to perform time cutback on constitutive failure");
78 560 : params.addParam<bool>("input_random_scaling_var",
79 560 : false,
80 : "Flag to input scaling variable: _Cijkl(0,0,0,0) when false");
81 560 : params.addParam<Real>("random_scaling_var",
82 560 : 1e9,
83 : "Random scaling variable: Large value can cause non-positive definiteness");
84 560 : params.addParam<unsigned int>(
85 : "random_seed",
86 560 : 2000,
87 : "Random integer used to generate random stress when constitutive failure occurs");
88 560 : params.addParam<unsigned int>(
89 560 : "maximum_substep_iteration", 1, "Maximum number of substep iteration");
90 560 : params.addParam<bool>("use_line_search", false, "Use line search in constitutive update");
91 560 : params.addParam<Real>("min_line_search_step_size", 0.01, "Minimum line search step size");
92 560 : params.addParam<Real>("line_search_tol", 0.5, "Line search bisection method tolerance");
93 560 : params.addParam<unsigned int>(
94 560 : "line_search_maxiter", 20, "Line search bisection method maximum number of iteration");
95 560 : MooseEnum line_search_method("CUT_HALF BISECTION", "CUT_HALF");
96 560 : params.addParam<MooseEnum>(
97 : "line_search_method", line_search_method, "The method used in line search");
98 :
99 280 : return params;
100 280 : }
101 :
102 210 : FiniteStrainCrystalPlasticity::FiniteStrainCrystalPlasticity(const InputParameters & parameters)
103 : : ComputeStressBase(parameters),
104 210 : _nss(getParam<int>("nss")),
105 420 : _gprops(getParam<std::vector<Real>>("gprops")),
106 420 : _hprops(getParam<std::vector<Real>>("hprops")),
107 420 : _flowprops(getParam<std::vector<Real>>("flowprops")),
108 420 : _slip_sys_file_name(getParam<FileName>("slip_sys_file_name")),
109 420 : _slip_sys_res_prop_file_name(getParam<FileName>("slip_sys_res_prop_file_name")),
110 420 : _slip_sys_flow_prop_file_name(getParam<FileName>("slip_sys_flow_prop_file_name")),
111 420 : _slip_sys_hard_prop_file_name(getParam<FileName>("slip_sys_hard_prop_file_name")),
112 420 : _rtol(getParam<Real>("rtol")),
113 420 : _abs_tol(getParam<Real>("abs_tol")),
114 420 : _gtol(getParam<Real>("gtol")),
115 420 : _slip_incr_tol(getParam<Real>("slip_incr_tol")),
116 420 : _maxiter(getParam<unsigned int>("maxiter")),
117 420 : _maxiterg(getParam<unsigned int>("maxitergss")),
118 420 : _num_slip_sys_flowrate_props(getParam<unsigned int>("num_slip_sys_flowrate_props")),
119 420 : _tan_mod_type(getParam<MooseEnum>("tan_mod_type")),
120 420 : _intvar_read_type(getParam<MooseEnum>("intvar_read_type")),
121 420 : _num_slip_sys_props(getParam<unsigned int>("num_slip_sys_props")),
122 420 : _gen_rndm_stress_flag(getParam<bool>("gen_random_stress_flag")),
123 420 : _input_rndm_scale_var(getParam<bool>("input_random_scaling_var")),
124 420 : _rndm_scale_var(getParam<Real>("random_scaling_var")),
125 420 : _rndm_seed(getParam<unsigned int>("random_seed")),
126 420 : _max_substep_iter(getParam<unsigned int>("maximum_substep_iteration")),
127 420 : _use_line_search(getParam<bool>("use_line_search")),
128 420 : _min_lsrch_step(getParam<Real>("min_line_search_step_size")),
129 420 : _lsrch_tol(getParam<Real>("line_search_tol")),
130 420 : _lsrch_max_iter(getParam<unsigned int>("line_search_maxiter")),
131 420 : _lsrch_method(getParam<MooseEnum>("line_search_method")),
132 210 : _fp(declareProperty<RankTwoTensor>("fp")), // Plastic deformation gradient
133 420 : _fp_old(getMaterialPropertyOld<RankTwoTensor>(
134 : "fp")), // Plastic deformation gradient of previous increment
135 210 : _pk2(declareProperty<RankTwoTensor>("pk2")), // 2nd Piola-Kirchoff Stress
136 420 : _pk2_old(getMaterialPropertyOld<RankTwoTensor>(
137 : "pk2")), // 2nd Piola Kirchoff Stress of previous increment
138 210 : _lag_e(declareProperty<RankTwoTensor>("lage")), // Lagrangian strain
139 210 : _lag_e_old(
140 210 : getMaterialPropertyOld<RankTwoTensor>("lage")), // Lagrangian strain of previous increment
141 210 : _gss(declareProperty<std::vector<Real>>("gss")), // Slip system resistances
142 420 : _gss_old(getMaterialPropertyOld<std::vector<Real>>(
143 : "gss")), // Slip system resistances of previous increment
144 210 : _acc_slip(declareProperty<Real>("acc_slip")), // Accumulated slip
145 210 : _acc_slip_old(
146 210 : getMaterialPropertyOld<Real>("acc_slip")), // Accumulated slip of previous increment
147 210 : _update_rot(declareProperty<RankTwoTensor>(
148 : "update_rot")), // Rotation tensor considering material rotation and crystal orientation
149 420 : _deformation_gradient(getMaterialProperty<RankTwoTensor>("deformation_gradient")),
150 420 : _deformation_gradient_old(getMaterialPropertyOld<RankTwoTensor>("deformation_gradient")),
151 210 : _elasticity_tensor_name(_base_name + "elasticity_tensor"),
152 210 : _elasticity_tensor(getMaterialPropertyByName<RankFourTensor>(_elasticity_tensor_name)),
153 420 : _crysrot(getMaterialProperty<RankTwoTensor>("crysrot")),
154 210 : _mo(_nss * LIBMESH_DIM),
155 210 : _no(_nss * LIBMESH_DIM),
156 210 : _slip_incr(_nss),
157 210 : _tau(_nss),
158 210 : _dslipdtau(_nss),
159 210 : _s0(_nss),
160 210 : _gss_tmp(_nss),
161 210 : _gss_tmp_old(_nss),
162 840 : _dgss_dsliprate(_nss, _nss)
163 : {
164 210 : _err_tol = false;
165 :
166 210 : if (_num_slip_sys_props > 0)
167 18 : _slip_sys_props.resize(_nss * _num_slip_sys_props);
168 :
169 : _pk2_tmp.zero();
170 : _delta_dfgrd.zero();
171 :
172 210 : _first_step_iter = false;
173 210 : _last_step_iter = false;
174 : // Initialize variables in the first iteration of substepping
175 210 : _first_substep = true;
176 :
177 210 : _read_from_slip_sys_file = false;
178 210 : if (_intvar_read_type == "slip_sys_file")
179 18 : _read_from_slip_sys_file = true;
180 :
181 210 : if (_read_from_slip_sys_file && !(_num_slip_sys_props > 0))
182 0 : mooseError("Crystal Plasticity Error: Specify number of internal variable's initial values to "
183 : "be read from slip system file");
184 :
185 210 : getSlipSystems();
186 :
187 210 : RankTwoTensor::initRandom(_rndm_seed);
188 210 : }
189 :
190 : void
191 896 : FiniteStrainCrystalPlasticity::initQpStatefulProperties()
192 : {
193 896 : _stress[_qp].zero();
194 :
195 896 : _fp[_qp].setToIdentity();
196 :
197 896 : _pk2[_qp].zero();
198 896 : _acc_slip[_qp] = 0.0;
199 896 : _lag_e[_qp].zero();
200 :
201 896 : _update_rot[_qp].setToIdentity();
202 :
203 896 : initSlipSysProps(); // Initializes slip system related properties
204 896 : initAdditionalProps();
205 896 : }
206 :
207 : void
208 896 : FiniteStrainCrystalPlasticity::initSlipSysProps()
209 : {
210 896 : switch (_intvar_read_type)
211 : {
212 64 : case 0:
213 64 : assignSlipSysRes();
214 64 : break;
215 64 : case 1:
216 64 : readFileInitSlipSysRes();
217 64 : break;
218 768 : default:
219 768 : getInitSlipSysRes();
220 : }
221 :
222 896 : if (_slip_sys_flow_prop_file_name.length() != 0)
223 64 : readFileFlowRateParams();
224 : else
225 832 : getFlowRateParams();
226 :
227 896 : if (_slip_sys_hard_prop_file_name.length() != 0)
228 0 : readFileHardnessParams();
229 : else
230 896 : getHardnessParams();
231 896 : }
232 :
233 : void
234 64 : FiniteStrainCrystalPlasticity::assignSlipSysRes()
235 : {
236 64 : _gss[_qp].resize(_nss);
237 :
238 832 : for (unsigned int i = 0; i < _nss; ++i)
239 768 : _gss[_qp][i] = _slip_sys_props(i);
240 64 : }
241 :
242 : // Read initial slip system resistances from .txt file. See test.
243 : void
244 64 : FiniteStrainCrystalPlasticity::readFileInitSlipSysRes()
245 : {
246 64 : _gss[_qp].resize(_nss);
247 :
248 64 : MooseUtils::checkFileReadable(_slip_sys_res_prop_file_name);
249 :
250 64 : std::ifstream file;
251 64 : file.open(_slip_sys_res_prop_file_name.c_str());
252 :
253 832 : for (unsigned int i = 0; i < _nss; ++i)
254 768 : if (!(file >> _gss[_qp][i]))
255 0 : mooseError("Error FiniteStrainCrystalPlasticity: Premature end of slip_sys_res_prop file");
256 :
257 64 : file.close();
258 64 : }
259 :
260 : // Read initial slip system resistances from .i file
261 : void
262 768 : FiniteStrainCrystalPlasticity::getInitSlipSysRes()
263 : {
264 768 : if (_gprops.size() <= 0)
265 0 : mooseError("FiniteStrainCrystalPLasticity: Error in reading slip system resistance properties: "
266 : "Specify input in .i file or in slip_sys_res_prop_file or in slip_sys_file");
267 :
268 768 : _gss[_qp].resize(_nss, 0.0);
269 :
270 : unsigned int num_data_grp = 3; // Number of data per group e.g. start_slip_sys, end_slip_sys,
271 : // value
272 :
273 3072 : for (unsigned int i = 0; i < _gprops.size() / num_data_grp; ++i)
274 : {
275 : Real vs, ve;
276 : unsigned int is, ie;
277 :
278 2304 : vs = _gprops[i * num_data_grp];
279 2304 : ve = _gprops[i * num_data_grp + 1];
280 :
281 2304 : if (vs <= 0 || ve <= 0)
282 0 : mooseError("FiniteStrainCrystalPLasticity: Indices in gss property read must be positive "
283 : "integers: is = ",
284 : vs,
285 : " ie = ",
286 : ve);
287 :
288 2304 : if (vs != floor(vs) || ve != floor(ve))
289 0 : mooseError("FiniteStrainCrystalPLasticity: Error in reading slip system resistances: Values "
290 : "specifying start and end number of slip system groups should be integer");
291 :
292 2304 : is = static_cast<unsigned int>(vs);
293 2304 : ie = static_cast<unsigned int>(ve);
294 :
295 2304 : if (is > ie)
296 0 : mooseError("FiniteStrainCrystalPLasticity: Start index is = ",
297 : is,
298 : " should be greater than end index ie = ",
299 : ie,
300 : " in slip system resistance property read");
301 :
302 11520 : for (unsigned int j = is; j <= ie; ++j)
303 9216 : _gss[_qp][j - 1] = _gprops[i * num_data_grp + 2];
304 : }
305 :
306 9984 : for (unsigned int i = 0; i < _nss; ++i)
307 9216 : if (_gss[_qp][i] <= 0.0)
308 0 : mooseError("FiniteStrainCrystalPLasticity: Value of resistance for slip system ",
309 0 : i + 1,
310 : " non positive");
311 768 : }
312 :
313 : // Read flow rate parameters from .txt file. See test.
314 : void
315 64 : FiniteStrainCrystalPlasticity::readFileFlowRateParams()
316 : {
317 64 : _a0.resize(_nss);
318 64 : _xm.resize(_nss);
319 :
320 64 : MooseUtils::checkFileReadable(_slip_sys_flow_prop_file_name);
321 :
322 64 : std::ifstream file;
323 64 : file.open(_slip_sys_flow_prop_file_name.c_str());
324 :
325 : std::vector<Real> vec;
326 64 : vec.resize(_num_slip_sys_flowrate_props);
327 :
328 832 : for (unsigned int i = 0; i < _nss; ++i)
329 : {
330 2304 : for (unsigned int j = 0; j < _num_slip_sys_flowrate_props; ++j)
331 1536 : if (!(file >> vec[j]))
332 0 : mooseError(
333 : "Error FiniteStrainCrystalPlasticity: Premature end of slip_sys_flow_rate_param file");
334 :
335 768 : _a0(i) = vec[0];
336 768 : _xm(i) = vec[1];
337 : }
338 :
339 64 : file.close();
340 64 : }
341 :
342 : // Read flow rate parameters from .i file
343 : void
344 832 : FiniteStrainCrystalPlasticity::getFlowRateParams()
345 : {
346 832 : if (_flowprops.size() <= 0)
347 0 : mooseError("FiniteStrainCrystalPLasticity: Error in reading flow rate properties: Specify "
348 : "input in .i file or a slip_sys_flow_prop_file_name");
349 :
350 832 : _a0.resize(_nss);
351 832 : _xm.resize(_nss);
352 :
353 832 : unsigned int num_data_grp = 2 + _num_slip_sys_flowrate_props; // Number of data per group e.g.
354 : // start_slip_sys, end_slip_sys,
355 : // value1, value2, ..
356 :
357 3328 : for (unsigned int i = 0; i < _flowprops.size() / num_data_grp; ++i)
358 : {
359 : Real vs, ve;
360 : unsigned int is, ie;
361 :
362 2496 : vs = _flowprops[i * num_data_grp];
363 2496 : ve = _flowprops[i * num_data_grp + 1];
364 :
365 2496 : if (vs <= 0 || ve <= 0)
366 0 : mooseError("FiniteStrainCrystalPLasticity: Indices in flow rate parameter read must be "
367 : "positive integers: is = ",
368 : vs,
369 : " ie = ",
370 : ve);
371 :
372 2496 : if (vs != floor(vs) || ve != floor(ve))
373 0 : mooseError("FiniteStrainCrystalPLasticity: Error in reading flow props: Values specifying "
374 : "start and end number of slip system groups should be integer");
375 :
376 2496 : is = static_cast<unsigned int>(vs);
377 2496 : ie = static_cast<unsigned int>(ve);
378 :
379 2496 : if (is > ie)
380 0 : mooseError("FiniteStrainCrystalPLasticity: Start index is = ",
381 : is,
382 : " should be greater than end index ie = ",
383 : ie,
384 : " in flow rate parameter read");
385 :
386 12480 : for (unsigned int j = is; j <= ie; ++j)
387 : {
388 9984 : _a0(j - 1) = _flowprops[i * num_data_grp + 2];
389 9984 : _xm(j - 1) = _flowprops[i * num_data_grp + 3];
390 : }
391 : }
392 :
393 10816 : for (unsigned int i = 0; i < _nss; ++i)
394 : {
395 9984 : if (!(_a0(i) > 0.0 && _xm(i) > 0.0))
396 : {
397 0 : mooseWarning(
398 : "FiniteStrainCrystalPlasticity: Non-positive flow rate parameters ", _a0(i), ",", _xm(i));
399 0 : break;
400 : }
401 : }
402 832 : }
403 :
404 : // Read hardness parameters from .txt file
405 : void
406 0 : FiniteStrainCrystalPlasticity::readFileHardnessParams()
407 : {
408 0 : }
409 :
410 : // Read hardness parameters from .i file
411 : void
412 896 : FiniteStrainCrystalPlasticity::getHardnessParams()
413 : {
414 896 : if (_hprops.size() <= 0)
415 0 : mooseError("FiniteStrainCrystalPLasticity: Error in reading hardness properties: Specify input "
416 : "in .i file or a slip_sys_hard_prop_file_name");
417 :
418 896 : _r = _hprops[0];
419 896 : _h0 = _hprops[1];
420 896 : _tau_init = _hprops[2];
421 896 : _tau_sat = _hprops[3];
422 896 : }
423 :
424 : // Read slip systems from file
425 : void
426 210 : FiniteStrainCrystalPlasticity::getSlipSystems()
427 : {
428 : Real vec[LIBMESH_DIM];
429 210 : std::ifstream fileslipsys;
430 :
431 210 : MooseUtils::checkFileReadable(_slip_sys_file_name);
432 :
433 210 : fileslipsys.open(_slip_sys_file_name.c_str());
434 :
435 2730 : for (unsigned int i = 0; i < _nss; ++i)
436 : {
437 : // Read the slip normal
438 10080 : for (const auto j : make_range(Moose::dim))
439 7560 : if (!(fileslipsys >> vec[j]))
440 0 : mooseError("Crystal Plasticity Error: Premature end of file reading slip system file \n");
441 :
442 : // Normalize the vectors
443 : Real mag;
444 2520 : mag = Utility::pow<2>(vec[0]) + Utility::pow<2>(vec[1]) + Utility::pow<2>(vec[2]);
445 2520 : mag = std::sqrt(mag);
446 :
447 10080 : for (unsigned j = 0; j < LIBMESH_DIM; ++j)
448 7560 : _no(i * LIBMESH_DIM + j) = vec[j] / mag;
449 :
450 : // Read the slip direction
451 10080 : for (const auto j : make_range(Moose::dim))
452 7560 : if (!(fileslipsys >> vec[j]))
453 0 : mooseError("Crystal Plasticity Error: Premature end of file reading slip system file \n");
454 :
455 : // Normalize the vectors
456 2520 : mag = Utility::pow<2>(vec[0]) + Utility::pow<2>(vec[1]) + Utility::pow<2>(vec[2]);
457 2520 : mag = std::sqrt(mag);
458 :
459 10080 : for (const auto j : make_range(Moose::dim))
460 7560 : _mo(i * LIBMESH_DIM + j) = vec[j] / mag;
461 :
462 : mag = 0.0;
463 10080 : for (const auto j : make_range(Moose::dim))
464 7560 : mag += _mo(i * LIBMESH_DIM + j) * _no(i * LIBMESH_DIM + j);
465 :
466 2520 : if (std::abs(mag) > 1e-8)
467 0 : mooseError(
468 : "Crystal Plasicity Error: Slip direction and normal not orthonormal, System number = ",
469 : i,
470 : "\n");
471 :
472 2520 : if (_read_from_slip_sys_file)
473 432 : for (unsigned int j = 0; j < _num_slip_sys_props; ++j)
474 216 : if (!(fileslipsys >> _slip_sys_props(i * _num_slip_sys_props + j)))
475 0 : mooseError("Crystal Plasticity Error: Premature end of file reading slip system file - "
476 : "check in slip system file read input options/values\n");
477 : }
478 :
479 210 : fileslipsys.close();
480 210 : }
481 :
482 : // Initialize addtional stateful material properties
483 : void
484 896 : FiniteStrainCrystalPlasticity::initAdditionalProps()
485 : {
486 896 : }
487 :
488 : /**
489 : * Solves stress residual equation using NR.
490 : * Updates slip system resistances iteratively.
491 : */
492 : void
493 67136 : FiniteStrainCrystalPlasticity::computeQpStress()
494 : {
495 : unsigned int substep_iter = 1; // Depth of substepping; Limited to maximum substep iteration
496 : unsigned int num_substep = 1; // Calculated from substep_iter as 2^substep_iter
497 67136 : Real dt_original = _dt; // Stores original _dt; Reset at the end of solve
498 67136 : _first_substep = true; // Initialize variables at substep_iter = 1
499 :
500 67136 : if (_max_substep_iter > 1)
501 : {
502 8912 : _dfgrd_tmp_old = _deformation_gradient_old[_qp];
503 8912 : if (_dfgrd_tmp_old.det() == 0)
504 0 : _dfgrd_tmp_old.addIa(1.0);
505 :
506 8912 : _delta_dfgrd = _deformation_gradient[_qp] - _dfgrd_tmp_old;
507 8912 : _err_tol = true; // Indicator to continue substepping
508 : }
509 :
510 : // Substepping loop
511 96224 : while (_err_tol && _max_substep_iter > 1)
512 : {
513 29088 : _dt = dt_original / num_substep;
514 :
515 92768 : for (unsigned int istep = 0; istep < num_substep; ++istep)
516 : {
517 83856 : _first_step_iter = false;
518 83856 : if (istep == 0)
519 29088 : _first_step_iter = true;
520 :
521 83856 : _last_step_iter = false;
522 83856 : if (istep == num_substep - 1)
523 16064 : _last_step_iter = true;
524 :
525 83856 : _dfgrd_scale_factor = (static_cast<Real>(istep) + 1) / num_substep;
526 :
527 83856 : preSolveQp();
528 83856 : solveQp();
529 :
530 83856 : if (_err_tol)
531 : {
532 20176 : substep_iter++;
533 20176 : num_substep *= 2;
534 20176 : break;
535 : }
536 : }
537 :
538 29088 : _first_substep = false; // Prevents reinitialization
539 29088 : _dt = dt_original; // Resets dt
540 :
541 : #ifdef DEBUG
542 : if (substep_iter > _max_substep_iter)
543 : mooseWarning("FiniteStrainCrystalPlasticity: Failure with substepping");
544 : #endif
545 :
546 29088 : if (!_err_tol || substep_iter > _max_substep_iter)
547 8912 : postSolveQp(); // Evaluate variables after successful solve or indicate failure
548 : }
549 :
550 : // No substepping
551 67136 : if (_max_substep_iter == 1)
552 : {
553 58224 : preSolveQp();
554 58224 : solveQp();
555 58224 : postSolveQp();
556 : }
557 67136 : }
558 :
559 : void
560 142080 : FiniteStrainCrystalPlasticity::preSolveQp()
561 : {
562 : // Initialize variable
563 142080 : if (_first_substep)
564 : {
565 67136 : _Jacobian_mult[_qp].zero(); // Initializes jacobian for preconditioner
566 67136 : calc_schmid_tensor();
567 : }
568 :
569 142080 : if (_max_substep_iter == 1)
570 58224 : _dfgrd_tmp = _deformation_gradient[_qp]; // Without substepping
571 : else
572 83856 : _dfgrd_tmp = _dfgrd_scale_factor * _delta_dfgrd + _dfgrd_tmp_old;
573 :
574 142080 : _err_tol = false;
575 142080 : }
576 :
577 : void
578 142080 : FiniteStrainCrystalPlasticity::solveQp()
579 : {
580 142080 : preSolveStatevar();
581 142080 : solveStatevar();
582 142080 : if (_err_tol)
583 : return;
584 112784 : postSolveStatevar();
585 : }
586 :
587 : void
588 67136 : FiniteStrainCrystalPlasticity::postSolveQp()
589 : {
590 67136 : if (_err_tol)
591 : {
592 9120 : _err_tol = false;
593 9120 : if (_gen_rndm_stress_flag)
594 : {
595 9120 : if (!_input_rndm_scale_var)
596 9120 : _rndm_scale_var = _elasticity_tensor[_qp](0, 0, 0, 0);
597 :
598 9120 : _stress[_qp] = RankTwoTensor::genRandomSymmTensor(_rndm_scale_var, 1.0);
599 : }
600 : else
601 0 : mooseError("FiniteStrainCrystalPlasticity: Constitutive failure");
602 : }
603 : else
604 : {
605 58016 : _stress[_qp] = _fe * _pk2[_qp] * _fe.transpose() / _fe.det();
606 :
607 58016 : _Jacobian_mult[_qp] += calcTangentModuli(); // Calculate jacobian for preconditioner
608 :
609 58016 : RankTwoTensor iden(RankTwoTensor::initIdentity);
610 :
611 58016 : _lag_e[_qp] = _deformation_gradient[_qp].transpose() * _deformation_gradient[_qp] - iden;
612 58016 : _lag_e[_qp] = _lag_e[_qp] * 0.5;
613 :
614 58016 : RankTwoTensor rot;
615 58016 : rot = get_current_rotation(_deformation_gradient[_qp]); // Calculate material rotation
616 58016 : _update_rot[_qp] = rot * _crysrot[_qp];
617 : }
618 67136 : }
619 :
620 : void
621 142080 : FiniteStrainCrystalPlasticity::preSolveStatevar()
622 : {
623 142080 : if (_max_substep_iter == 1) // No substepping
624 : {
625 58224 : _gss_tmp = _gss_old[_qp];
626 58224 : _accslip_tmp_old = _acc_slip_old[_qp];
627 : }
628 : else
629 : {
630 83856 : if (_first_step_iter)
631 : {
632 29088 : _gss_tmp = _gss_tmp_old = _gss_old[_qp];
633 29088 : _accslip_tmp_old = _acc_slip_old[_qp];
634 : }
635 : else
636 54768 : _gss_tmp = _gss_tmp_old;
637 : }
638 142080 : }
639 :
640 : void
641 58480 : FiniteStrainCrystalPlasticity::solveStatevar()
642 : {
643 : Real gmax, gdiff;
644 : unsigned int iterg;
645 58480 : std::vector<Real> gss_prev(_nss);
646 :
647 58480 : gmax = 1.1 * _gtol;
648 : iterg = 0;
649 :
650 121104 : while (gmax > _gtol && iterg < _maxiterg) // Check for slip system resistance update tolerance
651 : {
652 72480 : preSolveStress();
653 72480 : solveStress();
654 72480 : if (_err_tol)
655 : return;
656 62624 : postSolveStress();
657 :
658 62624 : gss_prev = _gss_tmp;
659 :
660 62624 : update_slip_system_resistance(); // Update slip system resistance
661 :
662 : gmax = 0.0;
663 814112 : for (unsigned i = 0; i < _nss; ++i)
664 : {
665 751488 : gdiff = std::abs(gss_prev[i] - _gss_tmp[i]); // Calculate increment size
666 :
667 751488 : if (gdiff > gmax)
668 : gmax = gdiff;
669 : }
670 62624 : iterg++;
671 : }
672 :
673 48624 : if (iterg == _maxiterg)
674 : {
675 : #ifdef DEBUG
676 : mooseWarning("FiniteStrainCrystalPLasticity: Hardness Integration error gmax", gmax, "\n");
677 : #endif
678 0 : _err_tol = true;
679 : }
680 : }
681 :
682 : void
683 112784 : FiniteStrainCrystalPlasticity::postSolveStatevar()
684 : {
685 112784 : if (_max_substep_iter == 1) // No substepping
686 : {
687 49104 : _gss[_qp] = _gss_tmp;
688 49104 : _acc_slip[_qp] = _accslip_tmp;
689 : }
690 : else
691 : {
692 63680 : if (_last_step_iter)
693 : {
694 8912 : _gss[_qp] = _gss_tmp;
695 8912 : _acc_slip[_qp] = _accslip_tmp;
696 : }
697 : else
698 : {
699 54768 : _gss_tmp_old = _gss_tmp;
700 54768 : _accslip_tmp_old = _accslip_tmp;
701 : }
702 : }
703 112784 : }
704 :
705 : void
706 156080 : FiniteStrainCrystalPlasticity::preSolveStress()
707 : {
708 156080 : if (_max_substep_iter == 1) // No substepping
709 : {
710 69856 : _pk2_tmp = _pk2_old[_qp];
711 69856 : _fp_old_inv = _fp_old[_qp].inverse();
712 69856 : _fp_inv = _fp_old_inv;
713 69856 : _fp_prev_inv = _fp_inv;
714 : }
715 : else
716 : {
717 86224 : if (_first_step_iter)
718 : {
719 30720 : _pk2_tmp = _pk2_tmp_old = _pk2_old[_qp];
720 30720 : _fp_old_inv = _fp_old[_qp].inverse();
721 : }
722 : else
723 55504 : _pk2_tmp = _pk2_tmp_old;
724 :
725 86224 : _fp_inv = _fp_old_inv;
726 86224 : _fp_prev_inv = _fp_inv;
727 : }
728 156080 : }
729 :
730 : void
731 72480 : FiniteStrainCrystalPlasticity::solveStress()
732 : {
733 : unsigned int iter = 0;
734 72480 : RankTwoTensor resid, dpk2;
735 72480 : RankFourTensor jac;
736 : Real rnorm, rnorm0, rnorm_prev;
737 :
738 72480 : calc_resid_jacob(resid, jac); // Calculate stress residual
739 72480 : if (_err_tol)
740 : {
741 : #ifdef DEBUG
742 : mooseWarning(
743 : "FiniteStrainCrystalPLasticity: Slip increment exceeds tolerance - Element number ",
744 : _current_elem->id(),
745 : " Gauss point = ",
746 : _qp);
747 : #endif
748 : return;
749 : }
750 :
751 72480 : rnorm = resid.L2norm();
752 : rnorm0 = rnorm;
753 :
754 377472 : while (rnorm > _rtol * rnorm0 && rnorm0 > _abs_tol &&
755 314848 : iter < _maxiter) // Check for stress residual tolerance
756 : {
757 314848 : dpk2 = -jac.invSymm() * resid; // Calculate stress increment
758 314848 : _pk2_tmp = _pk2_tmp + dpk2; // Update stress
759 314848 : calc_resid_jacob(resid, jac);
760 314848 : internalVariableUpdateNRiteration(); // update _fp_prev_inv
761 :
762 314848 : if (_err_tol)
763 : {
764 : #ifdef DEBUG
765 : mooseWarning(
766 : "FiniteStrainCrystalPLasticity: Slip increment exceeds tolerance - Element number ",
767 : _current_elem->id(),
768 : " Gauss point = ",
769 : _qp);
770 : #endif
771 : return;
772 : }
773 :
774 : rnorm_prev = rnorm;
775 304992 : rnorm = resid.L2norm();
776 :
777 304992 : if (_use_line_search && rnorm > rnorm_prev && !line_search_update(rnorm_prev, dpk2))
778 : {
779 : #ifdef DEBUG
780 : mooseWarning("FiniteStrainCrystalPLasticity: Failed with line search");
781 : #endif
782 0 : _err_tol = true;
783 0 : return;
784 : }
785 :
786 304992 : if (_use_line_search)
787 13312 : rnorm = resid.L2norm();
788 :
789 304992 : iter++;
790 : }
791 :
792 62624 : if (iter >= _maxiter)
793 : {
794 : #ifdef DEBUG
795 : mooseWarning("FiniteStrainCrystalPLasticity: Stress Integration error rmax = ", rnorm);
796 : #endif
797 0 : _err_tol = true;
798 : }
799 : }
800 :
801 : void
802 126784 : FiniteStrainCrystalPlasticity::postSolveStress()
803 : {
804 126784 : if (_max_substep_iter == 1) // No substepping
805 : {
806 60736 : _fp[_qp] = _fp_inv.inverse();
807 60736 : _pk2[_qp] = _pk2_tmp;
808 : }
809 : else
810 : {
811 66048 : if (_last_step_iter)
812 : {
813 10544 : _fp[_qp] = _fp_inv.inverse();
814 10544 : _pk2[_qp] = _pk2_tmp;
815 : }
816 : else
817 : {
818 55504 : _fp_old_inv = _fp_inv;
819 55504 : _pk2_tmp_old = _pk2_tmp;
820 : }
821 : }
822 126784 : }
823 :
824 : // Update slip system resistance. Overide to incorporate new slip system resistance laws
825 : void
826 694672 : FiniteStrainCrystalPlasticity::update_slip_system_resistance()
827 : {
828 694672 : updateGss();
829 694672 : }
830 :
831 : /**
832 : * Old function to update slip system resistances.
833 : * Kept to avoid code break at computeQpstress
834 : */
835 : void
836 694672 : FiniteStrainCrystalPlasticity::updateGss()
837 : {
838 694672 : DenseVector<Real> hb(_nss);
839 : Real qab;
840 :
841 694672 : Real a = _hprops[4]; // Kalidindi
842 :
843 694672 : _accslip_tmp = _accslip_tmp_old;
844 9030736 : for (unsigned int i = 0; i < _nss; ++i)
845 8336064 : _accslip_tmp += std::abs(_slip_incr(i));
846 :
847 : // Real val = std::cosh(_h0 * _accslip_tmp / (_tau_sat - _tau_init)); // Karthik
848 : // val = _h0 * std::pow(1.0/val,2.0); // Kalidindi
849 :
850 9030736 : for (unsigned int i = 0; i < _nss; ++i)
851 : // hb(i)=val;
852 8336064 : hb(i) = _h0 * std::pow(std::abs(1.0 - _gss_tmp[i] / _tau_sat), a) *
853 8336064 : std::copysign(1.0, 1.0 - _gss_tmp[i] / _tau_sat);
854 :
855 9030736 : for (unsigned int i = 0; i < _nss; ++i)
856 : {
857 8336064 : if (_max_substep_iter == 1) // No substepping
858 950016 : _gss_tmp[i] = _gss_old[_qp][i];
859 : else
860 7386048 : _gss_tmp[i] = _gss_tmp_old[i];
861 :
862 108368832 : for (unsigned int j = 0; j < _nss; ++j)
863 : {
864 : unsigned int iplane, jplane;
865 100032768 : iplane = i / 3;
866 100032768 : jplane = j / 3;
867 :
868 100032768 : if (iplane == jplane) // Kalidindi
869 : qab = 1.0;
870 : else
871 75024576 : qab = _r;
872 :
873 100032768 : _gss_tmp[i] += qab * hb(j) * std::abs(_slip_incr(j));
874 100032768 : _dgss_dsliprate(i, j) = qab * hb(j) * std::copysign(1.0, _slip_incr(j)) * _dt;
875 : }
876 : }
877 694672 : }
878 :
879 : // Calculates stress residual equation and jacobian
880 : void
881 387328 : FiniteStrainCrystalPlasticity::calc_resid_jacob(RankTwoTensor & resid, RankFourTensor & jac)
882 : {
883 387328 : calcResidual(resid);
884 387328 : if (_err_tol)
885 : return;
886 377472 : calcJacobian(jac);
887 : }
888 :
889 : void
890 387328 : FiniteStrainCrystalPlasticity::calcResidual(RankTwoTensor & resid)
891 : {
892 387328 : RankTwoTensor iden(RankTwoTensor::initIdentity), ce, ee, ce_pk2, eqv_slip_incr, pk2_new;
893 :
894 387328 : _fe = _dfgrd_tmp * _fp_prev_inv; // _fp_inv ==> _fp_prev_inv
895 :
896 387328 : ce = _fe.transpose() * _fe;
897 387328 : ce_pk2 = ce * _pk2_tmp;
898 387328 : ce_pk2 = ce_pk2 / _fe.det();
899 :
900 : // Calculate Schmid tensor and resolved shear stresses
901 5035264 : for (unsigned int i = 0; i < _nss; ++i)
902 4647936 : _tau(i) = ce_pk2.doubleContraction(_s0[i]);
903 :
904 387328 : getSlipIncrements(); // Calculate dslip,dslipdtau
905 :
906 387328 : if (_err_tol)
907 : return;
908 :
909 : eqv_slip_incr.zero();
910 4907136 : for (unsigned int i = 0; i < _nss; ++i)
911 4529664 : eqv_slip_incr += _s0[i] * _slip_incr(i);
912 :
913 377472 : eqv_slip_incr = iden - eqv_slip_incr;
914 377472 : _fp_inv = _fp_old_inv * eqv_slip_incr;
915 377472 : _fe = _dfgrd_tmp * _fp_inv;
916 :
917 377472 : ce = _fe.transpose() * _fe;
918 377472 : ee = ce - iden;
919 377472 : ee *= 0.5;
920 :
921 377472 : pk2_new = _elasticity_tensor[_qp] * ee;
922 :
923 377472 : resid = _pk2_tmp - pk2_new;
924 : }
925 :
926 : void
927 377472 : FiniteStrainCrystalPlasticity::calcJacobian(RankFourTensor & jac)
928 : {
929 377472 : RankFourTensor dfedfpinv, deedfe, dfpinvdpk2;
930 :
931 377472 : std::vector<RankTwoTensor> dtaudpk2(_nss), dfpinvdslip(_nss);
932 :
933 4907136 : for (unsigned int i = 0; i < _nss; ++i)
934 : {
935 4529664 : dtaudpk2[i] = _s0[i];
936 4529664 : dfpinvdslip[i] = -_fp_old_inv * _s0[i];
937 : }
938 :
939 1509888 : for (const auto i : make_range(Moose::dim))
940 4529664 : for (const auto j : make_range(Moose::dim))
941 13588992 : for (const auto k : make_range(Moose::dim))
942 10191744 : dfedfpinv(i, j, k, j) = _dfgrd_tmp(i, k);
943 :
944 1509888 : for (const auto i : make_range(Moose::dim))
945 4529664 : for (const auto j : make_range(Moose::dim))
946 13588992 : for (const auto k : make_range(Moose::dim))
947 : {
948 10191744 : deedfe(i, j, k, i) = deedfe(i, j, k, i) + _fe(k, j) * 0.5;
949 10191744 : deedfe(i, j, k, j) = deedfe(i, j, k, j) + _fe(k, i) * 0.5;
950 : }
951 :
952 4907136 : for (unsigned int i = 0; i < _nss; ++i)
953 4529664 : dfpinvdpk2 += (dfpinvdslip[i] * _dslipdtau(i)).outerProduct(dtaudpk2[i]);
954 :
955 : jac =
956 754944 : RankFourTensor::IdentityFour() - (_elasticity_tensor[_qp] * deedfe * dfedfpinv * dfpinvdpk2);
957 377472 : }
958 :
959 : // Calculate slip increment,dslipdtau. Override to modify.
960 : void
961 1019376 : FiniteStrainCrystalPlasticity::getSlipIncrements()
962 : {
963 12900336 : for (unsigned int i = 0; i < _nss; ++i)
964 : {
965 11910256 : _slip_incr(i) = _a0(i) * std::pow(std::abs(_tau(i) / _gss_tmp[i]), 1.0 / _xm(i)) *
966 11910256 : std::copysign(1.0, _tau(i)) * _dt;
967 11910256 : if (std::abs(_slip_incr(i)) > _slip_incr_tol)
968 : {
969 29296 : _err_tol = true;
970 : #ifdef DEBUG
971 : mooseWarning("Maximum allowable slip increment exceeded ", std::abs(_slip_incr(i)));
972 : #endif
973 29296 : return;
974 : }
975 : }
976 :
977 12871040 : for (unsigned int i = 0; i < _nss; ++i)
978 11880960 : _dslipdtau(i) = _a0(i) / _xm(i) *
979 11880960 : std::pow(std::abs(_tau(i) / _gss_tmp[i]), 1.0 / _xm(i) - 1.0) / _gss_tmp[i] *
980 11880960 : _dt;
981 : }
982 :
983 : // Calls getMatRot to perform RU factorization of a tensor.
984 : RankTwoTensor
985 58016 : FiniteStrainCrystalPlasticity::get_current_rotation(const RankTwoTensor & a)
986 : {
987 58016 : return getMatRot(a);
988 : }
989 :
990 : // Performs RU factorization of a tensor
991 : RankTwoTensor
992 58016 : FiniteStrainCrystalPlasticity::getMatRot(const RankTwoTensor & a)
993 : {
994 58016 : RankTwoTensor rot;
995 58016 : RankTwoTensor c, diag, evec;
996 : PetscScalar cmat[LIBMESH_DIM][LIBMESH_DIM], work[10];
997 : PetscReal w[LIBMESH_DIM];
998 58016 : PetscBLASInt nd = LIBMESH_DIM, lwork = 10, info;
999 :
1000 58016 : c = a.transpose() * a;
1001 :
1002 232064 : for (const auto i : make_range(Moose::dim))
1003 696192 : for (const auto j : make_range(Moose::dim))
1004 522144 : cmat[i][j] = c(i, j);
1005 :
1006 58016 : LAPACKsyev_("V", "U", &nd, &cmat[0][0], &nd, w, work, &lwork, &info);
1007 :
1008 58016 : if (info != 0)
1009 0 : mooseError("FiniteStrainCrystalPLasticity: DSYEV function call in getMatRot function failed");
1010 :
1011 : diag.zero();
1012 :
1013 232064 : for (const auto i : make_range(Moose::dim))
1014 174048 : diag(i, i) = std::sqrt(w[i]);
1015 :
1016 232064 : for (const auto i : make_range(Moose::dim))
1017 696192 : for (const auto j : make_range(Moose::dim))
1018 522144 : evec(i, j) = cmat[i][j];
1019 :
1020 58016 : rot = a * ((evec.transpose() * diag * evec).inverse());
1021 :
1022 58016 : return rot;
1023 : }
1024 :
1025 : // Calculates tangent moduli which is used for global solve
1026 : void
1027 0 : FiniteStrainCrystalPlasticity::computeQpElasticityTensor()
1028 : {
1029 0 : }
1030 :
1031 : RankFourTensor
1032 58016 : FiniteStrainCrystalPlasticity::calcTangentModuli()
1033 : {
1034 58016 : RankFourTensor tan_mod;
1035 :
1036 58016 : switch (_tan_mod_type)
1037 : {
1038 53392 : case 0:
1039 53392 : tan_mod = elastoPlasticTangentModuli();
1040 53392 : break;
1041 4624 : default:
1042 4624 : tan_mod = elasticTangentModuli();
1043 : }
1044 :
1045 58016 : return tan_mod;
1046 : }
1047 :
1048 : void
1049 67136 : FiniteStrainCrystalPlasticity::calc_schmid_tensor()
1050 : {
1051 67136 : DenseVector<Real> mo(LIBMESH_DIM * _nss), no(LIBMESH_DIM * _nss);
1052 :
1053 : // Update slip direction and normal with crystal orientation
1054 872768 : for (unsigned int i = 0; i < _nss; ++i)
1055 : {
1056 3222528 : for (const auto j : make_range(Moose::dim))
1057 : {
1058 2416896 : mo(i * LIBMESH_DIM + j) = 0.0;
1059 9667584 : for (const auto k : make_range(Moose::dim))
1060 7250688 : mo(i * LIBMESH_DIM + j) =
1061 7250688 : mo(i * LIBMESH_DIM + j) + _crysrot[_qp](j, k) * _mo(i * LIBMESH_DIM + k);
1062 : }
1063 :
1064 3222528 : for (const auto j : make_range(Moose::dim))
1065 : {
1066 2416896 : no(i * LIBMESH_DIM + j) = 0.0;
1067 9667584 : for (const auto k : make_range(Moose::dim))
1068 7250688 : no(i * LIBMESH_DIM + j) =
1069 7250688 : no(i * LIBMESH_DIM + j) + _crysrot[_qp](j, k) * _no(i * LIBMESH_DIM + k);
1070 : }
1071 : }
1072 :
1073 : // Calculate Schmid tensor and resolved shear stresses
1074 872768 : for (unsigned int i = 0; i < _nss; ++i)
1075 3222528 : for (const auto j : make_range(Moose::dim))
1076 9667584 : for (const auto k : make_range(Moose::dim))
1077 7250688 : _s0[i](j, k) = mo(i * LIBMESH_DIM + j) * no(i * LIBMESH_DIM + k);
1078 67136 : }
1079 :
1080 : RankFourTensor
1081 53392 : FiniteStrainCrystalPlasticity::elastoPlasticTangentModuli()
1082 : {
1083 53392 : RankFourTensor tan_mod;
1084 53392 : RankTwoTensor pk2fet, fepk2;
1085 53392 : RankFourTensor deedfe, dsigdpk2dfe;
1086 :
1087 : // Fill in the matrix stiffness material property
1088 :
1089 213568 : for (const auto i : make_range(Moose::dim))
1090 640704 : for (const auto j : make_range(Moose::dim))
1091 1922112 : for (const auto k : make_range(Moose::dim))
1092 : {
1093 1441584 : deedfe(i, j, k, i) = deedfe(i, j, k, i) + _fe(k, j) * 0.5;
1094 1441584 : deedfe(i, j, k, j) = deedfe(i, j, k, j) + _fe(k, i) * 0.5;
1095 : }
1096 :
1097 : usingTensorIndices(i_, j_, k_, l_);
1098 53392 : dsigdpk2dfe = _fe.times<i_, k_, j_, l_>(_fe) * _elasticity_tensor[_qp] * deedfe;
1099 :
1100 53392 : pk2fet = _pk2_tmp * _fe.transpose();
1101 53392 : fepk2 = _fe * _pk2_tmp;
1102 :
1103 213568 : for (const auto i : make_range(Moose::dim))
1104 640704 : for (const auto j : make_range(Moose::dim))
1105 1922112 : for (const auto l : make_range(Moose::dim))
1106 : {
1107 1441584 : tan_mod(i, j, i, l) = tan_mod(i, j, i, l) + pk2fet(l, j);
1108 1441584 : tan_mod(i, j, j, l) = tan_mod(i, j, j, l) + fepk2(i, l);
1109 : }
1110 :
1111 53392 : tan_mod += dsigdpk2dfe;
1112 :
1113 53392 : Real je = _fe.det();
1114 53392 : if (je > 0.0)
1115 53392 : tan_mod /= je;
1116 :
1117 53392 : return tan_mod;
1118 : }
1119 :
1120 : RankFourTensor
1121 4624 : FiniteStrainCrystalPlasticity::elasticTangentModuli()
1122 : {
1123 4624 : return _elasticity_tensor[_qp]; // update jacobian_mult
1124 : }
1125 :
1126 : bool
1127 0 : FiniteStrainCrystalPlasticity::line_search_update(const Real rnorm_prev, const RankTwoTensor dpk2)
1128 : {
1129 0 : if (_lsrch_method == "CUT_HALF")
1130 : {
1131 : Real rnorm;
1132 0 : RankTwoTensor resid;
1133 0 : Real step = 1.0;
1134 :
1135 : do
1136 : {
1137 0 : _pk2_tmp = _pk2_tmp - step * dpk2;
1138 0 : step /= 2.0;
1139 0 : _pk2_tmp = _pk2_tmp + step * dpk2;
1140 :
1141 0 : calcResidual(resid);
1142 0 : rnorm = resid.L2norm();
1143 0 : } while (rnorm > rnorm_prev && step > _min_lsrch_step);
1144 :
1145 0 : if (rnorm > rnorm_prev && step <= _min_lsrch_step)
1146 : return false;
1147 :
1148 0 : return true;
1149 : }
1150 0 : else if (_lsrch_method == "BISECTION")
1151 : {
1152 : unsigned int count = 0;
1153 : Real step_a = 0.0;
1154 : Real step_b = 1.0;
1155 0 : Real step = 1.0;
1156 : Real s_m = 1000.0;
1157 : Real rnorm = 1000.0;
1158 :
1159 0 : RankTwoTensor resid;
1160 0 : calcResidual(resid);
1161 0 : Real s_b = resid.doubleContraction(dpk2);
1162 0 : Real rnorm1 = resid.L2norm();
1163 0 : _pk2_tmp = _pk2_tmp - dpk2;
1164 0 : calcResidual(resid);
1165 0 : Real s_a = resid.doubleContraction(dpk2);
1166 0 : Real rnorm0 = resid.L2norm();
1167 0 : _pk2_tmp = _pk2_tmp + dpk2;
1168 :
1169 0 : if ((rnorm1 / rnorm0) < _lsrch_tol || s_a * s_b > 0)
1170 : {
1171 0 : calcResidual(resid);
1172 0 : return true;
1173 : }
1174 :
1175 0 : while ((rnorm / rnorm0) > _lsrch_tol && count < _lsrch_max_iter)
1176 : {
1177 0 : _pk2_tmp = _pk2_tmp - step * dpk2;
1178 0 : step = 0.5 * (step_b + step_a);
1179 0 : _pk2_tmp = _pk2_tmp + step * dpk2;
1180 0 : calcResidual(resid);
1181 0 : s_m = resid.doubleContraction(dpk2);
1182 0 : rnorm = resid.L2norm();
1183 :
1184 0 : if (s_m * s_a < 0.0)
1185 : {
1186 : step_b = step;
1187 : s_b = s_m;
1188 : }
1189 0 : if (s_m * s_b < 0.0)
1190 : {
1191 : step_a = step;
1192 : s_a = s_m;
1193 : }
1194 0 : count++;
1195 : }
1196 :
1197 0 : if ((rnorm / rnorm0) < _lsrch_tol && count < _lsrch_max_iter)
1198 : return true;
1199 :
1200 : return false;
1201 : }
1202 : else
1203 : {
1204 0 : mooseError("Line search meothod is not provided.");
1205 : return false;
1206 : }
1207 : }
1208 :
1209 : void
1210 314848 : FiniteStrainCrystalPlasticity::internalVariableUpdateNRiteration()
1211 : {
1212 314848 : _fp_prev_inv = _fp_inv; // update _fp_prev_inv
1213 314848 : }
|
