www.mooseframework.org
MultiPlasticityDebugger.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 "libmesh/utility.h"
9 
10 template <>
11 InputParameters
13 {
14  InputParameters params = validParams<MultiPlasticityLinearSystem>();
15  MooseEnum debug_fspb_type("none crash jacobian jacobian_and_linear_system", "none");
16  params.addParam<MooseEnum>("debug_fspb",
17  debug_fspb_type,
18  "Debug types for use by developers when creating new "
19  "plasticity models, not for general use. 2 = debug Jacobian "
20  "entries, 3 = check the entire Jacobian, and check Ax=b");
21  params.addParam<RealTensorValue>("debug_jac_at_stress",
22  RealTensorValue(),
23  "Debug Jacobian entries at this stress. For use by developers");
24  params.addParam<std::vector<Real>>("debug_jac_at_pm",
25  "Debug Jacobian entries at these plastic multipliers");
26  params.addParam<std::vector<Real>>("debug_jac_at_intnl",
27  "Debug Jacobian entries at these internal parameters");
28  params.addParam<Real>(
29  "debug_stress_change", 1.0, "Debug finite differencing parameter for the stress");
30  params.addParam<std::vector<Real>>(
31  "debug_pm_change", "Debug finite differencing parameters for the plastic multipliers");
32  params.addParam<std::vector<Real>>(
33  "debug_intnl_change", "Debug finite differencing parameters for the internal parameters");
34  return params;
35 }
36 
37 MultiPlasticityDebugger::MultiPlasticityDebugger(const MooseObject * moose_object)
38  : MultiPlasticityLinearSystem(moose_object),
39  _fspb_debug(_params.get<MooseEnum>("debug_fspb")),
40  _fspb_debug_stress(_params.get<RealTensorValue>("debug_jac_at_stress")),
41  _fspb_debug_pm(_params.get<std::vector<Real>>("debug_jac_at_pm")),
42  _fspb_debug_intnl(_params.get<std::vector<Real>>("debug_jac_at_intnl")),
43  _fspb_debug_stress_change(_params.get<Real>("debug_stress_change")),
44  _fspb_debug_pm_change(_params.get<std::vector<Real>>("debug_pm_change")),
45  _fspb_debug_intnl_change(_params.get<std::vector<Real>>("debug_intnl_change"))
46 {
47 }
48 
49 void
51 {
52  Moose::err << "Debug Parameters are as follows\n";
53  Moose::err << "stress = \n";
54  _fspb_debug_stress.print();
55 
56  if (_fspb_debug_pm.size() != _num_surfaces || _fspb_debug_intnl.size() != _num_models ||
59  mooseError("The debug parameters have the wrong size\n");
60 
61  Moose::err << "plastic multipliers =\n";
62  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
63  Moose::err << _fspb_debug_pm[surface] << "\n";
64 
65  Moose::err << "internal parameters =\n";
66  for (unsigned model = 0; model < _num_models; ++model)
67  Moose::err << _fspb_debug_intnl[model] << "\n";
68 
69  Moose::err << "finite-differencing parameter for stress-changes:\n"
70  << _fspb_debug_stress_change << "\n";
71  Moose::err << "finite-differencing parameter(s) for plastic-multiplier(s):\n";
72  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
73  Moose::err << _fspb_debug_pm_change[surface] << "\n";
74  Moose::err << "finite-differencing parameter(s) for internal-parameter(s):\n";
75  for (unsigned model = 0; model < _num_models; ++model)
76  Moose::err << _fspb_debug_intnl_change[model] << "\n";
77 }
78 
79 void
81 {
82  Moose::err
83  << "\n\n++++++++++++++++++++++++\nChecking the derivatives\n++++++++++++++++++++++++\n";
85 
86  std::vector<bool> act;
87  act.assign(_num_surfaces, true);
88 
89  Moose::err << "\ndyieldFunction_dstress. Relative L2 norms.\n";
90  std::vector<RankTwoTensor> df_dstress;
91  std::vector<RankTwoTensor> fddf_dstress;
94  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
95  {
96  Moose::err << "surface = " << surface << " Relative L2norm = "
97  << 2 * (df_dstress[surface] - fddf_dstress[surface]).L2norm() /
98  (df_dstress[surface] + fddf_dstress[surface]).L2norm()
99  << "\n";
100  Moose::err << "Coded:\n";
101  df_dstress[surface].print();
102  Moose::err << "Finite difference:\n";
103  fddf_dstress[surface].print();
104  }
105 
106  Moose::err << "\ndyieldFunction_dintnl.\n";
107  std::vector<Real> df_dintnl;
109  Moose::err << "Coded:\n";
110  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
111  Moose::err << df_dintnl[surface] << " ";
112  Moose::err << "\n";
113  std::vector<Real> fddf_dintnl;
115  Moose::err << "Finite difference:\n";
116  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
117  Moose::err << fddf_dintnl[surface] << " ";
118  Moose::err << "\n";
119 
120  Moose::err << "\ndflowPotential_dstress. Relative L2 norms.\n";
121  std::vector<RankFourTensor> dr_dstress;
122  std::vector<RankFourTensor> fddr_dstress;
125  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
126  {
127  Moose::err << "surface = " << surface << " Relative L2norm = "
128  << 2 * (dr_dstress[surface] - fddr_dstress[surface]).L2norm() /
129  (dr_dstress[surface] + fddr_dstress[surface]).L2norm()
130  << "\n";
131  Moose::err << "Coded:\n";
132  dr_dstress[surface].print();
133  Moose::err << "Finite difference:\n";
134  fddr_dstress[surface].print();
135  }
136 
137  Moose::err << "\ndflowPotential_dintnl. Relative L2 norms.\n";
138  std::vector<RankTwoTensor> dr_dintnl;
139  std::vector<RankTwoTensor> fddr_dintnl;
142  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
143  {
144  Moose::err << "surface = " << surface << " Relative L2norm = "
145  << 2 * (dr_dintnl[surface] - fddr_dintnl[surface]).L2norm() /
146  (dr_dintnl[surface] + fddr_dintnl[surface]).L2norm()
147  << "\n";
148  Moose::err << "Coded:\n";
149  dr_dintnl[surface].print();
150  Moose::err << "Finite difference:\n";
151  fddr_dintnl[surface].print();
152  }
153 }
154 
155 void
156 MultiPlasticityDebugger::checkJacobian(const RankFourTensor & E_inv,
157  const std::vector<Real> & intnl_old)
158 {
159  Moose::err << "\n\n+++++++++++++++++++++\nChecking the Jacobian\n+++++++++++++++++++++\n";
161 
162  std::vector<bool> act;
163  act.assign(_num_surfaces, true);
164  std::vector<bool> deactivated_due_to_ld;
165  deactivated_due_to_ld.assign(_num_surfaces, false);
166 
167  RankTwoTensor delta_dp = -E_inv * _fspb_debug_stress;
168 
169  std::vector<std::vector<Real>> jac;
170  calculateJacobian(_fspb_debug_stress,
173  E_inv,
174  act,
175  deactivated_due_to_ld,
176  jac);
177 
178  std::vector<std::vector<Real>> fdjac;
179  fdJacobian(_fspb_debug_stress,
180  intnl_old,
183  delta_dp,
184  E_inv,
185  false,
186  fdjac);
187 
188  Real L2_numer = 0;
189  Real L2_denom = 0;
190  for (unsigned row = 0; row < jac.size(); ++row)
191  for (unsigned col = 0; col < jac.size(); ++col)
192  {
193  L2_numer += Utility::pow<2>(jac[row][col] - fdjac[row][col]);
194  L2_denom += Utility::pow<2>(jac[row][col] + fdjac[row][col]);
195  }
196  Moose::err << "\nRelative L2norm = " << std::sqrt(L2_numer / L2_denom) / 0.5 << "\n";
197 
198  Moose::err << "\nHand-coded Jacobian:\n";
199  for (unsigned row = 0; row < jac.size(); ++row)
200  {
201  for (unsigned col = 0; col < jac.size(); ++col)
202  Moose::err << jac[row][col] << " ";
203  Moose::err << "\n";
204  }
205 
206  Moose::err << "Finite difference Jacobian:\n";
207  for (unsigned row = 0; row < fdjac.size(); ++row)
208  {
209  for (unsigned col = 0; col < fdjac.size(); ++col)
210  Moose::err << fdjac[row][col] << " ";
211  Moose::err << "\n";
212  }
213 }
214 
215 void
216 MultiPlasticityDebugger::fdJacobian(const RankTwoTensor & stress,
217  const std::vector<Real> & intnl_old,
218  const std::vector<Real> & intnl,
219  const std::vector<Real> & pm,
220  const RankTwoTensor & delta_dp,
221  const RankFourTensor & E_inv,
222  bool eliminate_ld,
223  std::vector<std::vector<Real>> & jac)
224 {
225  std::vector<bool> active;
226  active.assign(_num_surfaces, true);
227 
228  std::vector<bool> deactivated_due_to_ld;
229  std::vector<bool> deactivated_due_to_ld_ep;
230 
231  std::vector<Real> orig_rhs;
232  calculateRHS(stress,
233  intnl_old,
234  intnl,
235  pm,
236  delta_dp,
237  orig_rhs,
238  active,
239  eliminate_ld,
240  deactivated_due_to_ld); // this calculates RHS, and also set deactivated_due_to_ld.
241  // The latter stays fixed for the rest of this routine
242 
243  unsigned int whole_system_size = 6 + _num_surfaces + _num_models;
244  unsigned int system_size =
245  orig_rhs.size(); // will be = whole_system_size if eliminate_ld = false, since all active=true
246  jac.resize(system_size);
247  for (unsigned row = 0; row < system_size; ++row)
248  jac[row].assign(system_size, 0);
249 
250  std::vector<Real> rhs_ep;
251  unsigned col = 0;
252 
253  RankTwoTensor stressep;
254  RankTwoTensor delta_dpep;
255  Real ep = _fspb_debug_stress_change;
256  for (unsigned i = 0; i < 3; ++i)
257  for (unsigned j = 0; j <= i; ++j)
258  {
259  stressep = stress;
260  stressep(i, j) += ep;
261  if (i != j)
262  stressep(j, i) += ep;
263  delta_dpep = delta_dp;
264  for (unsigned k = 0; k < 3; ++k)
265  for (unsigned l = 0; l < 3; ++l)
266  {
267  delta_dpep(k, l) -= E_inv(k, l, i, j) * ep;
268  if (i != j)
269  delta_dpep(k, l) -= E_inv(k, l, j, i) * ep;
270  }
271  active.assign(_num_surfaces, true);
272  calculateRHS(stressep,
273  intnl_old,
274  intnl,
275  pm,
276  delta_dpep,
277  rhs_ep,
278  active,
279  false,
280  deactivated_due_to_ld_ep);
281  unsigned row = 0;
282  for (unsigned dof = 0; dof < whole_system_size; ++dof)
283  if (dof_included(dof, deactivated_due_to_ld))
284  {
285  jac[row][col] =
286  -(rhs_ep[dof] - orig_rhs[row]) / ep; // remember jacobian = -d(rhs)/d(something)
287  row++;
288  }
289  col++; // all of the first 6 columns are dof_included since they're stresses
290  }
291 
292  std::vector<Real> pmep;
293  pmep.resize(_num_surfaces);
294  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
295  pmep[surface] = pm[surface];
296  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
297  {
298  if (!dof_included(6 + surface, deactivated_due_to_ld))
299  continue;
300  ep = _fspb_debug_pm_change[surface];
301  pmep[surface] += ep;
302  active.assign(_num_surfaces, true);
303  calculateRHS(
304  stress, intnl_old, intnl, pmep, delta_dp, rhs_ep, active, false, deactivated_due_to_ld_ep);
305  unsigned row = 0;
306  for (unsigned dof = 0; dof < whole_system_size; ++dof)
307  if (dof_included(dof, deactivated_due_to_ld))
308  {
309  jac[row][col] =
310  -(rhs_ep[dof] - orig_rhs[row]) / ep; // remember jacobian = -d(rhs)/d(something)
311  row++;
312  }
313  pmep[surface] -= ep;
314  col++;
315  }
316 
317  std::vector<Real> intnlep;
318  intnlep.resize(_num_models);
319  for (unsigned model = 0; model < _num_models; ++model)
320  intnlep[model] = intnl[model];
321  for (unsigned model = 0; model < _num_models; ++model)
322  {
323  if (!dof_included(6 + _num_surfaces + model, deactivated_due_to_ld))
324  continue;
325  ep = _fspb_debug_intnl_change[model];
326  intnlep[model] += ep;
327  active.assign(_num_surfaces, true);
328  calculateRHS(
329  stress, intnl_old, intnlep, pm, delta_dp, rhs_ep, active, false, deactivated_due_to_ld_ep);
330  unsigned row = 0;
331  for (unsigned dof = 0; dof < whole_system_size; ++dof)
332  if (dof_included(dof, deactivated_due_to_ld))
333  {
334  jac[row][col] =
335  -(rhs_ep[dof] - orig_rhs[row]) / ep; // remember jacobian = -d(rhs)/d(something)
336  row++;
337  }
338  intnlep[model] -= ep;
339  col++;
340  }
341 }
342 
343 bool
345  const std::vector<bool> & deactivated_due_to_ld)
346 {
347  if (dof < unsigned(6))
348  // these are the stress components
349  return true;
350  unsigned eff_dof = dof - 6;
351  if (eff_dof < _num_surfaces)
352  // these are the plastic multipliers, pm
353  return !deactivated_due_to_ld[eff_dof];
354  eff_dof -= _num_surfaces; // now we know the dof is an intnl parameter
355  std::vector<bool> active_surface(_num_surfaces);
356  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
357  active_surface[surface] = !deactivated_due_to_ld[surface];
358  return anyActiveSurfaces(eff_dof, active_surface);
359 }
360 
361 void
362 MultiPlasticityDebugger::checkSolution(const RankFourTensor & E_inv)
363 {
364  Moose::err << "\n\n+++++++++++++++++++++\nChecking the Solution\n";
365  Moose::err << "(Ie, checking Ax = b)\n+++++++++++++++++++++\n";
367 
368  std::vector<bool> act;
369  act.assign(_num_surfaces, true);
370  std::vector<bool> deactivated_due_to_ld;
371  deactivated_due_to_ld.assign(_num_surfaces, false);
372 
373  RankTwoTensor delta_dp = -E_inv * _fspb_debug_stress;
374 
375  std::vector<Real> orig_rhs;
376  calculateRHS(_fspb_debug_stress,
380  delta_dp,
381  orig_rhs,
382  act,
383  true,
384  deactivated_due_to_ld);
385 
386  Moose::err << "\nb = ";
387  for (unsigned i = 0; i < orig_rhs.size(); ++i)
388  Moose::err << orig_rhs[i] << " ";
389  Moose::err << "\n\n";
390 
391  std::vector<std::vector<Real>> jac_coded;
392  calculateJacobian(_fspb_debug_stress,
395  E_inv,
396  act,
397  deactivated_due_to_ld,
398  jac_coded);
399 
400  Moose::err
401  << "Before checking Ax=b is correct, check that the Jacobians given below are equal.\n";
402  Moose::err
403  << "The hand-coded Jacobian is used in calculating the solution 'x', given 'b' above.\n";
404  Moose::err << "Note that this only includes degrees of freedom that aren't deactivated due to "
405  "linear dependence.\n";
406  Moose::err << "Hand-coded Jacobian:\n";
407  for (unsigned row = 0; row < jac_coded.size(); ++row)
408  {
409  for (unsigned col = 0; col < jac_coded.size(); ++col)
410  Moose::err << jac_coded[row][col] << " ";
411  Moose::err << "\n";
412  }
413 
414  deactivated_due_to_ld.assign(_num_surfaces,
415  false); // this potentially gets changed by nrStep, below
416  RankTwoTensor dstress;
417  std::vector<Real> dpm;
418  std::vector<Real> dintnl;
419  nrStep(_fspb_debug_stress,
423  E_inv,
424  delta_dp,
425  dstress,
426  dpm,
427  dintnl,
428  act,
429  deactivated_due_to_ld);
430 
431  std::vector<bool> active_not_deact(_num_surfaces);
432  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
433  active_not_deact[surface] = !deactivated_due_to_ld[surface];
434 
435  std::vector<Real> x;
436  x.assign(orig_rhs.size(), 0);
437  unsigned ind = 0;
438  for (unsigned i = 0; i < 3; ++i)
439  for (unsigned j = 0; j <= i; ++j)
440  x[ind++] = dstress(i, j);
441  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
442  if (active_not_deact[surface])
443  x[ind++] = dpm[surface];
444  for (unsigned model = 0; model < _num_models; ++model)
445  if (anyActiveSurfaces(model, active_not_deact))
446  x[ind++] = dintnl[model];
447 
448  mooseAssert(ind == orig_rhs.size(),
449  "Incorrect extracting of changes from NR solution in the "
450  "finite-difference checking of nrStep");
451 
452  Moose::err << "\nThis yields x =";
453  for (unsigned i = 0; i < orig_rhs.size(); ++i)
454  Moose::err << x[i] << " ";
455  Moose::err << "\n";
456 
457  std::vector<std::vector<Real>> jac_fd;
458  fdJacobian(_fspb_debug_stress,
462  delta_dp,
463  E_inv,
464  true,
465  jac_fd);
466 
467  Moose::err << "\nThe finite-difference Jacobian is used to multiply by this 'x',\n";
468  Moose::err << "in order to check that the solution is correct\n";
469  Moose::err << "Finite-difference Jacobian:\n";
470  for (unsigned row = 0; row < jac_fd.size(); ++row)
471  {
472  for (unsigned col = 0; col < jac_fd.size(); ++col)
473  Moose::err << jac_fd[row][col] << " ";
474  Moose::err << "\n";
475  }
476 
477  Real L2_numer = 0;
478  Real L2_denom = 0;
479  for (unsigned row = 0; row < jac_coded.size(); ++row)
480  for (unsigned col = 0; col < jac_coded.size(); ++col)
481  {
482  L2_numer += Utility::pow<2>(jac_coded[row][col] - jac_fd[row][col]);
483  L2_denom += Utility::pow<2>(jac_coded[row][col] + jac_fd[row][col]);
484  }
485  Moose::err << "Relative L2norm of the hand-coded and finite-difference Jacobian is "
486  << std::sqrt(L2_numer / L2_denom) / 0.5 << "\n";
487 
488  std::vector<Real> fd_times_x;
489  fd_times_x.assign(orig_rhs.size(), 0);
490  for (unsigned row = 0; row < orig_rhs.size(); ++row)
491  for (unsigned col = 0; col < orig_rhs.size(); ++col)
492  fd_times_x[row] += jac_fd[row][col] * x[col];
493 
494  Moose::err << "\n(Finite-difference Jacobian)*x =\n";
495  for (unsigned i = 0; i < orig_rhs.size(); ++i)
496  Moose::err << fd_times_x[i] << " ";
497  Moose::err << "\n";
498  Moose::err << "Recall that b = \n";
499  for (unsigned i = 0; i < orig_rhs.size(); ++i)
500  Moose::err << orig_rhs[i] << " ";
501  Moose::err << "\n";
502 
503  L2_numer = 0;
504  L2_denom = 0;
505  for (unsigned i = 0; i < orig_rhs.size(); ++i)
506  {
507  L2_numer += Utility::pow<2>(orig_rhs[i] - fd_times_x[i]);
508  L2_denom += Utility::pow<2>(orig_rhs[i] + fd_times_x[i]);
509  }
510  Moose::err << "\nRelative L2norm of these is " << std::sqrt(L2_numer / L2_denom) / 0.5 << "\n";
511 }
512 
513 void
515  const std::vector<Real> & intnl,
516  std::vector<RankTwoTensor> & df_dstress)
517 {
518  df_dstress.assign(_num_surfaces, RankTwoTensor());
519 
520  std::vector<bool> act;
521  act.assign(_num_surfaces, true);
522 
523  Real ep = _fspb_debug_stress_change;
524  RankTwoTensor stressep;
525  std::vector<Real> fep, fep_minus;
526  for (unsigned i = 0; i < 3; ++i)
527  for (unsigned j = 0; j < 3; ++j)
528  {
529  stressep = stress;
530  // do a central difference to attempt to capture discontinuities
531  // such as those encountered in tensile and Mohr-Coulomb
532  stressep(i, j) += ep / 2.0;
533  yieldFunction(stressep, intnl, act, fep);
534  stressep(i, j) -= ep;
535  yieldFunction(stressep, intnl, act, fep_minus);
536  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
537  df_dstress[surface](i, j) = (fep[surface] - fep_minus[surface]) / ep;
538  }
539 }
540 
541 void
543  const std::vector<Real> & intnl,
544  std::vector<Real> & df_dintnl)
545 {
546  df_dintnl.resize(_num_surfaces);
547 
548  std::vector<bool> act;
549  act.assign(_num_surfaces, true);
550 
551  std::vector<Real> origf;
552  yieldFunction(stress, intnl, act, origf);
553 
554  std::vector<Real> intnlep;
555  intnlep.resize(_num_models);
556  for (unsigned model = 0; model < _num_models; ++model)
557  intnlep[model] = intnl[model];
558  Real ep;
559  std::vector<Real> fep;
560  unsigned int model;
561  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
562  {
563  model = modelNumber(surface);
564  ep = _fspb_debug_intnl_change[model];
565  intnlep[model] += ep;
566  yieldFunction(stress, intnlep, act, fep);
567  df_dintnl[surface] = (fep[surface] - origf[surface]) / ep;
568  intnlep[model] -= ep;
569  }
570 }
571 
572 void
574  const std::vector<Real> & intnl,
575  std::vector<RankFourTensor> & dr_dstress)
576 {
577  dr_dstress.assign(_num_surfaces, RankFourTensor());
578 
579  std::vector<bool> act;
580  act.assign(_num_surfaces, true);
581 
582  Real ep = _fspb_debug_stress_change;
583  RankTwoTensor stressep;
584  std::vector<RankTwoTensor> rep, rep_minus;
585  for (unsigned i = 0; i < 3; ++i)
586  for (unsigned j = 0; j < 3; ++j)
587  {
588  stressep = stress;
589  // do a central difference
590  stressep(i, j) += ep / 2.0;
591  flowPotential(stressep, intnl, act, rep);
592  stressep(i, j) -= ep;
593  flowPotential(stressep, intnl, act, rep_minus);
594  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
595  for (unsigned k = 0; k < 3; ++k)
596  for (unsigned l = 0; l < 3; ++l)
597  dr_dstress[surface](k, l, i, j) = (rep[surface](k, l) - rep_minus[surface](k, l)) / ep;
598  }
599 }
600 
601 void
603  const std::vector<Real> & intnl,
604  std::vector<RankTwoTensor> & dr_dintnl)
605 {
606  dr_dintnl.resize(_num_surfaces);
607 
608  std::vector<bool> act;
609  act.assign(_num_surfaces, true);
610 
611  std::vector<RankTwoTensor> origr;
612  flowPotential(stress, intnl, act, origr);
613 
614  std::vector<Real> intnlep;
615  intnlep.resize(_num_models);
616  for (unsigned model = 0; model < _num_models; ++model)
617  intnlep[model] = intnl[model];
618  Real ep;
619  std::vector<RankTwoTensor> rep;
620  unsigned int model;
621  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
622  {
623  model = modelNumber(surface);
624  ep = _fspb_debug_intnl_change[model];
625  intnlep[model] += ep;
626  flowPotential(stress, intnlep, act, rep);
627  dr_dintnl[surface] = (rep[surface] - origr[surface]) / ep;
628  intnlep[model] -= ep;
629  }
630 }
void fddyieldFunction_dintnl(const RankTwoTensor &stress, const std::vector< Real > &intnl, std::vector< Real > &df_dintnl)
The finite-difference derivative of yield function(s) with respect to internal parameter(s) ...
virtual void fdJacobian(const RankTwoTensor &stress, const std::vector< Real > &intnl_old, const std::vector< Real > &intnl, const std::vector< Real > &pm, const RankTwoTensor &delta_dp, const RankFourTensor &E_inv, bool eliminate_ld, std::vector< std::vector< Real >> &jac)
The Jacobian calculated using finite differences.
void fddyieldFunction_dstress(const RankTwoTensor &stress, const std::vector< Real > &intnl, std::vector< RankTwoTensor > &df_dstress)
The finite-difference derivative of yield function(s) with respect to stress.
virtual void calculateRHS(const RankTwoTensor &stress, const std::vector< Real > &intnl_old, const std::vector< Real > &intnl, const std::vector< Real > &pm, const RankTwoTensor &delta_dp, std::vector< Real > &rhs, const std::vector< bool > &active, bool eliminate_ld, std::vector< bool > &deactivated_due_to_ld)
Calculate the RHS which is rhs = -(epp(0,0), epp(1,0), epp(1,1), epp(2,0), epp(2,1), epp(2,2), f[0], f[1], ..., f[num_f], ic[0], ic[1], ..., ic[num_ic])
virtual void dyieldFunction_dstress(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< RankTwoTensor > &df_dstress)
The derivative of the active yield function(s) with respect to stress.
unsigned int modelNumber(unsigned int surface)
returns the model number, given the surface number
virtual void calculateJacobian(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< Real > &pm, const RankFourTensor &E_inv, const std::vector< bool > &active, const std::vector< bool > &deactivated_due_to_ld, std::vector< std::vector< Real >> &jac)
d(rhs)/d(dof)
Real _fspb_debug_stress_change
Debug finite-differencing parameter for the stress.
virtual void fddflowPotential_dintnl(const RankTwoTensor &stress, const std::vector< Real > &intnl, std::vector< RankTwoTensor > &dr_dintnl)
The finite-difference derivative of the flow potentials with respect to internal parameters.
InputParameters validParams< MultiPlasticityDebugger >()
std::vector< Real > _fspb_debug_intnl_change
Debug finite-differencing parameters for the internal parameters.
virtual void yieldFunction(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< Real > &f)
The active yield function(s)
virtual void fddflowPotential_dstress(const RankTwoTensor &stress, const std::vector< Real > &intnl, std::vector< RankFourTensor > &dr_dstress)
The finite-difference derivative of the flow potential(s) with respect to stress. ...
std::vector< Real > _fspb_debug_intnl
Debug the Jacobian entires at these internal parameters.
MultiPlasticityLinearSystem computes the linear system and handles linear-dependence removal for use ...
unsigned int _num_surfaces
Number of surfaces within the plastic models.
void outputAndCheckDebugParameters()
Outputs the debug parameters: _fspb_debug_stress, _fspd_debug_pm, etc and checks that they are sized ...
virtual void dflowPotential_dintnl(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< RankTwoTensor > &dr_dintnl)
The derivative of the active flow potentials with respect to the active internal parameters The UserO...
void checkSolution(const RankFourTensor &E_inv)
Checks that Ax does equal b in the NR procedure.
void checkDerivatives()
Checks the derivatives, eg dyieldFunction_dstress by using finite difference approximations.
RankTwoTensor _fspb_debug_stress
Debug the Jacobian entries at this stress.
bool anyActiveSurfaces(int model, const std::vector< bool > &active)
returns true if any internal surfaces of the given model are active according to &#39;active&#39; ...
virtual void flowPotential(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< RankTwoTensor > &r)
The active flow potential(s) - one for each yield function.
bool dof_included(unsigned int dof, const std::vector< bool > &deactivated_due_to_ld)
std::vector< Real > _fspb_debug_pm_change
Debug finite-differencing parameters for the plastic multipliers.
unsigned int _num_models
Number of plastic models for this material.
virtual void dflowPotential_dstress(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< RankFourTensor > &dr_dstress)
The derivative of the active flow potential(s) with respect to stress.
std::vector< Real > _fspb_debug_pm
Debug the Jacobian entires at these plastic multipliers.
Real L2norm(const RankTwoTensor &r2tensor)
InputParameters validParams< MultiPlasticityLinearSystem >()
virtual void dyieldFunction_dintnl(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< Real > &df_dintnl)
The derivative of active yield function(s) with respect to their internal parameters (the user object...
MultiPlasticityDebugger(const MooseObject *moose_object)
virtual void nrStep(const RankTwoTensor &stress, const std::vector< Real > &intnl_old, const std::vector< Real > &intnl, const std::vector< Real > &pm, const RankFourTensor &E_inv, const RankTwoTensor &delta_dp, RankTwoTensor &dstress, std::vector< Real > &dpm, std::vector< Real > &dintnl, const std::vector< bool > &active, std::vector< bool > &deactivated_due_to_ld)
Performs one Newton-Raphson step.
void checkJacobian(const RankFourTensor &E_inv, const std::vector< Real > &intnl_old)
Checks the full Jacobian, which is just certain linear combinations of the dyieldFunction_dstress, etc, by using finite difference approximations.