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