www.mooseframework.org
MultiParameterPlasticityStressUpdate.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 
9 #include "Conversion.h" // for stringify
10 #include "libmesh/utility.h" // for Utility::pow
11 
12 template <>
13 InputParameters
15 {
16  InputParameters params = validParams<StressUpdateBase>();
17  params.addClassDescription("Return-map and Jacobian algorithms for plastic models where the "
18  "yield function and flow directions depend on multiple functions of "
19  "stress");
20  params.addParam<std::string>("base_name",
21  "Optional parameter that allows the user to define "
22  "multiple plastic models on the same block, and the "
23  "plastic_internal_parameter, plastic_yield_function, "
24  "plastic_NR_iterations and plastic_linesearch_needed Material "
25  "Properties will be prepended by this string");
26  params.addRangeCheckedParam<unsigned int>(
27  "max_NR_iterations",
28  20,
29  "max_NR_iterations>0",
30  "Maximum number of Newton-Raphson iterations allowed during the return-map algorithm");
31  params.addParam<bool>("perform_finite_strain_rotations",
32  false,
33  "Tensors are correctly rotated "
34  "in finite-strain simulations. "
35  "For optimal performance you can "
36  "set this to 'false' if you are "
37  "only ever using small strains");
38  params.addRequiredParam<Real>(
39  "smoothing_tol",
40  "Intersections of the yield surfaces will be smoothed by this amount (this "
41  "is measured in units of stress). Often this is related to other physical "
42  "parameters (eg, 0.1*cohesion) but it is important to set this small enough "
43  "so that the individual yield surfaces do not mix together in the smoothing "
44  "process to produce a result where no stress is admissible (for example, "
45  "mixing together tensile and compressive failure envelopes).");
46  params.addRequiredParam<Real>("yield_function_tol",
47  "The return-map process will be deemed to have converged if all "
48  "yield functions are within yield_function_tol of zero. If this "
49  "is set very low then precision-loss might be encountered: if the "
50  "code detects precision loss then it also deems the return-map "
51  "process has converged.");
52  MooseEnum tangent_operator("elastic nonlinear", "nonlinear");
53  params.addParam<Real>("min_step_size",
54  1.0,
55  "In order to help the Newton-Raphson procedure, the applied strain "
56  "increment may be applied in sub-increments of size greater than this "
57  "value. Usually it is better for Moose's nonlinear convergence to "
58  "increase max_NR_iterations rather than decrease this parameter.");
59  params.addParam<bool>("warn_about_precision_loss",
60  false,
61  "Output a message to the console every "
62  "time precision-loss is encountered "
63  "during the Newton-Raphson process");
64  params.addParam<std::vector<Real>>("admissible_stress",
65  "A single admissible value of the value of the stress "
66  "parameters for internal parameters = 0. This is used "
67  "to initialize the return-mapping algorithm during the first "
68  "nonlinear iteration. If not given then it is assumed that "
69  "stress parameters = 0 is admissible.");
70  return params;
71 }
72 
74  const InputParameters & parameters, unsigned num_sp, unsigned num_yf, unsigned num_intnl)
75  : StressUpdateBase(parameters),
76  _num_sp(num_sp),
77  _definitely_ok_sp(isParamValid("admissible_stress")
78  ? getParam<std::vector<Real>>("admissible_stress")
79  : std::vector<Real>(_num_sp, 0.0)),
80  _Eij(num_sp, std::vector<Real>(num_sp)),
81  _En(1.0),
82  _Cij(num_sp, std::vector<Real>(num_sp)),
83  _num_yf(num_yf),
84  _num_intnl(num_intnl),
85  _base_name(isParamValid("base_name") ? getParam<std::string>("base_name") + "_" : ""),
86  _max_nr_its(getParam<unsigned>("max_NR_iterations")),
87  _perform_finite_strain_rotations(getParam<bool>("perform_finite_strain_rotations")),
88  _smoothing_tol(getParam<Real>("smoothing_tol")),
89  _f_tol(getParam<Real>("yield_function_tol")),
90  _f_tol2(Utility::pow<2>(getParam<Real>("yield_function_tol"))),
91  _min_step_size(getParam<Real>("min_step_size")),
92  _step_one(declareRestartableData<bool>("step_one", true)),
93  _warn_about_precision_loss(getParam<bool>("warn_about_precision_loss")),
94 
95  _plastic_strain(declareProperty<RankTwoTensor>(_base_name + "plastic_strain")),
96  _plastic_strain_old(getMaterialPropertyOld<RankTwoTensor>(_base_name + "plastic_strain")),
97  _intnl(declareProperty<std::vector<Real>>(_base_name + "plastic_internal_parameter")),
98  _intnl_old(
99  getMaterialPropertyOld<std::vector<Real>>(_base_name + "plastic_internal_parameter")),
100  _yf(declareProperty<std::vector<Real>>(_base_name + "plastic_yield_function")),
101  _iter(declareProperty<Real>(_base_name +
102  "plastic_NR_iterations")), // this is really an unsigned int, but
103  // for visualisation i convert it to Real
104  _linesearch_needed(
105  declareProperty<Real>(_base_name + "plastic_linesearch_needed")), // this is really a
106  // boolean, but for
107  // visualisation i
108  // convert it to Real
109  _trial_sp(num_sp),
110  _stress_trial(RankTwoTensor()),
111  _rhs(num_sp + 1),
112  _dvar_dtrial(num_sp + 1, std::vector<Real>(num_sp, 0.0)),
113  _ok_sp(num_sp),
114  _ok_intnl(num_intnl),
115  _del_stress_params(num_sp),
116  _current_sp(num_sp),
117  _current_intnl(num_intnl)
118 {
119  if (_definitely_ok_sp.size() != _num_sp)
120  mooseError("MultiParameterPlasticityStressUpdate: admissible_stress parameter must consist of ",
121  _num_sp,
122  " numbers");
123 }
124 
125 void
127 {
128  _plastic_strain[_qp].zero();
129  _intnl[_qp].assign(_num_intnl, 0);
130  _yf[_qp].assign(_num_yf, 0);
131  _iter[_qp] = 0.0;
132  _linesearch_needed[_qp] = 0.0;
133 }
134 
135 void
137 {
139  for (unsigned i = 0; i < _num_intnl; ++i)
140  _intnl[_qp][i] = _intnl_old[_qp][i];
141 }
142 
143 void
144 MultiParameterPlasticityStressUpdate::updateState(RankTwoTensor & strain_increment,
145  RankTwoTensor & inelastic_strain_increment,
146  const RankTwoTensor & rotation_increment,
147  RankTwoTensor & stress_new,
148  const RankTwoTensor & stress_old,
149  const RankFourTensor & elasticity_tensor,
150  const RankTwoTensor & /*elastic_strain_old*/,
151  bool compute_full_tangent_operator,
152  RankFourTensor & tangent_operator)
153 {
155 
156  if (_t_step >= 2)
157  _step_one = false;
158 
159  // initially assume an elastic deformation
160  std::copy(_intnl_old[_qp].begin(), _intnl_old[_qp].end(), _intnl[_qp].begin());
161  _iter[_qp] = 0.0;
162  _linesearch_needed[_qp] = 0.0;
163 
164  computeStressParams(stress_new, _trial_sp);
166 
167  if (yieldF(_yf[_qp]) <= _f_tol)
168  {
170  inelastic_strain_increment.zero();
171  if (_fe_problem.currentlyComputingJacobian())
172  tangent_operator = elasticity_tensor;
173  return;
174  }
175 
176  _stress_trial = stress_new;
177  /* The trial stress must be inadmissible
178  * so we need to return to the yield surface. The following
179  * equations must be satisfied.
180  *
181  * 0 = rhs[0] = S[0] - S[0]^trial + ga * E[0, i] * dg/dS[i]
182  * 0 = rhs[1] = S[1] - S[1]^trial + ga * E[1, i] * dg/dS[i]
183  * ...
184  * 0 = rhs[N-1] = S[N-1] - S[N-1]^trial + ga * E[N-1, i] * dg/dS[i]
185  * 0 = rhs[N] = f(S, intnl)
186  *
187  * as well as equations defining intnl parameters as functions of
188  * stress_params, trial_stress_params and intnl_old
189  *
190  * The unknowns are S[0], ..., S[N-1], gaE, and the intnl parameters.
191  * Here gaE = ga * _En (the _En serves to make gaE similar magnitude to S)
192  * I find it convenient to solve the first N+1 equations for p, q and gaE,
193  * while substituting the "intnl parameters" equations into these during the solve process
194  */
195 
196  for (auto & deriv : _dvar_dtrial)
197  deriv.assign(_num_sp, 0.0);
198 
199  preReturnMapV(_trial_sp, stress_new, _intnl_old[_qp], _yf[_qp], elasticity_tensor);
200 
201  setEffectiveElasticity(elasticity_tensor);
202 
203  if (_step_one)
204  std::copy(_definitely_ok_sp.begin(), _definitely_ok_sp.end(), _ok_sp.begin());
205  else
206  computeStressParams(stress_old, _ok_sp);
207  std::copy(_intnl_old[_qp].begin(), _intnl_old[_qp].end(), _ok_intnl.begin());
208 
209  // Return-map problem: must apply the following changes in stress_params,
210  // and find the returned stress_params and gaE
211  for (unsigned i = 0; i < _num_sp; ++i)
212  _del_stress_params[i] = _trial_sp[i] - _ok_sp[i];
213 
214  Real step_taken = 0.0; // amount of del_stress_params that we've applied and the return-map
215  // problem has succeeded
216  Real step_size = 1.0; // potentially can apply del_stress_params in substeps
217  Real gaE_total = 0.0;
218 
219  // current values of the yield function, derivatives, etc
220  yieldAndFlow smoothed_q;
221 
222  // In the following sub-stepping procedure it is possible that
223  // the last step is an elastic step, and therefore smoothed_q won't
224  // be computed on the last step, so we have to compute it.
225  bool smoothed_q_calculated = false;
226 
227  while (step_taken < 1.0 && step_size >= _min_step_size)
228  {
229  if (1.0 - step_taken < step_size)
230  // prevent over-shoots of substepping
231  step_size = 1.0 - step_taken;
232 
233  // trial variables in terms of admissible variables
234  for (unsigned i = 0; i < _num_sp; ++i)
235  _trial_sp[i] = _ok_sp[i] + step_size * _del_stress_params[i];
236 
237  // initialize variables that are to be found via Newton-Raphson
239  Real gaE = 0.0;
240 
241  // flags indicating failure of Newton-Raphson and line-search
242  int nr_failure = 0;
243  int ls_failure = 0;
244 
245  // NR iterations taken in this substep
246  unsigned step_iter = 0;
247 
248  // The residual-squared for the line-search
249  Real res2 = 0.0;
250 
251  if (step_size < 1.0 && yieldF(_trial_sp, _ok_intnl) <= _f_tol)
252  // This is an elastic step
253  // The "step_size < 1.0" in above condition is for efficiency: we definitely
254  // know that this is a plastic step if step_size = 1.0
255  smoothed_q_calculated = false;
256  else
257  {
258  // this is a plastic step
259 
260  // initialize current_sp, gaE and current_intnl based on the non-smoothed situation
262  // and find the smoothed yield function, flow potential and derivatives
264  smoothed_q_calculated = true;
265  calculateRHS(_trial_sp, _current_sp, gaE, smoothed_q, _rhs);
266  res2 = calculateRes2(_rhs);
267 
268  // Perform a Newton-Raphson with linesearch to get current_sp, gaE, and also smoothed_q
269  while (res2 > _f_tol2 && step_iter < _max_nr_its && nr_failure == 0 && ls_failure == 0)
270  {
271  // solve the linear system and store the answer (the "updates") in rhs
272  nr_failure = nrStep(smoothed_q, _trial_sp, _current_sp, _current_intnl, gaE, _rhs);
273  if (nr_failure != 0)
274  break;
275 
276  // handle precision loss
277  if (precisionLoss(_rhs, _current_sp, gaE))
278  {
280  {
281  Moose::err << "MultiParameterPlasticityStressUpdate: precision-loss in element "
282  << _current_elem->id() << " quadpoint=" << _qp << ". Stress_params =";
283  for (unsigned i = 0; i < _num_sp; ++i)
284  Moose::err << " " << _current_sp[i];
285  Moose::err << " gaE = " << gaE << "\n";
286  }
287  res2 = 0.0;
288  break;
289  }
290 
291  // apply (parts of) the updates, re-calculate smoothed_q, and res2
292  ls_failure = lineSearch(res2,
293  _current_sp,
294  gaE,
295  _trial_sp,
296  smoothed_q,
297  _ok_intnl,
299  _rhs,
300  _linesearch_needed[_qp]);
301  step_iter++;
302  }
303  }
304 
305  if (res2 <= _f_tol2 && step_iter < _max_nr_its && nr_failure == 0 && ls_failure == 0 &&
306  gaE >= 0.0)
307  {
308  // this Newton-Raphson worked fine, or this was an elastic step
309  std::copy(_current_sp.begin(), _current_sp.end(), _ok_sp.begin());
310  gaE_total += gaE;
311  step_taken += step_size;
313  std::copy(_intnl[_qp].begin(), _intnl[_qp].end(), _ok_intnl.begin());
314  // calculate dp/dp_trial, dp/dq_trial, etc, for Jacobian
315  dVardTrial(!smoothed_q_calculated,
316  _trial_sp,
317  _ok_sp,
318  gaE,
319  _ok_intnl,
320  smoothed_q,
321  step_size,
322  compute_full_tangent_operator,
323  _dvar_dtrial);
324  if (static_cast<Real>(step_iter) > _iter[_qp])
325  _iter[_qp] = static_cast<Real>(step_iter);
326  step_size *= 1.1;
327  }
328  else
329  {
330  // Newton-Raphson + line-search process failed
331  step_size *= 0.5;
332  }
333  }
334 
335  if (step_size < _min_step_size)
336  errorHandler("MultiParameterPlasticityStressUpdate: Minimum step-size violated");
337 
338  // success!
339  finalizeReturnProcess(rotation_increment);
340  yieldFunctionValuesV(_ok_sp, _intnl[_qp], _yf[_qp]);
341 
342  if (!smoothed_q_calculated)
343  smoothed_q = smoothAllQuantities(_ok_sp, _intnl[_qp]);
344 
346  _stress_trial, _ok_sp, gaE_total, _intnl[_qp], smoothed_q, elasticity_tensor, stress_new);
347 
349  gaE_total,
350  smoothed_q,
351  elasticity_tensor,
352  stress_new,
353  inelastic_strain_increment);
354 
355  strain_increment = strain_increment - inelastic_strain_increment;
356  _plastic_strain[_qp] = _plastic_strain_old[_qp] + inelastic_strain_increment;
357 
358  if (_fe_problem.currentlyComputingJacobian())
359  // for efficiency, do not compute the tangent operator if not currently computing Jacobian
361  _trial_sp,
362  stress_new,
363  _ok_sp,
364  gaE_total,
365  smoothed_q,
366  elasticity_tensor,
367  compute_full_tangent_operator,
368  _dvar_dtrial,
369  tangent_operator);
370 }
371 
373 MultiParameterPlasticityStressUpdate::smoothAllQuantities(const std::vector<Real> & stress_params,
374  const std::vector<Real> & intnl) const
375 {
376  std::vector<yieldAndFlow> all_q(_num_yf, yieldAndFlow(_num_sp, _num_intnl));
377  computeAllQV(stress_params, intnl, all_q);
378 
379  /* This routine holds the key to my smoothing strategy. It
380  * may be proved that this smoothing strategy produces a
381  * yield surface that is both C2 differentiable and convex,
382  * assuming the individual yield functions are C2 and
383  * convex too.
384  * Of course all the derivatives must also be smoothed.
385  * Also, I assume that d(flow potential)/dstress gets smoothed
386  * by the Yield Function (which produces a C2 flow potential).
387  * See the line identified in the loop below.
388  * Only time will tell whether this is a good strategy, but it
389  * works well in all tests so far. Convexity is irrelevant
390  * for the non-associated case, but at least the return-map
391  * problem should always have a unique solution.
392  * For two yield functions+flows, labelled 1 and 2, we
393  * should have
394  * d(g1 - g2) . d(f1 - f2) >= 0
395  * If not then the return-map problem for even the
396  * multi-surface plasticity with no smoothing won't have a
397  * unique solution. If the multi-surface plasticity has
398  * a unique solution then the smoothed version defined
399  * below will too.
400  */
401 
402  // res_f is the index that contains the smoothed yieldAndFlow
403  std::size_t res_f = 0;
404 
405  for (std::size_t a = 1; a < all_q.size(); ++a)
406  {
407  if (all_q[res_f].f >= all_q[a].f + _smoothing_tol)
408  // no smoothing is needed: res_f is already indexes the largest yield function
409  continue;
410  else if (all_q[a].f >= all_q[res_f].f + _smoothing_tol)
411  {
412  // no smoothing is needed, and res_f needs to index to all_q[a]
413  res_f = a;
414  continue;
415  }
416  else
417  {
418  // smoothing is required
419  const Real f_diff = all_q[res_f].f - all_q[a].f;
420  const Real ism = ismoother(f_diff);
421  const Real sm = smoother(f_diff);
422  const Real dsm = dsmoother(f_diff);
423  // we want: all_q[res_f].f = 0.5 * all_q[res_f].f + all_q[a].f + _smoothing_tol) + ism,
424  // but we have to do the derivatives first
425  for (unsigned i = 0; i < _num_sp; ++i)
426  {
427  for (unsigned j = 0; j < _num_sp; ++j)
428  all_q[res_f].d2g[i][j] =
429  0.5 * (all_q[res_f].d2g[i][j] + all_q[a].d2g[i][j]) +
430  dsm * (all_q[res_f].df[j] - all_q[a].df[j]) * (all_q[res_f].dg[i] - all_q[a].dg[i]) +
431  sm * (all_q[res_f].d2g[i][j] - all_q[a].d2g[i][j]);
432  for (unsigned j = 0; j < _num_intnl; ++j)
433  all_q[res_f].d2g_di[i][j] = 0.5 * (all_q[res_f].d2g_di[i][j] + all_q[a].d2g_di[i][j]) +
434  dsm * (all_q[res_f].df_di[j] - all_q[a].df_di[j]) *
435  (all_q[res_f].dg[i] - all_q[a].dg[i]) +
436  sm * (all_q[res_f].d2g_di[i][j] - all_q[a].d2g_di[i][j]);
437  }
438  for (unsigned i = 0; i < _num_sp; ++i)
439  {
440  all_q[res_f].df[i] = 0.5 * (all_q[res_f].df[i] + all_q[a].df[i]) +
441  sm * (all_q[res_f].df[i] - all_q[a].df[i]);
442  // whether the following (smoothing g with f's smoother) is a good strategy remains to be
443  // seen...
444  all_q[res_f].dg[i] = 0.5 * (all_q[res_f].dg[i] + all_q[a].dg[i]) +
445  sm * (all_q[res_f].dg[i] - all_q[a].dg[i]);
446  }
447  for (unsigned i = 0; i < _num_intnl; ++i)
448  all_q[res_f].df_di[i] = 0.5 * (all_q[res_f].df_di[i] + all_q[a].df_di[i]) +
449  sm * (all_q[res_f].df_di[i] - all_q[a].df_di[i]);
450  all_q[res_f].f = 0.5 * (all_q[res_f].f + all_q[a].f + _smoothing_tol) + ism;
451  }
452  }
453  return all_q[res_f];
454 }
455 
456 Real
458 {
459  if (std::abs(f_diff) >= _smoothing_tol)
460  return 0.0;
461  return -_smoothing_tol / M_PI * std::cos(0.5 * M_PI * f_diff / _smoothing_tol);
462 }
463 
464 Real
466 {
467  if (std::abs(f_diff) >= _smoothing_tol)
468  return 0.0;
469  return 0.5 * std::sin(f_diff * M_PI * 0.5 / _smoothing_tol);
470 }
471 
472 Real
474 {
475  if (std::abs(f_diff) >= _smoothing_tol)
476  return 0.0;
477  return 0.25 * M_PI / _smoothing_tol * std::cos(f_diff * M_PI * 0.5 / _smoothing_tol);
478 }
479 
480 int
482  std::vector<Real> & stress_params,
483  Real & gaE,
484  const std::vector<Real> & trial_stress_params,
485  yieldAndFlow & smoothed_q,
486  const std::vector<Real> & intnl_ok,
487  std::vector<Real> & intnl,
488  std::vector<Real> & rhs,
489  Real & linesearch_needed) const
490 {
491  const Real res2_old = res2;
492  const std::vector<Real> sp_params_old = stress_params;
493  const Real gaE_old = gaE;
494  const std::vector<Real> delta_nr_params = rhs;
495 
496  Real lam = 1.0; // line-search parameter
497  const Real lam_min = 1E-10; // minimum value of lam allowed
498  const Real slope = -2.0 * res2_old; // "Numerical Recipes" uses -b*A*x, in order to check for
499  // roundoff, but i hope the nrStep would warn if there were
500  // problems
501  Real tmp_lam; // cached value of lam used in quadratic & cubic line search
502  Real f2 = res2_old; // cached value of f = residual2 used in the cubic in the line search
503  Real lam2 = lam; // cached value of lam used in the cubic in the line search
504 
505  while (true)
506  {
507  // update variables using the current line-search parameter
508  for (unsigned i = 0; i < _num_sp; ++i)
509  stress_params[i] = sp_params_old[i] + lam * delta_nr_params[i];
510  gaE = gaE_old + lam * delta_nr_params[_num_sp];
511 
512  // and internal parameters
513  setIntnlValuesV(trial_stress_params, stress_params, intnl_ok, intnl);
514 
515  smoothed_q = smoothAllQuantities(stress_params, intnl);
516 
517  // update rhs for next-time through
518  calculateRHS(trial_stress_params, stress_params, gaE, smoothed_q, rhs);
519  res2 = calculateRes2(rhs);
520 
521  // do the line-search
522  if (res2 < res2_old + 1E-4 * lam * slope)
523  break;
524  else if (lam < lam_min)
525  return 1;
526  else if (lam == 1.0)
527  {
528  // model as a quadratic
529  tmp_lam = -0.5 * slope / (res2 - res2_old - slope);
530  }
531  else
532  {
533  // model as a cubic
534  const Real rhs1 = res2 - res2_old - lam * slope;
535  const Real rhs2 = f2 - res2_old - lam2 * slope;
536  const Real a = (rhs1 / Utility::pow<2>(lam) - rhs2 / Utility::pow<2>(lam2)) / (lam - lam2);
537  const Real b =
538  (-lam2 * rhs1 / Utility::pow<2>(lam) + lam * rhs2 / Utility::pow<2>(lam2)) / (lam - lam2);
539  if (a == 0.0)
540  tmp_lam = -slope / (2.0 * b);
541  else
542  {
543  const Real disc = Utility::pow<2>(b) - 3.0 * a * slope;
544  if (disc < 0)
545  tmp_lam = 0.5 * lam;
546  else if (b <= 0)
547  tmp_lam = (-b + std::sqrt(disc)) / (3.0 * a);
548  else
549  tmp_lam = -slope / (b + std::sqrt(disc));
550  }
551  if (tmp_lam > 0.5 * lam)
552  tmp_lam = 0.5 * lam;
553  }
554  lam2 = lam;
555  f2 = res2;
556  lam = std::max(tmp_lam, 0.1 * lam);
557  }
558 
559  if (lam < 1.0)
560  linesearch_needed = 1.0;
561  return 0;
562 }
563 
564 int
566  const std::vector<Real> & trial_stress_params,
567  const std::vector<Real> & stress_params,
568  const std::vector<Real> & intnl,
569  Real gaE,
570  std::vector<Real> & rhs) const
571 {
572  std::vector<std::vector<Real>> dintnl(_num_intnl, std::vector<Real>(_num_sp));
573  setIntnlDerivativesV(trial_stress_params, stress_params, intnl, dintnl);
574 
575  std::vector<double> jac((_num_sp + 1) * (_num_sp + 1));
576  dnRHSdVar(smoothed_q, dintnl, stress_params, gaE, jac);
577 
578  // use LAPACK to solve the linear system
579  const int nrhs = 1;
580  std::vector<int> ipiv(_num_sp + 1);
581  int info;
582  const int gesv_num_rhs = _num_sp + 1;
583  LAPACKgesv_(
584  &gesv_num_rhs, &nrhs, &jac[0], &gesv_num_rhs, &ipiv[0], &rhs[0], &gesv_num_rhs, &info);
585  return info;
586 }
587 
588 void
589 MultiParameterPlasticityStressUpdate::errorHandler(const std::string & message) const
590 {
591  throw MooseException(message);
592 }
593 
594 void
596 {
597 }
598 
599 void
601  const RankTwoTensor & /*rotation_increment*/)
602 {
603 }
604 
605 void
607  const std::vector<Real> & /*trial_stress_params*/,
608  const RankTwoTensor & /*stress_trial*/,
609  const std::vector<Real> & /*intnl_old*/,
610  const std::vector<Real> & /*yf*/,
611  const RankFourTensor & /*Eijkl*/)
612 {
613 }
614 
615 Real
616 MultiParameterPlasticityStressUpdate::yieldF(const std::vector<Real> & stress_params,
617  const std::vector<Real> & intnl) const
618 {
619  std::vector<Real> yfs(_num_yf);
620  yieldFunctionValuesV(stress_params, intnl, yfs);
621  return yieldF(yfs);
622 }
623 
624 Real
625 MultiParameterPlasticityStressUpdate::yieldF(const std::vector<Real> & yfs) const
626 {
627  Real yf = yfs[0];
628  for (std::size_t i = 1; i < yfs.size(); ++i)
629  if (yf >= yfs[i] + _smoothing_tol)
630  // no smoothing is needed, and yf is the biggest yield function
631  continue;
632  else if (yfs[i] >= yf + _smoothing_tol)
633  // no smoothing is needed, and yfs[i] is the biggest yield function
634  yf = yfs[i];
635  else
636  yf = 0.5 * (yf + yfs[i] + _smoothing_tol) + ismoother(yf - yfs[i]);
637  return yf;
638 }
639 
640 void
641 MultiParameterPlasticityStressUpdate::initializeVarsV(const std::vector<Real> & trial_stress_params,
642  const std::vector<Real> & intnl_old,
643  std::vector<Real> & stress_params,
644  Real & gaE,
645  std::vector<Real> & intnl) const
646 {
647  gaE = 0.0;
648  std::copy(trial_stress_params.begin(), trial_stress_params.end(), stress_params.begin());
649  std::copy(intnl_old.begin(), intnl_old.end(), intnl.begin());
650 }
651 
652 void
654  const RankTwoTensor & stress_trial,
655  const std::vector<Real> & /*trial_stress_params*/,
656  const RankTwoTensor & stress,
657  const std::vector<Real> & /*stress_params*/,
658  Real gaE,
659  const yieldAndFlow & smoothed_q,
660  const RankFourTensor & elasticity_tensor,
661  bool compute_full_tangent_operator,
662  const std::vector<std::vector<Real>> & dvar_dtrial,
663  RankFourTensor & cto)
664 {
665  cto = elasticity_tensor;
666  if (!compute_full_tangent_operator)
667  return;
668 
669  const Real ga = gaE / _En;
670 
671  const std::vector<RankTwoTensor> dsp = dstress_param_dstress(stress);
672  const std::vector<RankTwoTensor> dsp_trial = dstress_param_dstress(stress_trial);
673 
674  for (unsigned a = 0; a < _num_sp; ++a)
675  {
676  for (unsigned b = 0; b < _num_sp; ++b)
677  {
678  const RankTwoTensor t = elasticity_tensor * dsp_trial[a];
679  RankTwoTensor s = _Cij[b][a] * dsp[b];
680  for (unsigned c = 0; c < _num_sp; ++c)
681  s -= dsp[b] * _Cij[b][c] * dvar_dtrial[c][a];
682  s = elasticity_tensor * s;
683  cto -= s.outerProduct(t);
684  }
685  }
686 
687  const std::vector<RankFourTensor> d2sp = d2stress_param_dstress(stress);
688  RankFourTensor Tijab;
689  for (unsigned i = 0; i < _num_sp; ++i)
690  Tijab += smoothed_q.dg[i] * d2sp[i];
691  Tijab = ga * elasticity_tensor * Tijab;
692 
693  RankFourTensor inv = RankFourTensor(RankFourTensor::initIdentityFour) + Tijab;
694  try
695  {
696  inv = inv.transposeMajor().invSymm();
697  }
698  catch (const MooseException & e)
699  {
700  // Cannot form the inverse, so probably at some degenerate place in stress space.
701  // Just return with the "best estimate" of the cto.
702  mooseWarning("MultiParameterPlasticityStressUpdate: Cannot invert 1+T in consistent tangent "
703  "operator computation at quadpoint ",
704  _qp,
705  " of element ",
706  _current_elem->id());
707  return;
708  }
709 
710  cto = (cto.transposeMajor() * inv).transposeMajor();
711 }
712 
713 void
715  const RankTwoTensor & /*stress_trial*/,
716  Real gaE,
717  const yieldAndFlow & smoothed_q,
718  const RankFourTensor & /*elasticity_tensor*/,
719  const RankTwoTensor & returned_stress,
720  RankTwoTensor & inelastic_strain_increment) const
721 {
722  const std::vector<RankTwoTensor> dsp_dstress = dstress_param_dstress(returned_stress);
723  inelastic_strain_increment = RankTwoTensor();
724  for (unsigned i = 0; i < _num_sp; ++i)
725  inelastic_strain_increment += smoothed_q.dg[i] * dsp_dstress[i];
726  inelastic_strain_increment *= gaE / _En;
727 }
728 
729 Real
730 MultiParameterPlasticityStressUpdate::calculateRes2(const std::vector<Real> & rhs) const
731 {
732  Real res2 = 0.0;
733  for (const auto & r : rhs)
734  res2 += r * r;
735  return res2;
736 }
737 
738 void
739 MultiParameterPlasticityStressUpdate::calculateRHS(const std::vector<Real> & trial_stress_params,
740  const std::vector<Real> & stress_params,
741  Real gaE,
742  const yieldAndFlow & smoothed_q,
743  std::vector<Real> & rhs) const
744 {
745  const Real ga = gaE / _En;
746  for (unsigned i = 0; i < _num_sp; ++i)
747  {
748  rhs[i] = stress_params[i] - trial_stress_params[i];
749  for (unsigned j = 0; j < _num_sp; ++j)
750  rhs[i] += ga * _Eij[i][j] * smoothed_q.dg[j];
751  }
752  rhs[_num_sp] = smoothed_q.f;
753 }
754 
755 void
757  const std::vector<std::vector<Real>> & dintnl,
758  const std::vector<Real> & /*stress_params*/,
759  Real gaE,
760  std::vector<double> & jac) const
761 {
762  for (auto & jac_entry : jac)
763  jac_entry = 0.0;
764 
765  const Real ga = gaE / _En;
766 
767  unsigned ind = 0;
768  for (unsigned var = 0; var < _num_sp; ++var)
769  {
770  for (unsigned rhs = 0; rhs < _num_sp; ++rhs)
771  {
772  if (var == rhs)
773  jac[ind] -= 1.0;
774  for (unsigned j = 0; j < _num_sp; ++j)
775  {
776  jac[ind] -= ga * _Eij[rhs][j] * smoothed_q.d2g[j][var];
777  for (unsigned k = 0; k < _num_intnl; ++k)
778  jac[ind] -= ga * _Eij[rhs][j] * smoothed_q.d2g_di[j][k] * dintnl[k][var];
779  }
780  ind++;
781  }
782  // now rhs = _num_sp (that is, the yield function)
783  jac[ind] -= smoothed_q.df[var];
784  for (unsigned k = 0; k < _num_intnl; ++k)
785  jac[ind] -= smoothed_q.df_di[k] * dintnl[k][var];
786  ind++;
787  }
788 
789  // now var = _num_sp (that is, gaE)
790  for (unsigned rhs = 0; rhs < _num_sp; ++rhs)
791  {
792  for (unsigned j = 0; j < _num_sp; ++j)
793  jac[ind] -= (1.0 / _En) * _Eij[rhs][j] * smoothed_q.dg[j];
794  ind++;
795  }
796  // now rhs = _num_sp (that is, the yield function)
797  jac[ind] = 0.0;
798 }
799 
800 void
802  const std::vector<Real> & trial_stress_params,
803  const std::vector<Real> & stress_params,
804  Real gaE,
805  const std::vector<Real> & intnl,
806  const yieldAndFlow & smoothed_q,
807  Real step_size,
808  bool compute_full_tangent_operator,
809  std::vector<std::vector<Real>> & dvar_dtrial) const
810 {
811  if (!_fe_problem.currentlyComputingJacobian())
812  return;
813 
814  if (!compute_full_tangent_operator)
815  return;
816 
817  if (elastic_only)
818  {
819  // no change to gaE, and all off-diag stuff remains unchanged from previous step
820  for (unsigned v = 0; v < _num_sp; ++v)
821  dvar_dtrial[v][v] += step_size;
822  return;
823  }
824 
825  const Real ga = gaE / _En;
826 
827  std::vector<std::vector<Real>> dintnl(_num_intnl, std::vector<Real>(_num_sp));
828  setIntnlDerivativesV(trial_stress_params, stress_params, intnl, dintnl);
829 
830  // rhs is described elsewhere, the following are changes in rhs wrt the trial_stress_param
831  // values
832  // In the following we use d(intnl)/d(trial variable) = - d(intnl)/d(variable)
833  std::vector<Real> rhs_cto((_num_sp + 1) * _num_sp);
834 
835  unsigned ind = 0;
836  for (unsigned a = 0; a < _num_sp; ++a)
837  {
838  // change in RHS[b] wrt changes in stress_param_trial[a]
839  for (unsigned b = 0; b < _num_sp; ++b)
840  {
841  if (a == b)
842  rhs_cto[ind] -= 1.0;
843  for (unsigned j = 0; j < _num_sp; ++j)
844  for (unsigned k = 0; k < _num_intnl; ++k)
845  rhs_cto[ind] -= ga * _Eij[b][j] * smoothed_q.d2g_di[j][k] * dintnl[k][a];
846  ind++;
847  }
848  // now b = _num_sp (that is, the yield function)
849  for (unsigned k = 0; k < _num_intnl; ++k)
850  rhs_cto[ind] -= smoothed_q.df_di[k] * dintnl[k][a];
851  ind++;
852  }
853 
854  // jac = d(-rhs)/d(var)
855  std::vector<double> jac((_num_sp + 1) * (_num_sp + 1));
856  dnRHSdVar(smoothed_q, dintnl, stress_params, gaE, jac);
857 
858  std::vector<int> ipiv(_num_sp + 1);
859  int info;
860  const int gesv_num_rhs = _num_sp + 1;
861  const int gesv_num_pq = _num_sp;
862  LAPACKgesv_(&gesv_num_rhs,
863  &gesv_num_pq,
864  &jac[0],
865  &gesv_num_rhs,
866  &ipiv[0],
867  &rhs_cto[0],
868  &gesv_num_rhs,
869  &info);
870  if (info != 0)
871  errorHandler("MultiParameterPlasticityStressUpdate: PETSC LAPACK gsev routine returned with "
872  "error code " +
873  Moose::stringify(info));
874 
875  ind = 0;
876  std::vector<std::vector<Real>> dvarn_dtrialn(_num_sp + 1, std::vector<Real>(_num_sp, 0.0));
877  for (unsigned spt = 0; spt < _num_sp; ++spt) // loop over trial stress-param variables
878  {
879  for (unsigned v = 0; v < _num_sp; ++v) // loop over variables in NR procedure
880  {
881  dvarn_dtrialn[v][spt] = rhs_cto[ind];
882  ind++;
883  }
884  // the final NR variable is gaE
885  dvarn_dtrialn[_num_sp][spt] = rhs_cto[ind];
886  ind++;
887  }
888 
889  const std::vector<std::vector<Real>> dvar_dtrial_old = dvar_dtrial;
890 
891  for (unsigned v = 0; v < _num_sp; ++v) // loop over variables in NR procedure
892  {
893  for (unsigned spt = 0; spt < _num_sp; ++spt) // loop over trial stress-param variables
894  {
895  dvar_dtrial[v][spt] = step_size * dvarn_dtrialn[v][spt];
896  for (unsigned a = 0; a < _num_sp; ++a)
897  dvar_dtrial[v][spt] += dvarn_dtrialn[v][a] * dvar_dtrial_old[a][spt];
898  }
899  }
900  // for gaE the formulae are a little different
901  const unsigned v = _num_sp;
902  for (unsigned spt = 0; spt < _num_sp; ++spt)
903  {
904  dvar_dtrial[v][spt] += step_size * dvarn_dtrialn[v][spt]; // note +=
905  for (unsigned a = 0; a < _num_sp; ++a)
906  dvar_dtrial[v][spt] += dvarn_dtrialn[v][a] * dvar_dtrial_old[a][spt];
907  }
908 }
909 
910 bool
911 MultiParameterPlasticityStressUpdate::precisionLoss(const std::vector<Real> & solution,
912  const std::vector<Real> & stress_params,
913  Real gaE) const
914 {
915  if (std::abs(solution[_num_sp]) > 1E-13 * std::abs(gaE))
916  return false;
917  for (unsigned i = 0; i < _num_sp; ++i)
918  if (std::abs(solution[i]) > 1E-13 * std::abs(stress_params[i]))
919  return false;
920  return true;
921 }
int nrStep(const yieldAndFlow &smoothed_q, const std::vector< Real > &trial_stress_params, const std::vector< Real > &stress_params, const std::vector< Real > &intnl, Real gaE, std::vector< Real > &rhs) const
Performs a Newton-Raphson step to attempt to zero rhs Upon return, rhs will contain the solution...
virtual void initializeVarsV(const std::vector< Real > &trial_stress_params, const std::vector< Real > &intnl_old, std::vector< Real > &stress_params, Real &gaE, std::vector< Real > &intnl) const
Sets (stress_params, intnl) at "good guesses" of the solution to the Return-Map algorithm.
MaterialProperty< std::vector< Real > > & _yf
yield functions
std::vector< Real > _trial_sp
"Trial" value of stress_params that initializes the return-map process This is derived from stress = ...
virtual void computeAllQV(const std::vector< Real > &stress_params, const std::vector< Real > &intnl, std::vector< yieldAndFlow > &all_q) const =0
Completely fills all_q with correct values.
MultiParameterPlasticityStressUpdate(const InputParameters &parameters, unsigned num_sp, unsigned num_yf, unsigned num_intnl)
virtual void setEffectiveElasticity(const RankFourTensor &Eijkl)=0
Sets _Eij and _En and _Cij.
const bool _warn_about_precision_loss
Output a warning message if precision loss is encountered during the return-map process.
Real smoother(Real f_diff) const
Derivative of ismoother.
StressUpdateBase is a material that is not called by MOOSE because of the compute=false flag set in t...
const unsigned _num_intnl
Number of internal parameters.
Struct designed to hold info about a single yield function and its derivatives, as well as the flow d...
void calculateRHS(const std::vector< Real > &trial_stress_params, const std::vector< Real > &stress_params, Real gaE, const yieldAndFlow &smoothed_q, std::vector< Real > &rhs) const
Calculates the RHS in the following 0 = rhs[0] = S[0] - S[0]^trial + ga * E[0, j] * dg/dS[j] 0 = rhs[...
std::vector< std::vector< Real > > _Cij
_Cij[i, j] * _Eij[j, k] = 1 iff j == k
virtual void consistentTangentOperatorV(const RankTwoTensor &stress_trial, const std::vector< Real > &trial_stress_params, const RankTwoTensor &stress, const std::vector< Real > &stress_params, Real gaE, const yieldAndFlow &smoothed_q, const RankFourTensor &Eijkl, bool compute_full_tangent_operator, const std::vector< std::vector< Real >> &dvar_dtrial, RankFourTensor &cto)
Calculates the consistent tangent operator.
virtual void setStressAfterReturnV(const RankTwoTensor &stress_trial, const std::vector< Real > &stress_params, Real gaE, const std::vector< Real > &intnl, const yieldAndFlow &smoothed_q, const RankFourTensor &Eijkl, RankTwoTensor &stress) const =0
Sets stress from the admissible parameters.
virtual void yieldFunctionValuesV(const std::vector< Real > &stress_params, const std::vector< Real > &intnl, std::vector< Real > &yf) const =0
Computes the values of the yield functions, given stress_params and intnl parameters.
std::vector< Real > _ok_intnl
The state (ok_sp, ok_intnl) is known to be admissible.
const Real _f_tol2
Square of the yield-function tolerance.
const Real _min_step_size
In order to help the Newton-Raphson procedure, the applied strain increment may be applied in sub-inc...
virtual void updateState(RankTwoTensor &strain_increment, RankTwoTensor &inelastic_strain_increment, const RankTwoTensor &rotation_increment, RankTwoTensor &stress_new, const RankTwoTensor &stress_old, const RankFourTensor &elasticity_tensor, const RankTwoTensor &elastic_strain_old, bool compute_full_tangent_operator, RankFourTensor &tangent_operator) override
Given a strain increment that results in a trial stress, perform some procedure (such as an iterative...
MaterialProperty< std::vector< Real > > & _intnl
internal parameters
std::vector< std::vector< Real > > _Eij
E[i, j] in the system of equations to be solved.
const Real _f_tol
The yield-function tolerance.
const unsigned _num_yf
Number of yield functions.
int lineSearch(Real &res2, std::vector< Real > &stress_params, Real &gaE, const std::vector< Real > &trial_stress_params, yieldAndFlow &smoothed_q, const std::vector< Real > &intnl_ok, std::vector< Real > &intnl, std::vector< Real > &rhs, Real &linesearch_needed) const
Performs a line-search to find stress_params and gaE Upon entry:
Real ismoother(Real f_diff) const
Smooths yield functions.
const std::vector< Real > _definitely_ok_sp
An admissible value of stress_params at the initial time.
const MaterialProperty< std::vector< Real > > & _intnl_old
old values of internal parameters
RankTwoTensor _stress_trial
"Trial" value of stress that is set at the beginning of the return-map process.
virtual void propagateQpStatefulProperties() override
If updateState is not called during a timestep, this will be.
virtual void setInelasticStrainIncrementAfterReturn(const RankTwoTensor &stress_trial, Real gaE, const yieldAndFlow &smoothed_q, const RankFourTensor &elasticity_tensor, const RankTwoTensor &returned_stress, RankTwoTensor &inelastic_strain_increment) const
Sets inelastic strain increment from the returned configuration This is called after the return-map p...
yieldAndFlow smoothAllQuantities(const std::vector< Real > &stress_params, const std::vector< Real > &intnl) const
Calculates all yield functions and derivatives, and then performs the smoothing scheme.
virtual void computeStressParams(const RankTwoTensor &stress, std::vector< Real > &stress_params) const =0
Computes stress_params, given stress.
InputParameters validParams< StressUpdateBase >()
MaterialProperty< Real > & _iter
Number of Newton-Raphson iterations used in the return-map.
MaterialProperty< RankTwoTensor > & _plastic_strain
plastic strain
const Real _smoothing_tol
Smoothing tolerance: edges of the yield surface get smoothed by this amount.
Real dsmoother(Real f_diff) const
Derivative of smoother.
virtual std::vector< RankTwoTensor > dstress_param_dstress(const RankTwoTensor &stress) const =0
d(stress_param[i])/d(stress) at given stress
std::vector< Real > _current_intnl
The current values of the internal params during the Newton-Raphson.
virtual void finalizeReturnProcess(const RankTwoTensor &rotation_increment)
Derived classes may use this to perform calculations after the return-map process has completed succe...
InputParameters validParams< MultiParameterPlasticityStressUpdate >()
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)
Real yieldF(const std::vector< Real > &stress_params, const std::vector< Real > &intnl) const
Computes the smoothed yield function.
virtual std::vector< RankFourTensor > d2stress_param_dstress(const RankTwoTensor &stress) const =0
d2(stress_param[i])/d(stress)/d(stress) at given stress
virtual void setIntnlDerivativesV(const std::vector< Real > &trial_stress_params, const std::vector< Real > &current_stress_params, const std::vector< Real > &intnl, std::vector< std::vector< Real >> &dintnl) const =0
Sets the derivatives of internal parameters, based on the trial values of stress_params, their current values, and the current values of the internal parameters.
bool precisionLoss(const std::vector< Real > &solution, const std::vector< Real > &stress_params, Real gaE) const
Check whether precision loss has occurred.
MaterialProperty< Real > & _linesearch_needed
Whether a line-search was needed in the latest Newton-Raphson process (1 if true, 0 otherwise) ...
std::vector< std::vector< Real > > _dvar_dtrial
d({stress_param[i], gaE})/d(trial_stress_param[j])
Real calculateRes2(const std::vector< Real > &rhs) const
Calculates the residual-squared for the Newton-Raphson + line-search.
virtual void setIntnlValuesV(const std::vector< Real > &trial_stress_params, const std::vector< Real > &current_stress_params, const std::vector< Real > &intnl_old, std::vector< Real > &intnl) const =0
Sets the internal parameters based on the trial values of stress_params, their current values...
std::vector< Real > _current_sp
The current values of the stress params during the Newton-Raphson.
bool _step_one
handles case of initial_stress that is inadmissible being supplied
virtual void errorHandler(const std::string &message) const
Performs any necessary cleaning-up, then throw MooseException(message)
void dnRHSdVar(const yieldAndFlow &smoothed_q, const std::vector< std::vector< Real >> &dintnl, const std::vector< Real > &stress_params, Real gaE, std::vector< double > &jac) const
Derivative of -RHS with respect to the stress_params and gaE, placed into an array ready for solving ...
virtual void preReturnMapV(const std::vector< Real > &trial_stress_params, const RankTwoTensor &stress_trial, const std::vector< Real > &intnl_old, const std::vector< Real > &yf, const RankFourTensor &Eijkl)
Derived classes may employ this function to record stuff or do other computations prior to the return...
std::vector< Real > _del_stress_params
_del_stress_params = trial_stress_params - ok_sp This is fixed at the beginning of the return-map pro...
const unsigned _max_nr_its
Maximum number of Newton-Raphson iterations allowed in the return-map process.
virtual void initializeReturnProcess()
Derived classes may use this to perform calculations before any return-map process is performed...
const MaterialProperty< RankTwoTensor > & _plastic_strain_old
Old value of plastic strain.
std::vector< Real > _rhs
0 = rhs[0] = S[0] - S[0]^trial + ga * E[0, i] * dg/dS[i] 0 = rhs[1] = S[1] - S[1]^trial + ga * E[1...
std::vector< Real > _ok_sp
The state (ok_sp, ok_intnl) is known to be admissible, so ok_sp are stress_params that are "OK"...
void dVardTrial(bool elastic_only, const std::vector< Real > &trial_stress_params, const std::vector< Real > &stress_params, Real gaE, const std::vector< Real > &intnl, const yieldAndFlow &smoothed_q, Real step_size, bool compute_full_tangent_operator, std::vector< std::vector< Real >> &dvar_dtrial) const
Calculates derivatives of the stress_params and gaE with repect to the trial values of the stress_par...
const unsigned _num_sp
Number of stress parameters.