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