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 328 : FiniteStrainCrystalPlasticity::validParams()
23 : {
24 328 : InputParameters params = ComputeStressBase::validParams();
25 328 : 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 656 : params.addRequiredParam<int>("nss", "Number of slip systems");
29 656 : params.addParam<std::vector<Real>>("gprops", {}, "Initial values of slip system resistances");
30 656 : params.addParam<std::vector<Real>>("hprops", {}, "Hardening properties");
31 656 : params.addParam<std::vector<Real>>("flowprops", {}, "Parameters used in slip rate equations");
32 656 : params.addRequiredParam<FileName>("slip_sys_file_name",
33 : "Name of the file containing the slip system");
34 656 : 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 656 : 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 656 : params.addParam<FileName>(
43 : "slip_sys_hard_prop_file_name",
44 : "",
45 : "Name of the file containing the values of hardness evolution parameters");
46 656 : params.addParam<Real>("rtol", 1e-6, "Constitutive stress residue relative tolerance");
47 656 : params.addParam<Real>("abs_tol", 1e-6, "Constitutive stress residue absolute tolerance");
48 656 : params.addParam<Real>("gtol", 1e2, "Constitutive slip system resistance residual tolerance");
49 656 : params.addParam<Real>("slip_incr_tol", 2e-2, "Maximum allowable slip in an increment");
50 656 : params.addParam<unsigned int>("maxiter", 100, "Maximum number of iterations for stress update");
51 656 : params.addParam<unsigned int>(
52 656 : "maxitergss", 100, "Maximum number of iterations for slip system resistance update");
53 656 : params.addParam<unsigned int>(
54 : "num_slip_sys_flowrate_props",
55 656 : 2,
56 : "Number of flow rate properties for a slip system"); // Used for reading flow rate parameters
57 656 : params.addParam<UserObjectName>("read_prop_user_object",
58 : "The ElementReadPropertyFile "
59 : "GeneralUserObject to read element "
60 : "specific property values from file");
61 656 : MooseEnum tan_mod_options("exact none", "none"); // Type of read
62 656 : params.addParam<MooseEnum>("tan_mod_type",
63 : tan_mod_options,
64 : "Type of tangent moduli for preconditioner: default elastic");
65 656 : MooseEnum intvar_read_options("slip_sys_file slip_sys_res_file none", "none");
66 656 : 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 656 : params.addParam<unsigned int>("num_slip_sys_props",
71 656 : 0,
72 : "Number of slip system specific properties provided in the file "
73 : "containing slip system normals and directions");
74 656 : params.addParam<bool>(
75 : "gen_random_stress_flag",
76 656 : false,
77 : "Flag to generate random stress to perform time cutback on constitutive failure");
78 656 : params.addParam<bool>("input_random_scaling_var",
79 656 : false,
80 : "Flag to input scaling variable: _Cijkl(0,0,0,0) when false");
81 656 : params.addParam<Real>("random_scaling_var",
82 656 : 1e9,
83 : "Random scaling variable: Large value can cause non-positive definiteness");
84 656 : params.addParam<unsigned int>(
85 : "random_seed",
86 656 : 2000,
87 : "Random integer used to generate random stress when constitutive failure occurs");
88 656 : params.addParam<unsigned int>(
89 656 : "maximum_substep_iteration", 1, "Maximum number of substep iteration");
90 656 : params.addParam<bool>("use_line_search", false, "Use line search in constitutive update");
91 656 : params.addParam<Real>("min_line_search_step_size", 0.01, "Minimum line search step size");
92 656 : params.addParam<Real>("line_search_tol", 0.5, "Line search bisection method tolerance");
93 656 : params.addParam<unsigned int>(
94 656 : "line_search_maxiter", 20, "Line search bisection method maximum number of iteration");
95 656 : MooseEnum line_search_method("CUT_HALF BISECTION", "CUT_HALF");
96 656 : params.addParam<MooseEnum>(
97 : "line_search_method", line_search_method, "The method used in line search");
98 :
99 328 : return params;
100 328 : }
101 :
102 246 : FiniteStrainCrystalPlasticity::FiniteStrainCrystalPlasticity(const InputParameters & parameters)
103 : : ComputeStressBase(parameters),
104 246 : _nss(getParam<int>("nss")),
105 492 : _gprops(getParam<std::vector<Real>>("gprops")),
106 492 : _hprops(getParam<std::vector<Real>>("hprops")),
107 492 : _flowprops(getParam<std::vector<Real>>("flowprops")),
108 492 : _slip_sys_file_name(getParam<FileName>("slip_sys_file_name")),
109 492 : _slip_sys_res_prop_file_name(getParam<FileName>("slip_sys_res_prop_file_name")),
110 492 : _slip_sys_flow_prop_file_name(getParam<FileName>("slip_sys_flow_prop_file_name")),
111 492 : _slip_sys_hard_prop_file_name(getParam<FileName>("slip_sys_hard_prop_file_name")),
112 492 : _rtol(getParam<Real>("rtol")),
113 492 : _abs_tol(getParam<Real>("abs_tol")),
114 492 : _gtol(getParam<Real>("gtol")),
115 492 : _slip_incr_tol(getParam<Real>("slip_incr_tol")),
116 492 : _maxiter(getParam<unsigned int>("maxiter")),
117 492 : _maxiterg(getParam<unsigned int>("maxitergss")),
118 492 : _num_slip_sys_flowrate_props(getParam<unsigned int>("num_slip_sys_flowrate_props")),
119 492 : _tan_mod_type(getParam<MooseEnum>("tan_mod_type")),
120 492 : _intvar_read_type(getParam<MooseEnum>("intvar_read_type")),
121 492 : _num_slip_sys_props(getParam<unsigned int>("num_slip_sys_props")),
122 492 : _gen_rndm_stress_flag(getParam<bool>("gen_random_stress_flag")),
123 492 : _input_rndm_scale_var(getParam<bool>("input_random_scaling_var")),
124 492 : _rndm_scale_var(getParam<Real>("random_scaling_var")),
125 492 : _rndm_seed(getParam<unsigned int>("random_seed")),
126 492 : _max_substep_iter(getParam<unsigned int>("maximum_substep_iteration")),
127 492 : _use_line_search(getParam<bool>("use_line_search")),
128 492 : _min_lsrch_step(getParam<Real>("min_line_search_step_size")),
129 492 : _lsrch_tol(getParam<Real>("line_search_tol")),
130 492 : _lsrch_max_iter(getParam<unsigned int>("line_search_maxiter")),
131 492 : _lsrch_method(getParam<MooseEnum>("line_search_method")),
132 246 : _fp(declareProperty<RankTwoTensor>("fp")), // Plastic deformation gradient
133 492 : _fp_old(getMaterialPropertyOld<RankTwoTensor>(
134 : "fp")), // Plastic deformation gradient of previous increment
135 246 : _pk2(declareProperty<RankTwoTensor>("pk2")), // 2nd Piola-Kirchoff Stress
136 492 : _pk2_old(getMaterialPropertyOld<RankTwoTensor>(
137 : "pk2")), // 2nd Piola Kirchoff Stress of previous increment
138 246 : _lag_e(declareProperty<RankTwoTensor>("lage")), // Lagrangian strain
139 246 : _lag_e_old(
140 246 : getMaterialPropertyOld<RankTwoTensor>("lage")), // Lagrangian strain of previous increment
141 246 : _gss(declareProperty<std::vector<Real>>("gss")), // Slip system resistances
142 492 : _gss_old(getMaterialPropertyOld<std::vector<Real>>(
143 : "gss")), // Slip system resistances of previous increment
144 246 : _acc_slip(declareProperty<Real>("acc_slip")), // Accumulated slip
145 246 : _acc_slip_old(
146 246 : getMaterialPropertyOld<Real>("acc_slip")), // Accumulated slip of previous increment
147 246 : _update_rot(declareProperty<RankTwoTensor>(
148 : "update_rot")), // Rotation tensor considering material rotation and crystal orientation
149 492 : _deformation_gradient(getMaterialProperty<RankTwoTensor>("deformation_gradient")),
150 492 : _deformation_gradient_old(getMaterialPropertyOld<RankTwoTensor>("deformation_gradient")),
151 246 : _elasticity_tensor_name(_base_name + "elasticity_tensor"),
152 246 : _elasticity_tensor(getMaterialPropertyByName<RankFourTensor>(_elasticity_tensor_name)),
153 492 : _crysrot(getMaterialProperty<RankTwoTensor>("crysrot")),
154 246 : _mo(_nss * LIBMESH_DIM),
155 246 : _no(_nss * LIBMESH_DIM),
156 246 : _slip_incr(_nss),
157 246 : _tau(_nss),
158 246 : _dslipdtau(_nss),
159 246 : _s0(_nss),
160 246 : _gss_tmp(_nss),
161 246 : _gss_tmp_old(_nss),
162 984 : _dgss_dsliprate(_nss, _nss)
163 : {
164 246 : _err_tol = false;
165 :
166 246 : if (_num_slip_sys_props > 0)
167 21 : _slip_sys_props.resize(_nss * _num_slip_sys_props);
168 :
169 : _pk2_tmp.zero();
170 : _delta_dfgrd.zero();
171 :
172 246 : _first_step_iter = false;
173 246 : _last_step_iter = false;
174 : // Initialize variables in the first iteration of substepping
175 246 : _first_substep = true;
176 :
177 246 : _read_from_slip_sys_file = false;
178 246 : if (_intvar_read_type == "slip_sys_file")
179 21 : _read_from_slip_sys_file = true;
180 :
181 246 : 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 246 : getSlipSystems();
186 :
187 246 : RankTwoTensor::initRandom(_rndm_seed);
188 246 : }
189 :
190 : void
191 1120 : FiniteStrainCrystalPlasticity::initQpStatefulProperties()
192 : {
193 1120 : _stress[_qp].zero();
194 :
195 1120 : _fp[_qp].setToIdentity();
196 :
197 1120 : _pk2[_qp].zero();
198 1120 : _acc_slip[_qp] = 0.0;
199 1120 : _lag_e[_qp].zero();
200 :
201 1120 : _update_rot[_qp].setToIdentity();
202 :
203 1120 : initSlipSysProps(); // Initializes slip system related properties
204 1120 : initAdditionalProps();
205 1120 : }
206 :
207 : void
208 1120 : FiniteStrainCrystalPlasticity::initSlipSysProps()
209 : {
210 1120 : switch (_intvar_read_type)
211 : {
212 80 : case 0:
213 80 : assignSlipSysRes();
214 80 : break;
215 80 : case 1:
216 80 : readFileInitSlipSysRes();
217 80 : break;
218 960 : default:
219 960 : getInitSlipSysRes();
220 : }
221 :
222 1120 : if (_slip_sys_flow_prop_file_name.length() != 0)
223 80 : readFileFlowRateParams();
224 : else
225 1040 : getFlowRateParams();
226 :
227 1120 : if (_slip_sys_hard_prop_file_name.length() != 0)
228 0 : readFileHardnessParams();
229 : else
230 1120 : getHardnessParams();
231 1120 : }
232 :
233 : void
234 80 : FiniteStrainCrystalPlasticity::assignSlipSysRes()
235 : {
236 80 : _gss[_qp].resize(_nss);
237 :
238 1040 : for (unsigned int i = 0; i < _nss; ++i)
239 960 : _gss[_qp][i] = _slip_sys_props(i);
240 80 : }
241 :
242 : // Read initial slip system resistances from .txt file. See test.
243 : void
244 80 : FiniteStrainCrystalPlasticity::readFileInitSlipSysRes()
245 : {
246 80 : _gss[_qp].resize(_nss);
247 :
248 80 : MooseUtils::checkFileReadable(_slip_sys_res_prop_file_name);
249 :
250 80 : std::ifstream file;
251 80 : file.open(_slip_sys_res_prop_file_name.c_str());
252 :
253 1040 : for (unsigned int i = 0; i < _nss; ++i)
254 960 : if (!(file >> _gss[_qp][i]))
255 0 : mooseError("Error FiniteStrainCrystalPlasticity: Premature end of slip_sys_res_prop file");
256 :
257 80 : file.close();
258 80 : }
259 :
260 : // Read initial slip system resistances from .i file
261 : void
262 960 : FiniteStrainCrystalPlasticity::getInitSlipSysRes()
263 : {
264 960 : 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 960 : _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 3840 : for (unsigned int i = 0; i < _gprops.size() / num_data_grp; ++i)
274 : {
275 : Real vs, ve;
276 : unsigned int is, ie;
277 :
278 2880 : vs = _gprops[i * num_data_grp];
279 2880 : ve = _gprops[i * num_data_grp + 1];
280 :
281 2880 : 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 2880 : 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 2880 : is = static_cast<unsigned int>(vs);
293 2880 : ie = static_cast<unsigned int>(ve);
294 :
295 2880 : 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 14400 : for (unsigned int j = is; j <= ie; ++j)
303 11520 : _gss[_qp][j - 1] = _gprops[i * num_data_grp + 2];
304 : }
305 :
306 12480 : for (unsigned int i = 0; i < _nss; ++i)
307 11520 : if (_gss[_qp][i] <= 0.0)
308 0 : mooseError("FiniteStrainCrystalPLasticity: Value of resistance for slip system ",
309 0 : i + 1,
310 : " non positive");
311 960 : }
312 :
313 : // Read flow rate parameters from .txt file. See test.
314 : void
315 80 : FiniteStrainCrystalPlasticity::readFileFlowRateParams()
316 : {
317 80 : _a0.resize(_nss);
318 80 : _xm.resize(_nss);
319 :
320 80 : MooseUtils::checkFileReadable(_slip_sys_flow_prop_file_name);
321 :
322 80 : std::ifstream file;
323 80 : file.open(_slip_sys_flow_prop_file_name.c_str());
324 :
325 : std::vector<Real> vec;
326 80 : vec.resize(_num_slip_sys_flowrate_props);
327 :
328 1040 : for (unsigned int i = 0; i < _nss; ++i)
329 : {
330 2880 : for (unsigned int j = 0; j < _num_slip_sys_flowrate_props; ++j)
331 1920 : if (!(file >> vec[j]))
332 0 : mooseError(
333 : "Error FiniteStrainCrystalPlasticity: Premature end of slip_sys_flow_rate_param file");
334 :
335 960 : _a0(i) = vec[0];
336 960 : _xm(i) = vec[1];
337 : }
338 :
339 80 : file.close();
340 80 : }
341 :
342 : // Read flow rate parameters from .i file
343 : void
344 1040 : FiniteStrainCrystalPlasticity::getFlowRateParams()
345 : {
346 1040 : 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 1040 : _a0.resize(_nss);
351 1040 : _xm.resize(_nss);
352 :
353 1040 : 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 4160 : for (unsigned int i = 0; i < _flowprops.size() / num_data_grp; ++i)
358 : {
359 : Real vs, ve;
360 : unsigned int is, ie;
361 :
362 3120 : vs = _flowprops[i * num_data_grp];
363 3120 : ve = _flowprops[i * num_data_grp + 1];
364 :
365 3120 : 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 3120 : 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 3120 : is = static_cast<unsigned int>(vs);
377 3120 : ie = static_cast<unsigned int>(ve);
378 :
379 3120 : 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 15600 : for (unsigned int j = is; j <= ie; ++j)
387 : {
388 12480 : _a0(j - 1) = _flowprops[i * num_data_grp + 2];
389 12480 : _xm(j - 1) = _flowprops[i * num_data_grp + 3];
390 : }
391 : }
392 :
393 13520 : for (unsigned int i = 0; i < _nss; ++i)
394 : {
395 12480 : 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 1040 : }
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 1120 : FiniteStrainCrystalPlasticity::getHardnessParams()
413 : {
414 1120 : 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 1120 : _r = _hprops[0];
419 1120 : _h0 = _hprops[1];
420 1120 : _tau_init = _hprops[2];
421 1120 : _tau_sat = _hprops[3];
422 1120 : }
423 :
424 : // Read slip systems from file
425 : void
426 246 : FiniteStrainCrystalPlasticity::getSlipSystems()
427 : {
428 : Real vec[LIBMESH_DIM];
429 246 : std::ifstream fileslipsys;
430 :
431 246 : MooseUtils::checkFileReadable(_slip_sys_file_name);
432 :
433 246 : fileslipsys.open(_slip_sys_file_name.c_str());
434 :
435 3198 : for (unsigned int i = 0; i < _nss; ++i)
436 : {
437 : // Read the slip normal
438 11808 : for (const auto j : make_range(Moose::dim))
439 8856 : 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 2952 : mag = Utility::pow<2>(vec[0]) + Utility::pow<2>(vec[1]) + Utility::pow<2>(vec[2]);
445 2952 : mag = std::sqrt(mag);
446 :
447 11808 : for (unsigned j = 0; j < LIBMESH_DIM; ++j)
448 8856 : _no(i * LIBMESH_DIM + j) = vec[j] / mag;
449 :
450 : // Read the slip direction
451 11808 : for (const auto j : make_range(Moose::dim))
452 8856 : 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 2952 : mag = Utility::pow<2>(vec[0]) + Utility::pow<2>(vec[1]) + Utility::pow<2>(vec[2]);
457 2952 : mag = std::sqrt(mag);
458 :
459 11808 : for (const auto j : make_range(Moose::dim))
460 8856 : _mo(i * LIBMESH_DIM + j) = vec[j] / mag;
461 :
462 : mag = 0.0;
463 11808 : for (const auto j : make_range(Moose::dim))
464 8856 : mag += _mo(i * LIBMESH_DIM + j) * _no(i * LIBMESH_DIM + j);
465 :
466 2952 : 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 2952 : if (_read_from_slip_sys_file)
473 504 : for (unsigned int j = 0; j < _num_slip_sys_props; ++j)
474 252 : 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 246 : fileslipsys.close();
480 246 : }
481 :
482 : // Initialize addtional stateful material properties
483 : void
484 1120 : FiniteStrainCrystalPlasticity::initAdditionalProps()
485 : {
486 1120 : }
487 :
488 : /**
489 : * Solves stress residual equation using NR.
490 : * Updates slip system resistances iteratively.
491 : */
492 : void
493 85816 : 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 85816 : Real dt_original = _dt; // Stores original _dt; Reset at the end of solve
498 85816 : _first_substep = true; // Initialize variables at substep_iter = 1
499 :
500 85816 : if (_max_substep_iter > 1)
501 : {
502 11408 : _dfgrd_tmp_old = _deformation_gradient_old[_qp];
503 11408 : if (_dfgrd_tmp_old.det() == 0)
504 0 : _dfgrd_tmp_old.addIa(1.0);
505 :
506 11408 : _delta_dfgrd = _deformation_gradient[_qp] - _dfgrd_tmp_old;
507 11408 : _err_tol = true; // Indicator to continue substepping
508 : }
509 :
510 : // Substepping loop
511 123080 : while (_err_tol && _max_substep_iter > 1)
512 : {
513 37264 : _dt = dt_original / num_substep;
514 :
515 119064 : for (unsigned int istep = 0; istep < num_substep; ++istep)
516 : {
517 107656 : _first_step_iter = false;
518 107656 : if (istep == 0)
519 37264 : _first_step_iter = true;
520 :
521 107656 : _last_step_iter = false;
522 107656 : if (istep == num_substep - 1)
523 20560 : _last_step_iter = true;
524 :
525 107656 : _dfgrd_scale_factor = (static_cast<Real>(istep) + 1) / num_substep;
526 :
527 107656 : preSolveQp();
528 107656 : solveQp();
529 :
530 107656 : if (_err_tol)
531 : {
532 25856 : substep_iter++;
533 25856 : num_substep *= 2;
534 25856 : break;
535 : }
536 : }
537 :
538 37264 : _first_substep = false; // Prevents reinitialization
539 37264 : _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 37264 : if (!_err_tol || substep_iter > _max_substep_iter)
547 11408 : postSolveQp(); // Evaluate variables after successful solve or indicate failure
548 : }
549 :
550 : // No substepping
551 85816 : if (_max_substep_iter == 1)
552 : {
553 74408 : preSolveQp();
554 74408 : solveQp();
555 74408 : postSolveQp();
556 : }
557 85816 : }
558 :
559 : void
560 182064 : FiniteStrainCrystalPlasticity::preSolveQp()
561 : {
562 : // Initialize variable
563 182064 : if (_first_substep)
564 : {
565 85816 : _Jacobian_mult[_qp].zero(); // Initializes jacobian for preconditioner
566 85816 : calc_schmid_tensor();
567 : }
568 :
569 182064 : if (_max_substep_iter == 1)
570 74408 : _dfgrd_tmp = _deformation_gradient[_qp]; // Without substepping
571 : else
572 107656 : _dfgrd_tmp = _dfgrd_scale_factor * _delta_dfgrd + _dfgrd_tmp_old;
573 :
574 182064 : _err_tol = false;
575 182064 : }
576 :
577 : void
578 182064 : FiniteStrainCrystalPlasticity::solveQp()
579 : {
580 182064 : preSolveStatevar();
581 182064 : solveStatevar();
582 182064 : if (_err_tol)
583 : return;
584 144344 : postSolveStatevar();
585 : }
586 :
587 : void
588 85816 : FiniteStrainCrystalPlasticity::postSolveQp()
589 : {
590 85816 : if (_err_tol)
591 : {
592 11864 : _err_tol = false;
593 11864 : if (_gen_rndm_stress_flag)
594 : {
595 11864 : if (!_input_rndm_scale_var)
596 11864 : _rndm_scale_var = _elasticity_tensor[_qp](0, 0, 0, 0);
597 :
598 11864 : _stress[_qp] = RankTwoTensor::genRandomSymmTensor(_rndm_scale_var, 1.0);
599 : }
600 : else
601 0 : mooseError("FiniteStrainCrystalPlasticity: Constitutive failure");
602 : }
603 : else
604 : {
605 73952 : _stress[_qp] = _fe * _pk2[_qp] * _fe.transpose() / _fe.det();
606 :
607 73952 : _Jacobian_mult[_qp] += calcTangentModuli(); // Calculate jacobian for preconditioner
608 :
609 73952 : RankTwoTensor iden(RankTwoTensor::initIdentity);
610 :
611 73952 : _lag_e[_qp] = _deformation_gradient[_qp].transpose() * _deformation_gradient[_qp] - iden;
612 73952 : _lag_e[_qp] = _lag_e[_qp] * 0.5;
613 :
614 73952 : RankTwoTensor rot;
615 73952 : rot = get_current_rotation(_deformation_gradient[_qp]); // Calculate material rotation
616 73952 : _update_rot[_qp] = rot * _crysrot[_qp];
617 : }
618 85816 : }
619 :
620 : void
621 182064 : FiniteStrainCrystalPlasticity::preSolveStatevar()
622 : {
623 182064 : if (_max_substep_iter == 1) // No substepping
624 : {
625 74408 : _gss_tmp = _gss_old[_qp];
626 74408 : _accslip_tmp_old = _acc_slip_old[_qp];
627 : }
628 : else
629 : {
630 107656 : if (_first_step_iter)
631 : {
632 37264 : _gss_tmp = _gss_tmp_old = _gss_old[_qp];
633 37264 : _accslip_tmp_old = _acc_slip_old[_qp];
634 : }
635 : else
636 70392 : _gss_tmp = _gss_tmp_old;
637 : }
638 182064 : }
639 :
640 : void
641 74760 : FiniteStrainCrystalPlasticity::solveStatevar()
642 : {
643 : Real gmax, gdiff;
644 : unsigned int iterg;
645 74760 : std::vector<Real> gss_prev(_nss);
646 :
647 74760 : gmax = 1.1 * _gtol;
648 : iterg = 0;
649 :
650 154744 : while (gmax > _gtol && iterg < _maxiterg) // Check for slip system resistance update tolerance
651 : {
652 92800 : preSolveStress();
653 92800 : solveStress();
654 92800 : if (_err_tol)
655 : return;
656 79984 : postSolveStress();
657 :
658 79984 : gss_prev = _gss_tmp;
659 :
660 79984 : update_slip_system_resistance(); // Update slip system resistance
661 :
662 : gmax = 0.0;
663 1039792 : for (unsigned i = 0; i < _nss; ++i)
664 : {
665 959808 : gdiff = std::abs(gss_prev[i] - _gss_tmp[i]); // Calculate increment size
666 :
667 959808 : if (gdiff > gmax)
668 : gmax = gdiff;
669 : }
670 79984 : iterg++;
671 : }
672 :
673 61944 : 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 74760 : }
681 :
682 : void
683 144344 : FiniteStrainCrystalPlasticity::postSolveStatevar()
684 : {
685 144344 : if (_max_substep_iter == 1) // No substepping
686 : {
687 62544 : _gss[_qp] = _gss_tmp;
688 62544 : _acc_slip[_qp] = _accslip_tmp;
689 : }
690 : else
691 : {
692 81800 : if (_last_step_iter)
693 : {
694 11408 : _gss[_qp] = _gss_tmp;
695 11408 : _acc_slip[_qp] = _accslip_tmp;
696 : }
697 : else
698 : {
699 70392 : _gss_tmp_old = _gss_tmp;
700 70392 : _accslip_tmp_old = _accslip_tmp;
701 : }
702 : }
703 144344 : }
704 :
705 : void
706 200104 : FiniteStrainCrystalPlasticity::preSolveStress()
707 : {
708 200104 : if (_max_substep_iter == 1) // No substepping
709 : {
710 89312 : _pk2_tmp = _pk2_old[_qp];
711 89312 : _fp_old_inv = _fp_old[_qp].inverse();
712 89312 : _fp_inv = _fp_old_inv;
713 89312 : _fp_prev_inv = _fp_inv;
714 : }
715 : else
716 : {
717 110792 : if (_first_step_iter)
718 : {
719 39448 : _pk2_tmp = _pk2_tmp_old = _pk2_old[_qp];
720 39448 : _fp_old_inv = _fp_old[_qp].inverse();
721 : }
722 : else
723 71344 : _pk2_tmp = _pk2_tmp_old;
724 :
725 110792 : _fp_inv = _fp_old_inv;
726 110792 : _fp_prev_inv = _fp_inv;
727 : }
728 200104 : }
729 :
730 : void
731 92800 : FiniteStrainCrystalPlasticity::solveStress()
732 : {
733 : unsigned int iter = 0;
734 92800 : RankTwoTensor resid, dpk2;
735 92800 : RankFourTensor jac;
736 : Real rnorm, rnorm0, rnorm_prev;
737 :
738 92800 : calc_resid_jacob(resid, jac); // Calculate stress residual
739 92800 : 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 92800 : rnorm = resid.L2norm();
752 : rnorm0 = rnorm;
753 :
754 485348 : while (rnorm > _rtol * rnorm0 && rnorm0 > _abs_tol &&
755 405364 : iter < _maxiter) // Check for stress residual tolerance
756 : {
757 405364 : dpk2 = -jac.invSymm() * resid; // Calculate stress increment
758 405364 : _pk2_tmp = _pk2_tmp + dpk2; // Update stress
759 405364 : calc_resid_jacob(resid, jac);
760 405364 : internalVariableUpdateNRiteration(); // update _fp_prev_inv
761 :
762 405364 : 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 392548 : rnorm = resid.L2norm();
776 :
777 392548 : 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 392548 : if (_use_line_search)
787 16816 : rnorm = resid.L2norm();
788 :
789 392548 : iter++;
790 : }
791 :
792 79984 : 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 162384 : FiniteStrainCrystalPlasticity::postSolveStress()
803 : {
804 162384 : if (_max_substep_iter == 1) // No substepping
805 : {
806 77448 : _fp[_qp] = _fp_inv.inverse();
807 77448 : _pk2[_qp] = _pk2_tmp;
808 : }
809 : else
810 : {
811 84936 : if (_last_step_iter)
812 : {
813 13592 : _fp[_qp] = _fp_inv.inverse();
814 13592 : _pk2[_qp] = _pk2_tmp;
815 : }
816 : else
817 : {
818 71344 : _fp_old_inv = _fp_inv;
819 71344 : _pk2_tmp_old = _pk2_tmp;
820 : }
821 : }
822 162384 : }
823 :
824 : // Update slip system resistance. Overide to incorporate new slip system resistance laws
825 : void
826 892560 : FiniteStrainCrystalPlasticity::update_slip_system_resistance()
827 : {
828 892560 : updateGss();
829 892560 : }
830 :
831 : /**
832 : * Old function to update slip system resistances.
833 : * Kept to avoid code break at computeQpstress
834 : */
835 : void
836 892560 : FiniteStrainCrystalPlasticity::updateGss()
837 : {
838 892560 : DenseVector<Real> hb(_nss);
839 : Real qab;
840 :
841 892560 : Real a = _hprops[4]; // Kalidindi
842 :
843 892560 : _accslip_tmp = _accslip_tmp_old;
844 11603280 : for (unsigned int i = 0; i < _nss; ++i)
845 10710720 : _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 11603280 : for (unsigned int i = 0; i < _nss; ++i)
851 : // hb(i)=val;
852 10710720 : hb(i) = _h0 * std::pow(std::abs(1.0 - _gss_tmp[i] / _tau_sat), a) *
853 10710720 : std::copysign(1.0, 1.0 - _gss_tmp[i] / _tau_sat);
854 :
855 11603280 : for (unsigned int i = 0; i < _nss; ++i)
856 : {
857 10710720 : if (_max_substep_iter == 1) // No substepping
858 1218432 : _gss_tmp[i] = _gss_old[_qp][i];
859 : else
860 9492288 : _gss_tmp[i] = _gss_tmp_old[i];
861 :
862 139239360 : for (unsigned int j = 0; j < _nss; ++j)
863 : {
864 : unsigned int iplane, jplane;
865 128528640 : iplane = i / 3;
866 128528640 : jplane = j / 3;
867 :
868 128528640 : if (iplane == jplane) // Kalidindi
869 : qab = 1.0;
870 : else
871 96396480 : qab = _r;
872 :
873 128528640 : _gss_tmp[i] += qab * hb(j) * std::abs(_slip_incr(j));
874 128528640 : _dgss_dsliprate(i, j) = qab * hb(j) * std::copysign(1.0, _slip_incr(j)) * _dt;
875 : }
876 : }
877 892560 : }
878 :
879 : // Calculates stress residual equation and jacobian
880 : void
881 498164 : FiniteStrainCrystalPlasticity::calc_resid_jacob(RankTwoTensor & resid, RankFourTensor & jac)
882 : {
883 498164 : calcResidual(resid);
884 498164 : if (_err_tol)
885 : return;
886 485348 : calcJacobian(jac);
887 : }
888 :
889 : void
890 498164 : FiniteStrainCrystalPlasticity::calcResidual(RankTwoTensor & resid)
891 : {
892 498164 : RankTwoTensor iden(RankTwoTensor::initIdentity), ce, ee, ce_pk2, eqv_slip_incr, pk2_new;
893 :
894 498164 : _fe = _dfgrd_tmp * _fp_prev_inv; // _fp_inv ==> _fp_prev_inv
895 :
896 498164 : ce = _fe.transpose() * _fe;
897 498164 : ce_pk2 = ce * _pk2_tmp;
898 498164 : ce_pk2 = ce_pk2 / _fe.det();
899 :
900 : // Calculate Schmid tensor and resolved shear stresses
901 6476132 : for (unsigned int i = 0; i < _nss; ++i)
902 5977968 : _tau(i) = ce_pk2.doubleContraction(_s0[i]);
903 :
904 498164 : getSlipIncrements(); // Calculate dslip,dslipdtau
905 :
906 498164 : if (_err_tol)
907 : return;
908 :
909 : eqv_slip_incr.zero();
910 6309524 : for (unsigned int i = 0; i < _nss; ++i)
911 5824176 : eqv_slip_incr += _s0[i] * _slip_incr(i);
912 :
913 485348 : eqv_slip_incr = iden - eqv_slip_incr;
914 485348 : _fp_inv = _fp_old_inv * eqv_slip_incr;
915 485348 : _fe = _dfgrd_tmp * _fp_inv;
916 :
917 485348 : ce = _fe.transpose() * _fe;
918 485348 : ee = ce - iden;
919 485348 : ee *= 0.5;
920 :
921 485348 : pk2_new = _elasticity_tensor[_qp] * ee;
922 :
923 485348 : resid = _pk2_tmp - pk2_new;
924 : }
925 :
926 : void
927 485348 : FiniteStrainCrystalPlasticity::calcJacobian(RankFourTensor & jac)
928 : {
929 485348 : RankFourTensor dfedfpinv, deedfe, dfpinvdpk2;
930 :
931 485348 : std::vector<RankTwoTensor> dtaudpk2(_nss), dfpinvdslip(_nss);
932 :
933 6309524 : for (unsigned int i = 0; i < _nss; ++i)
934 : {
935 5824176 : dtaudpk2[i] = _s0[i];
936 5824176 : dfpinvdslip[i] = -_fp_old_inv * _s0[i];
937 : }
938 :
939 1941392 : for (const auto i : make_range(Moose::dim))
940 5824176 : for (const auto j : make_range(Moose::dim))
941 17472528 : for (const auto k : make_range(Moose::dim))
942 13104396 : dfedfpinv(i, j, k, j) = _dfgrd_tmp(i, k);
943 :
944 1941392 : for (const auto i : make_range(Moose::dim))
945 5824176 : for (const auto j : make_range(Moose::dim))
946 17472528 : for (const auto k : make_range(Moose::dim))
947 : {
948 13104396 : deedfe(i, j, k, i) = deedfe(i, j, k, i) + _fe(k, j) * 0.5;
949 13104396 : deedfe(i, j, k, j) = deedfe(i, j, k, j) + _fe(k, i) * 0.5;
950 : }
951 :
952 6309524 : for (unsigned int i = 0; i < _nss; ++i)
953 5824176 : dfpinvdpk2 += (dfpinvdslip[i] * _dslipdtau(i)).outerProduct(dtaudpk2[i]);
954 :
955 : jac =
956 970696 : RankFourTensor::IdentityFour() - (_elasticity_tensor[_qp] * deedfe * dfedfpinv * dfpinvdpk2);
957 485348 : }
958 :
959 : // Calculate slip increment,dslipdtau. Override to modify.
960 : void
961 1310740 : FiniteStrainCrystalPlasticity::getSlipIncrements()
962 : {
963 16586980 : for (unsigned int i = 0; i < _nss; ++i)
964 : {
965 15313960 : _slip_incr(i) = _a0(i) * std::pow(std::abs(_tau(i) / _gss_tmp[i]), 1.0 / _xm(i)) *
966 15313960 : std::copysign(1.0, _tau(i)) * _dt;
967 15313960 : if (std::abs(_slip_incr(i)) > _slip_incr_tol)
968 : {
969 37720 : _err_tol = true;
970 : #ifdef DEBUG
971 : mooseWarning("Maximum allowable slip increment exceeded ", std::abs(_slip_incr(i)));
972 : #endif
973 37720 : return;
974 : }
975 : }
976 :
977 16549260 : for (unsigned int i = 0; i < _nss; ++i)
978 15276240 : _dslipdtau(i) = _a0(i) / _xm(i) *
979 15276240 : std::pow(std::abs(_tau(i) / _gss_tmp[i]), 1.0 / _xm(i) - 1.0) / _gss_tmp[i] *
980 15276240 : _dt;
981 : }
982 :
983 : // Calls getMatRot to perform RU factorization of a tensor.
984 : RankTwoTensor
985 73952 : FiniteStrainCrystalPlasticity::get_current_rotation(const RankTwoTensor & a)
986 : {
987 73952 : return getMatRot(a);
988 : }
989 :
990 : // Performs RU factorization of a tensor
991 : RankTwoTensor
992 73952 : FiniteStrainCrystalPlasticity::getMatRot(const RankTwoTensor & a)
993 : {
994 73952 : RankTwoTensor rot;
995 73952 : RankTwoTensor c, diag, evec;
996 : PetscScalar cmat[LIBMESH_DIM][LIBMESH_DIM], work[10];
997 : PetscReal w[LIBMESH_DIM];
998 73952 : PetscBLASInt nd = LIBMESH_DIM, lwork = 10, info;
999 :
1000 73952 : c = a.transpose() * a;
1001 :
1002 295808 : for (const auto i : make_range(Moose::dim))
1003 887424 : for (const auto j : make_range(Moose::dim))
1004 665568 : cmat[i][j] = c(i, j);
1005 :
1006 73952 : LAPACKsyev_("V", "U", &nd, &cmat[0][0], &nd, w, work, &lwork, &info);
1007 :
1008 73952 : if (info != 0)
1009 0 : mooseError("FiniteStrainCrystalPLasticity: DSYEV function call in getMatRot function failed");
1010 :
1011 : diag.zero();
1012 :
1013 295808 : for (const auto i : make_range(Moose::dim))
1014 221856 : diag(i, i) = std::sqrt(w[i]);
1015 :
1016 295808 : for (const auto i : make_range(Moose::dim))
1017 887424 : for (const auto j : make_range(Moose::dim))
1018 665568 : evec(i, j) = cmat[i][j];
1019 :
1020 73952 : rot = a * ((evec.transpose() * diag * evec).inverse());
1021 :
1022 73952 : 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 73952 : FiniteStrainCrystalPlasticity::calcTangentModuli()
1033 : {
1034 73952 : RankFourTensor tan_mod;
1035 :
1036 73952 : switch (_tan_mod_type)
1037 : {
1038 68048 : case 0:
1039 68048 : tan_mod = elastoPlasticTangentModuli();
1040 68048 : break;
1041 5904 : default:
1042 5904 : tan_mod = elasticTangentModuli();
1043 : }
1044 :
1045 73952 : return tan_mod;
1046 : }
1047 :
1048 : void
1049 85816 : FiniteStrainCrystalPlasticity::calc_schmid_tensor()
1050 : {
1051 85816 : DenseVector<Real> mo(LIBMESH_DIM * _nss), no(LIBMESH_DIM * _nss);
1052 :
1053 : // Update slip direction and normal with crystal orientation
1054 1115608 : for (unsigned int i = 0; i < _nss; ++i)
1055 : {
1056 4119168 : for (const auto j : make_range(Moose::dim))
1057 : {
1058 3089376 : mo(i * LIBMESH_DIM + j) = 0.0;
1059 12357504 : for (const auto k : make_range(Moose::dim))
1060 9268128 : mo(i * LIBMESH_DIM + j) =
1061 9268128 : mo(i * LIBMESH_DIM + j) + _crysrot[_qp](j, k) * _mo(i * LIBMESH_DIM + k);
1062 : }
1063 :
1064 4119168 : for (const auto j : make_range(Moose::dim))
1065 : {
1066 3089376 : no(i * LIBMESH_DIM + j) = 0.0;
1067 12357504 : for (const auto k : make_range(Moose::dim))
1068 9268128 : no(i * LIBMESH_DIM + j) =
1069 9268128 : 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 1115608 : for (unsigned int i = 0; i < _nss; ++i)
1075 4119168 : for (const auto j : make_range(Moose::dim))
1076 12357504 : for (const auto k : make_range(Moose::dim))
1077 9268128 : _s0[i](j, k) = mo(i * LIBMESH_DIM + j) * no(i * LIBMESH_DIM + k);
1078 85816 : }
1079 :
1080 : RankFourTensor
1081 68048 : FiniteStrainCrystalPlasticity::elastoPlasticTangentModuli()
1082 : {
1083 68048 : RankFourTensor tan_mod;
1084 68048 : RankTwoTensor pk2fet, fepk2;
1085 68048 : RankFourTensor deedfe, dsigdpk2dfe;
1086 :
1087 : // Fill in the matrix stiffness material property
1088 :
1089 272192 : for (const auto i : make_range(Moose::dim))
1090 816576 : for (const auto j : make_range(Moose::dim))
1091 2449728 : for (const auto k : make_range(Moose::dim))
1092 : {
1093 1837296 : deedfe(i, j, k, i) = deedfe(i, j, k, i) + _fe(k, j) * 0.5;
1094 1837296 : deedfe(i, j, k, j) = deedfe(i, j, k, j) + _fe(k, i) * 0.5;
1095 : }
1096 :
1097 : usingTensorIndices(i_, j_, k_, l_);
1098 68048 : dsigdpk2dfe = _fe.times<i_, k_, j_, l_>(_fe) * _elasticity_tensor[_qp] * deedfe;
1099 :
1100 68048 : pk2fet = _pk2_tmp * _fe.transpose();
1101 68048 : fepk2 = _fe * _pk2_tmp;
1102 :
1103 272192 : for (const auto i : make_range(Moose::dim))
1104 816576 : for (const auto j : make_range(Moose::dim))
1105 2449728 : for (const auto l : make_range(Moose::dim))
1106 : {
1107 1837296 : tan_mod(i, j, i, l) = tan_mod(i, j, i, l) + pk2fet(l, j);
1108 1837296 : tan_mod(i, j, j, l) = tan_mod(i, j, j, l) + fepk2(i, l);
1109 : }
1110 :
1111 68048 : tan_mod += dsigdpk2dfe;
1112 :
1113 68048 : Real je = _fe.det();
1114 68048 : if (je > 0.0)
1115 68048 : tan_mod /= je;
1116 :
1117 68048 : return tan_mod;
1118 : }
1119 :
1120 : RankFourTensor
1121 5904 : FiniteStrainCrystalPlasticity::elasticTangentModuli()
1122 : {
1123 5904 : 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 405364 : FiniteStrainCrystalPlasticity::internalVariableUpdateNRiteration()
1211 : {
1212 405364 : _fp_prev_inv = _fp_inv; // update _fp_prev_inv
1213 405364 : }
|
