www.mooseframework.org
FiniteStrainUObasedCP.C
Go to the documentation of this file.
1 /****************************************************************/
2 /* MOOSE - Multiphysics Object Oriented Simulation Environment */
3 /* */
4 /* All contents are licensed under LGPL V2.1 */
5 /* See LICENSE for full restrictions */
6 /****************************************************************/
8 #include "petscblaslapack.h"
9 #include "MooseException.h"
14 
15 template <>
16 InputParameters
18 {
19  InputParameters params = validParams<ComputeStressBase>();
20  params.addClassDescription(
21  "Crystal Plasticity base class: FCC system with power law flow rule implemented");
22  params.addParam<Real>("rtol", 1e-6, "Constitutive stress residue relative tolerance");
23  params.addParam<Real>("abs_tol", 1e-6, "Constitutive stress residue absolute tolerance");
24  params.addParam<Real>(
25  "stol", 1e-2, "Constitutive slip system resistance relative residual tolerance");
26  params.addParam<Real>(
27  "zero_tol", 1e-12, "Tolerance for residual check when variable value is zero");
28  params.addParam<unsigned int>("maxiter", 100, "Maximum number of iterations for stress update");
29  params.addParam<unsigned int>(
30  "maxiter_state_variable", 100, "Maximum number of iterations for state variable update");
31  MooseEnum tan_mod_options("exact none", "none"); // Type of read
32  params.addParam<MooseEnum>("tan_mod_type",
33  tan_mod_options,
34  "Type of tangent moduli for preconditioner: default elastic");
35  params.addParam<unsigned int>(
36  "maximum_substep_iteration", 1, "Maximum number of substep iteration");
37  params.addParam<bool>("use_line_search", false, "Use line search in constitutive update");
38  params.addParam<Real>("min_line_search_step_size", 0.01, "Minimum line search step size");
39  params.addParam<Real>("line_search_tol", 0.5, "Line search bisection method tolerance");
40  params.addParam<unsigned int>(
41  "line_search_maxiter", 20, "Line search bisection method maximum number of iteration");
42  MooseEnum line_search_method("CUT_HALF BISECTION", "CUT_HALF");
43  params.addParam<MooseEnum>(
44  "line_search_method", line_search_method, "The method used in line search");
45  params.addRequiredParam<std::vector<UserObjectName>>(
46  "uo_slip_rates",
47  "List of names of user objects that define the slip rates for this material.");
48  params.addRequiredParam<std::vector<UserObjectName>>(
49  "uo_slip_resistances",
50  "List of names of user objects that define the slip resistances for this material.");
51  params.addRequiredParam<std::vector<UserObjectName>>(
52  "uo_state_vars",
53  "List of names of user objects that define the state variable for this material.");
54  params.addRequiredParam<std::vector<UserObjectName>>(
55  "uo_state_var_evol_rate_comps",
56  "List of names of user objects that define the state "
57  "variable evolution rate components for this material.");
58  return params;
59 }
60 
61 FiniteStrainUObasedCP::FiniteStrainUObasedCP(const InputParameters & parameters)
62  : ComputeStressBase(parameters),
63  _num_uo_slip_rates(parameters.get<std::vector<UserObjectName>>("uo_slip_rates").size()),
64  _num_uo_slip_resistances(
65  parameters.get<std::vector<UserObjectName>>("uo_slip_resistances").size()),
66  _num_uo_state_vars(parameters.get<std::vector<UserObjectName>>("uo_state_vars").size()),
67  _num_uo_state_var_evol_rate_comps(
68  parameters.get<std::vector<UserObjectName>>("uo_state_var_evol_rate_comps").size()),
69  _rtol(getParam<Real>("rtol")),
70  _abs_tol(getParam<Real>("abs_tol")),
71  _stol(getParam<Real>("stol")),
72  _zero_tol(getParam<Real>("zero_tol")),
73  _maxiter(getParam<unsigned int>("maxiter")),
74  _maxiterg(getParam<unsigned int>("maxiter_state_variable")),
75  _tan_mod_type(getParam<MooseEnum>("tan_mod_type")),
76  _max_substep_iter(getParam<unsigned int>("maximum_substep_iteration")),
77  _use_line_search(getParam<bool>("use_line_search")),
78  _min_lsrch_step(getParam<Real>("min_line_search_step_size")),
79  _lsrch_tol(getParam<Real>("line_search_tol")),
80  _lsrch_max_iter(getParam<unsigned int>("line_search_maxiter")),
81  _lsrch_method(getParam<MooseEnum>("line_search_method")),
82  _fp(declareProperty<RankTwoTensor>("fp")), // Plastic deformation gradient
83  _fp_old(getMaterialPropertyOld<RankTwoTensor>(
84  "fp")), // Plastic deformation gradient of previous increment
85  _pk2(declareProperty<RankTwoTensor>("pk2")), // 2nd Piola Kirchoff Stress
86  _pk2_old(getMaterialPropertyOld<RankTwoTensor>(
87  "pk2")), // 2nd Piola Kirchoff Stress of previous increment
88  _lag_e(declareProperty<RankTwoTensor>("lage")), // Lagrangian strain
89  _update_rot(declareProperty<RankTwoTensor>(
90  "update_rot")), // Rotation tensor considering material rotation and crystal orientation
91  _update_rot_old(getMaterialPropertyOld<RankTwoTensor>("update_rot")),
92  _deformation_gradient(getMaterialProperty<RankTwoTensor>("deformation_gradient")),
93  _deformation_gradient_old(getMaterialPropertyOld<RankTwoTensor>("deformation_gradient")),
94  _crysrot(getMaterialProperty<RankTwoTensor>("crysrot"))
95 {
96  _err_tol = false;
97 
98  _delta_dfgrd.zero();
99 
100  // resize the material properties for each userobject
106 
107  // resize the flow direction
109 
110  // resize local state variables
113 
114  // resize user objects
119 
120  // assign the user objects
121  for (unsigned int i = 0; i < _num_uo_slip_rates; ++i)
122  {
123  _uo_slip_rates[i] = &getUserObjectByName<CrystalPlasticitySlipRate>(
124  parameters.get<std::vector<UserObjectName>>("uo_slip_rates")[i]);
125  _mat_prop_slip_rates[i] = &declareProperty<std::vector<Real>>(
126  parameters.get<std::vector<UserObjectName>>("uo_slip_rates")[i]);
127  _flow_direction[i] = &declareProperty<std::vector<RankTwoTensor>>(
128  parameters.get<std::vector<UserObjectName>>("uo_slip_rates")[i] + "_flow_direction");
129  }
130 
131  for (unsigned int i = 0; i < _num_uo_slip_resistances; ++i)
132  {
133  _uo_slip_resistances[i] = &getUserObjectByName<CrystalPlasticitySlipResistance>(
134  parameters.get<std::vector<UserObjectName>>("uo_slip_resistances")[i]);
135  _mat_prop_slip_resistances[i] = &declareProperty<std::vector<Real>>(
136  parameters.get<std::vector<UserObjectName>>("uo_slip_resistances")[i]);
137  }
138 
139  for (unsigned int i = 0; i < _num_uo_state_vars; ++i)
140  {
141  _uo_state_vars[i] = &getUserObjectByName<CrystalPlasticityStateVariable>(
142  parameters.get<std::vector<UserObjectName>>("uo_state_vars")[i]);
143  _mat_prop_state_vars[i] = &declareProperty<std::vector<Real>>(
144  parameters.get<std::vector<UserObjectName>>("uo_state_vars")[i]);
145  _mat_prop_state_vars_old[i] = &getMaterialPropertyOld<std::vector<Real>>(
146  parameters.get<std::vector<UserObjectName>>("uo_state_vars")[i]);
147  }
148 
149  for (unsigned int i = 0; i < _num_uo_state_var_evol_rate_comps; ++i)
150  {
151  _uo_state_var_evol_rate_comps[i] = &getUserObjectByName<CrystalPlasticityStateVarRateComponent>(
152  parameters.get<std::vector<UserObjectName>>("uo_state_var_evol_rate_comps")[i]);
153  _mat_prop_state_var_evol_rate_comps[i] = &declareProperty<std::vector<Real>>(
154  parameters.get<std::vector<UserObjectName>>("uo_state_var_evol_rate_comps")[i]);
155  }
156 }
157 
158 void
160 {
161  for (unsigned int i = 0; i < _num_uo_slip_rates; ++i)
162  {
163  (*_mat_prop_slip_rates[i])[_qp].resize(_uo_slip_rates[i]->variableSize());
164  (*_flow_direction[i])[_qp].resize(_uo_slip_rates[i]->variableSize());
165  }
166 
167  for (unsigned int i = 0; i < _num_uo_slip_resistances; ++i)
168  (*_mat_prop_slip_resistances[i])[_qp].resize(_uo_slip_resistances[i]->variableSize());
169 
170  for (unsigned int i = 0; i < _num_uo_state_vars; ++i)
171  {
172  (*_mat_prop_state_vars[i])[_qp].resize(_uo_state_vars[i]->variableSize());
173  // TODO: remove this nasty const_cast if you can figure out how
174  const_cast<MaterialProperty<std::vector<Real>> &>(*_mat_prop_state_vars_old[i])[_qp].resize(
175  _uo_state_vars[i]->variableSize());
176  _state_vars_old[i].resize(_uo_state_vars[i]->variableSize());
177  _state_vars_prev[i].resize(_uo_state_vars[i]->variableSize());
178  }
179 
180  for (unsigned int i = 0; i < _num_uo_state_var_evol_rate_comps; ++i)
181  (*_mat_prop_state_var_evol_rate_comps[i])[_qp].resize(
182  _uo_state_var_evol_rate_comps[i]->variableSize());
183 
184  _stress[_qp].zero();
185 
186  _fp[_qp].zero();
187  _fp[_qp].addIa(1.0);
188 
189  _pk2[_qp].zero();
190 
191  _lag_e[_qp].zero();
192 
193  _update_rot[_qp].zero();
194  _update_rot[_qp].addIa(1.0);
195 
196  for (unsigned int i = 0; i < _num_uo_state_vars; ++i)
197  {
198  // Initializes slip system related properties
199  _uo_state_vars[i]->initSlipSysProps((*_mat_prop_state_vars[i])[_qp], _q_point[_qp]);
200  // TODO: remove this nasty const_cast if you can figure out how
201  const_cast<MaterialProperty<std::vector<Real>> &>(*_mat_prop_state_vars_old[i])[_qp] =
202  (*_mat_prop_state_vars[i])[_qp];
203  }
204 }
205 
210 void
212 {
213  // Userobject based crystal plasticity does not support face/boundary material property
214  // calculation.
215  if (isBoundaryMaterial())
216  return;
217  // Depth of substepping; Limited to maximum substep iteration
218  unsigned int substep_iter = 1;
219  // Calculated from substep_iter as 2^substep_iter
220  unsigned int num_substep = 1;
221  // Store original _dt; Reset at the end of solve
222  Real dt_original = _dt;
223 
225  if (_dfgrd_tmp_old.det() == 0)
226  _dfgrd_tmp_old.addIa(1.0);
227 
229 
230  // Saves the old stateful properties that is modified during sub stepping
231  for (unsigned int i = 0; i < _num_uo_state_vars; ++i)
233 
234  for (unsigned int i = 0; i < _num_uo_slip_rates; ++i)
235  _uo_slip_rates[i]->calcFlowDirection(_qp, (*_flow_direction[i])[_qp]);
236 
237  do
238  {
239  _err_tol = false;
240 
241  preSolveQp();
242 
243  _dt = dt_original / num_substep;
244 
245  for (unsigned int istep = 0; istep < num_substep; ++istep)
246  {
247  _dfgrd_tmp = (static_cast<Real>(istep) + 1) / num_substep * _delta_dfgrd + _dfgrd_tmp_old;
248 
249  solveQp();
250 
251  if (_err_tol)
252  {
253  substep_iter++;
254  num_substep *= 2;
255  break;
256  }
257  }
258  if (substep_iter > _max_substep_iter && _err_tol)
259  throw MooseException("FiniteStrainUObasedCP: Constitutive failure.");
260  } while (_err_tol);
261 
262  _dt = dt_original;
263 
264  postSolveQp();
265 }
266 
267 void
269 {
270  // TODO: remove this nasty const_cast if you can figure out how
271  for (unsigned int i = 0; i < _num_uo_state_vars; ++i)
272  (*_mat_prop_state_vars[i])[_qp] =
273  const_cast<MaterialProperty<std::vector<Real>> &>(*_mat_prop_state_vars_old[i])[_qp] =
274  _state_vars_old[i];
275 
276  _pk2[_qp] = _pk2_old[_qp];
277  _fp_old_inv = _fp_old[_qp].inverse();
278 }
279 
280 void
282 {
284  solveStatevar();
285  if (_err_tol)
286  return;
288 }
289 
290 void
292 {
293  // TODO: remove this nasty const_cast if you can figure out how
294  // Restores the the old stateful properties after a successful solve
295  for (unsigned int i = 0; i < _num_uo_state_vars; ++i)
296  const_cast<MaterialProperty<std::vector<Real>> &>(*_mat_prop_state_vars_old[i])[_qp] =
297  _state_vars_old[i];
298 
299  _stress[_qp] = _fe * _pk2[_qp] * _fe.transpose() / _fe.det();
300 
301  // Calculate jacobian for preconditioner
303 
304  RankTwoTensor iden;
305  iden.addIa(1.0);
306 
307  _lag_e[_qp] = _deformation_gradient[_qp].transpose() * _deformation_gradient[_qp] - iden;
308  _lag_e[_qp] = _lag_e[_qp] * 0.5;
309 
310  RankTwoTensor rot;
311  // Calculate material rotation
312  _deformation_gradient[_qp].getRUDecompositionRotation(rot);
313  _update_rot[_qp] = rot * _crysrot[_qp];
314 }
315 
316 void
318 {
319  for (unsigned int i = 0; i < _num_uo_state_vars; ++i)
320  (*_mat_prop_state_vars[i])[_qp] = (*_mat_prop_state_vars_old[i])[_qp];
321 
322  for (unsigned int i = 0; i < _num_uo_slip_resistances; ++i)
323  _uo_slip_resistances[i]->calcSlipResistance(_qp, (*_mat_prop_slip_resistances[i])[_qp]);
324 
326 }
327 
328 void
330 {
331  unsigned int iterg;
332  bool iter_flag = true;
333 
334  iterg = 0;
335  // Check for slip system resistance update tolerance
336  while (iter_flag && iterg < _maxiterg)
337  {
338  preSolveStress();
339  solveStress();
340  if (_err_tol)
341  return;
342  postSolveStress();
343 
344  // Update slip system resistance and state variable
346 
347  if (_err_tol)
348  return;
349 
350  iter_flag = isStateVariablesConverged();
351  iterg++;
352  }
353 
354  if (iterg == _maxiterg)
355  {
356 #ifdef DEBUG
357  mooseWarning("FiniteStrainUObasedCP: Hardness Integration error\n");
358 #endif
359  _err_tol = true;
360  }
361 }
362 
363 bool
365 {
366  Real diff;
367 
368  for (unsigned int i = 0; i < _num_uo_state_vars; ++i)
369  {
370  unsigned int n = (*_mat_prop_state_vars[i])[_qp].size();
371  for (unsigned j = 0; j < n; j++)
372  {
373  diff = std::abs((*_mat_prop_state_vars[i])[_qp][j] -
374  _state_vars_prev[i][j]); // Calculate increment size
375  if (std::abs((*_mat_prop_state_vars_old[i])[_qp][j]) < _zero_tol && diff > _zero_tol)
376  return true;
377  if (std::abs((*_mat_prop_state_vars_old[i])[_qp][j]) > _zero_tol &&
378  diff > _stol * std::abs((*_mat_prop_state_vars_old[i])[_qp][j]))
379  return true;
380  }
381  }
382  return false;
383 }
384 
385 void
387 {
388  // TODO: remove this nasty const_cast if you can figure out how
389  for (unsigned int i = 0; i < _num_uo_state_vars; ++i)
390  const_cast<MaterialProperty<std::vector<Real>> &>(*_mat_prop_state_vars_old[i])[_qp] =
391  (*_mat_prop_state_vars[i])[_qp];
392 
394 }
395 
396 void
398 {
399 }
400 
401 void
403 {
404  unsigned int iter = 0;
405  RankTwoTensor dpk2;
406  Real rnorm, rnorm0, rnorm_prev;
407 
408  // Calculate stress residual
409  calcResidJacob();
410  if (_err_tol)
411  {
412 #ifdef DEBUG
413  mooseWarning("FiniteStrainUObasedCP: Slip increment exceeds tolerance - Element number ",
414  _current_elem->id(),
415  " Gauss point = ",
416  _qp);
417 #endif
418  return;
419  }
420 
421  rnorm = _resid.L2norm();
422  rnorm0 = rnorm;
423 
424  // Check for stress residual tolerance
425  while (rnorm > _rtol * rnorm0 && rnorm0 > _abs_tol && iter < _maxiter)
426  {
427  // Calculate stress increment
428  dpk2 = -_jac.invSymm() * _resid;
429  _pk2[_qp] = _pk2[_qp] + dpk2;
430  calcResidJacob();
431 
432  if (_err_tol)
433  {
434 #ifdef DEBUG
435  mooseWarning("FiniteStrainUObasedCP: Slip increment exceeds tolerance - Element number ",
436  _current_elem->id(),
437  " Gauss point = ",
438  _qp);
439 #endif
440  return;
441  }
442 
443  rnorm_prev = rnorm;
444  rnorm = _resid.L2norm();
445 
446  if (_use_line_search && rnorm > rnorm_prev && !lineSearchUpdate(rnorm_prev, dpk2))
447  {
448 #ifdef DEBUG
449  mooseWarning("FiniteStrainUObasedCP: Failed with line search");
450 #endif
451  _err_tol = true;
452  return;
453  }
454 
455  if (_use_line_search)
456  rnorm = _resid.L2norm();
457 
458  iter++;
459  }
460 
461  if (iter >= _maxiter)
462  {
463 #ifdef DEBUG
464  mooseWarning("FiniteStrainUObasedCP: Stress Integration error rmax = ", rnorm);
465 #endif
466  _err_tol = true;
467  }
468 }
469 
470 void
472 {
473  _fp[_qp] = _fp_inv.inverse();
474 }
475 
476 void
478 {
479  for (unsigned int i = 0; i < _num_uo_state_vars; ++i)
480  _state_vars_prev[i] = (*_mat_prop_state_vars[i])[_qp];
481 
482  for (unsigned int i = 0; i < _num_uo_state_var_evol_rate_comps; ++i)
483  _uo_state_var_evol_rate_comps[i]->calcStateVariableEvolutionRateComponent(
484  _qp, (*_mat_prop_state_var_evol_rate_comps[i])[_qp]);
485 
486  for (unsigned int i = 0; i < _num_uo_state_vars; ++i)
487  {
488  if (!_uo_state_vars[i]->updateStateVariable(_qp, _dt, (*_mat_prop_state_vars[i])[_qp]))
489  _err_tol = true;
490  }
491 
492  for (unsigned int i = 0; i < _num_uo_slip_resistances; ++i)
493  _uo_slip_resistances[i]->calcSlipResistance(_qp, (*_mat_prop_slip_resistances[i])[_qp]);
494 }
495 
496 // Calculates stress residual equation and jacobian
497 void
499 {
500  calcResidual();
501  if (_err_tol)
502  return;
503  calcJacobian();
504 }
505 
506 void
508 {
509  for (unsigned int i = 0; i < _num_uo_slip_rates; ++i)
510  {
511  if (!_uo_slip_rates[i]->calcSlipRate(_qp, _dt, (*_mat_prop_slip_rates[i])[_qp]))
512  {
513  _err_tol = true;
514  return;
515  }
516  }
517 }
518 
519 void
521 {
522  RankTwoTensor iden, ce, ee, ce_pk2, eqv_slip_incr, pk2_new;
523 
524  iden.zero();
525  iden.addIa(1.0);
526 
527  eqv_slip_incr.zero();
528 
529  getSlipRates();
530 
531  if (_err_tol)
532  return;
533 
534  for (unsigned int i = 0; i < _num_uo_slip_rates; ++i)
535  for (unsigned int j = 0; j < _uo_slip_rates[i]->variableSize(); ++j)
536  eqv_slip_incr += (*_flow_direction[i])[_qp][j] * (*_mat_prop_slip_rates[i])[_qp][j] * _dt;
537 
538  eqv_slip_incr = iden - eqv_slip_incr;
539  _fp_inv = _fp_old_inv * eqv_slip_incr;
540  _fe = _dfgrd_tmp * _fp_inv;
541 
542  ce = _fe.transpose() * _fe;
543  ee = ce - iden;
544  ee *= 0.5;
545 
546  pk2_new = _elasticity_tensor[_qp] * ee;
547 
548  _resid = _pk2[_qp] - pk2_new;
549 }
550 
551 void
553 {
554  RankFourTensor dfedfpinv, deedfe, dfpinvdpk2;
555 
556  for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
557  for (unsigned int j = 0; j < LIBMESH_DIM; ++j)
558  for (unsigned int k = 0; k < LIBMESH_DIM; ++k)
559  dfedfpinv(i, j, k, j) = _dfgrd_tmp(i, k);
560 
561  for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
562  for (unsigned int j = 0; j < LIBMESH_DIM; ++j)
563  for (unsigned int k = 0; k < LIBMESH_DIM; ++k)
564  {
565  deedfe(i, j, k, i) = deedfe(i, j, k, i) + _fe(k, j) * 0.5;
566  deedfe(i, j, k, j) = deedfe(i, j, k, j) + _fe(k, i) * 0.5;
567  }
568 
569  for (unsigned int i = 0; i < _num_uo_slip_rates; ++i)
570  {
571  unsigned int nss = _uo_slip_rates[i]->variableSize();
572  std::vector<RankTwoTensor> dtaudpk2(nss), dfpinvdslip(nss);
573  std::vector<Real> dslipdtau;
574  dslipdtau.resize(nss);
575  _uo_slip_rates[i]->calcSlipRateDerivative(_qp, _dt, dslipdtau);
576  for (unsigned int j = 0; j < nss; j++)
577  {
578  dtaudpk2[j] = (*_flow_direction[i])[_qp][j];
579  dfpinvdslip[j] = -_fp_old_inv * (*_flow_direction[i])[_qp][j];
580  }
581 
582  for (unsigned int j = 0; j < nss; j++)
583  dfpinvdpk2 += (dfpinvdslip[j] * dslipdtau[j] * _dt).outerProduct(dtaudpk2[j]);
584  }
585  _jac =
586  RankFourTensor::IdentityFour() - (_elasticity_tensor[_qp] * deedfe * dfedfpinv * dfpinvdpk2);
587 }
588 
589 void
591 {
592  switch (_tan_mod_type)
593  {
594  case 0:
596  break;
597  default:
599  }
600 }
601 
602 void
604 {
605  RankFourTensor tan_mod;
606  RankTwoTensor pk2fet, fepk2;
607  RankFourTensor deedfe, dsigdpk2dfe, dfedf;
608 
609  // Fill in the matrix stiffness material property
610  for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
611  for (unsigned int j = 0; j < LIBMESH_DIM; ++j)
612  for (unsigned int k = 0; k < LIBMESH_DIM; ++k)
613  {
614  deedfe(i, j, k, i) = deedfe(i, j, k, i) + _fe(k, j) * 0.5;
615  deedfe(i, j, k, j) = deedfe(i, j, k, j) + _fe(k, i) * 0.5;
616  }
617 
618  dsigdpk2dfe = _fe.mixedProductIkJl(_fe) * _elasticity_tensor[_qp] * deedfe;
619 
620  pk2fet = _pk2[_qp] * _fe.transpose();
621  fepk2 = _fe * _pk2[_qp];
622 
623  for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
624  for (unsigned int j = 0; j < LIBMESH_DIM; ++j)
625  for (unsigned int l = 0; l < LIBMESH_DIM; ++l)
626  {
627  tan_mod(i, j, i, l) += pk2fet(l, j);
628  tan_mod(i, j, j, l) += fepk2(i, l);
629  }
630 
631  tan_mod += dsigdpk2dfe;
632 
633  Real je = _fe.det();
634  if (je > 0.0)
635  tan_mod /= je;
636 
637  for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
638  for (unsigned int j = 0; j < LIBMESH_DIM; ++j)
639  for (unsigned int l = 0; l < LIBMESH_DIM; ++l)
640  dfedf(i, j, i, l) = _fp_inv(l, j);
641 
642  _Jacobian_mult[_qp] = tan_mod * dfedf;
643 }
644 
645 void
647 {
648  // update jacobian_mult
650 }
651 
652 bool
653 FiniteStrainUObasedCP::lineSearchUpdate(const Real rnorm_prev, const RankTwoTensor dpk2)
654 {
655  switch (_lsrch_method)
656  {
657  case 0: // CUT_HALF
658  {
659  Real rnorm;
660  Real step = 1.0;
661 
662  do
663  {
664  _pk2[_qp] = _pk2[_qp] - step * dpk2;
665  step /= 2.0;
666  _pk2[_qp] = _pk2[_qp] + step * dpk2;
667 
668  calcResidual();
669  rnorm = _resid.L2norm();
670  } while (rnorm > rnorm_prev && step > _min_lsrch_step);
671 
672  // has norm improved or is the step still above minumum search step size?
673  return (rnorm <= rnorm_prev || step > _min_lsrch_step);
674  }
675 
676  case 1: // BISECTION
677  {
678  unsigned int count = 0;
679  Real step_a = 0.0;
680  Real step_b = 1.0;
681  Real step = 1.0;
682  Real s_m = 1000.0;
683  Real rnorm = 1000.0;
684 
685  calcResidual();
686  Real s_b = _resid.doubleContraction(dpk2);
687  Real rnorm1 = _resid.L2norm();
688  _pk2[_qp] = _pk2[_qp] - dpk2;
689  calcResidual();
690  Real s_a = _resid.doubleContraction(dpk2);
691  Real rnorm0 = _resid.L2norm();
692  _pk2[_qp] = _pk2[_qp] + dpk2;
693 
694  if ((rnorm1 / rnorm0) < _lsrch_tol || s_a * s_b > 0)
695  {
696  calcResidual();
697  return true;
698  }
699 
700  while ((rnorm / rnorm0) > _lsrch_tol && count < _lsrch_max_iter)
701  {
702  _pk2[_qp] = _pk2[_qp] - step * dpk2;
703  step = 0.5 * (step_b + step_a);
704  _pk2[_qp] = _pk2[_qp] + step * dpk2;
705  calcResidual();
706  s_m = _resid.doubleContraction(dpk2);
707  rnorm = _resid.L2norm();
708 
709  if (s_m * s_a < 0.0)
710  {
711  step_b = step;
712  s_b = s_m;
713  }
714  if (s_m * s_b < 0.0)
715  {
716  step_a = step;
717  s_a = s_m;
718  }
719  count++;
720  }
721 
722  // below tolerance and max iterations?
723  return ((rnorm / rnorm0) < _lsrch_tol && count < _lsrch_max_iter);
724  }
725 
726  default:
727  mooseError("Line search method is not provided.");
728  }
729 }
MaterialProperty< RankFourTensor > & _Jacobian_mult
derivative of stress w.r.t. strain (_dstress_dstrain)
virtual bool isStateVariablesConverged()
evaluates convergence of state variables.
const MaterialProperty< RankTwoTensor > & _pk2_old
std::vector< std::vector< Real > > _state_vars_old
Local state variable.
ComputeStressBase is the base class for stress tensors.
virtual void postSolveStatevar()
update internal variable after solve.
bool _err_tol
Flag to check whether convergence is achieved.
std::vector< MaterialProperty< std::vector< Real > > * > _mat_prop_slip_rates
Slip rates material property.
Real _stol
Internal variable update equation tolerance.
std::vector< MaterialProperty< std::vector< Real > > * > _mat_prop_state_var_evol_rate_comps
State variable evolution rate component material property.
std::vector< const CrystalPlasticitySlipResistance * > _uo_slip_resistances
User objects that define the slip resistance.
MooseEnum _tan_mod_type
Type of tangent moduli calculation.
unsigned int _lsrch_max_iter
Line search bisection method maximum iteration number.
Real _zero_tol
Residual tolerance when variable value is zero. Default 1e-12.
FiniteStrainUObasedCP(const InputParameters &parameters)
std::vector< MaterialProperty< std::vector< Real > > * > _mat_prop_slip_resistances
Slip resistance material property.
std::vector< const CrystalPlasticityStateVarRateComponent * > _uo_state_var_evol_rate_comps
User objects that define the state variable evolution rate component.
virtual void initQpStatefulProperties()
initializes the stateful properties such as stress, plastic deformation gradient, slip system resista...
RankTwoTensor _delta_dfgrd
Used for substepping; Uniformly divides the increment in deformation gradient.
MaterialProperty< RankTwoTensor > & _stress
std::vector< std::vector< Real > > _state_vars_prev
Local old state variable.
unsigned int _maxiter
Maximum number of iterations for stress update.
std::vector< const CrystalPlasticityStateVariable * > _uo_state_vars
User objects that define the state variable.
virtual void calcTangentModuli()
calculate the tangent moduli for preconditioner.
virtual void calcResidJacob()
calls the residual and jacobian functions used in the stress update algorithm.
virtual void updateSlipSystemResistanceAndStateVariable()
updates the slip system resistances and state variables.
virtual void calcResidual()
calculate stress residual.
virtual void computeQpStress()
updates the stress at a quadrature point.
virtual void postSolveQp()
update stress and internal variable after solve.
MaterialProperty< RankTwoTensor > & _fp
virtual void getSlipRates()
updates the slip rates.
Real _lsrch_tol
Line search bisection method tolerance.
virtual void postSolveStress()
update stress and plastic deformation gradient after solve.
virtual void preSolveStress()
set variables for stress solve.
unsigned int _num_uo_slip_rates
Number of slip rate user objects.
unsigned int _maxiterg
Maximum number of iterations for internal variable update.
const MaterialProperty< RankTwoTensor > & _deformation_gradient
std::vector< const CrystalPlasticitySlipRate * > _uo_slip_rates
User objects that define the slip rate.
RankFourTensor _jac
Jacobian tensor.
unsigned int _num_uo_state_vars
Number of state variable user objects.
unsigned int _num_uo_slip_resistances
Number of slip resistance user objects.
MaterialProperty< RankTwoTensor > & _lag_e
unsigned int _num_uo_state_var_evol_rate_comps
Number of state variable evolution rate component user objects.
Real _rtol
Stress residual equation relative tolerance.
bool lineSearchUpdate(const Real rnorm_prev, const RankTwoTensor)
performs the line search update
virtual void preSolveQp()
set variables for stress and internal variable solve.
virtual void elastoPlasticTangentModuli()
calculate the exact tangent moduli for preconditioner.
const MaterialProperty< RankTwoTensor > & _crysrot
Crystal rotation.
virtual void elasticTangentModuli()
calculate the elastic tangent moduli for preconditioner.
InputParameters validParams< ComputeStressBase >()
std::vector< const MaterialProperty< std::vector< Real > > * > _mat_prop_state_vars_old
Old state variable material property.
virtual void preSolveStatevar()
set variables for internal variable solve.
MaterialProperty< RankTwoTensor > & _update_rot
virtual void solveQp()
solve stress and internal variables.
RankTwoTensor _resid
Residual tensor.
const MaterialProperty< RankFourTensor > & _elasticity_tensor
bool _use_line_search
Flag to activate line serach.
const MaterialProperty< RankTwoTensor > & _fp_old
unsigned int _max_substep_iter
Maximum number of substep iterations.
Real _abs_tol
Stress residual equation absolute tolerance.
virtual void solveStatevar()
solve internal variables.
const MaterialProperty< RankTwoTensor > & _deformation_gradient_old
std::vector< MaterialProperty< std::vector< Real > > * > _mat_prop_state_vars
State variable material property.
std::vector< MaterialProperty< std::vector< RankTwoTensor > > * > _flow_direction
Real _min_lsrch_step
Minimum line search step size.
virtual void solveStress()
solves for stress, updates plastic deformation gradient.
MooseEnum _lsrch_method
Line search method.
InputParameters validParams< FiniteStrainUObasedCP >()
MaterialProperty< RankTwoTensor > & _pk2
virtual void calcJacobian()
calculate jacobian.