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 "MultiParameterPlasticityStressUpdate.h"
11 : #include "Conversion.h" // for stringify
12 : #include "MooseEnum.h" // for enum
13 :
14 : // libMesh includes
15 : #include "libmesh/utility.h" // for Utility::pow
16 :
17 : // PETSc includes
18 : #include <petscblaslapack.h> // LAPACKgesv_
19 :
20 : InputParameters
21 2410 : MultiParameterPlasticityStressUpdate::validParams()
22 : {
23 2410 : InputParameters params = StressUpdateBase::validParams();
24 2410 : params.addClassDescription("Return-map and Jacobian algorithms for plastic models where the "
25 : "yield function and flow directions depend on multiple functions of "
26 : "stress");
27 7230 : params.addRangeCheckedParam<unsigned int>(
28 : "max_NR_iterations",
29 4820 : 20,
30 : "max_NR_iterations>0",
31 : "Maximum number of Newton-Raphson iterations allowed during the return-map algorithm");
32 4820 : params.addParam<bool>("perform_finite_strain_rotations",
33 4820 : false,
34 : "Tensors are correctly rotated "
35 : "in finite-strain simulations. "
36 : "For optimal performance you can "
37 : "set this to 'false' if you are "
38 : "only ever using small strains");
39 4820 : params.addRequiredParam<Real>(
40 : "smoothing_tol",
41 : "Intersections of the yield surfaces will be smoothed by this amount (this "
42 : "is measured in units of stress). Often this is related to other physical "
43 : "parameters (eg, 0.1*cohesion) but it is important to set this small enough "
44 : "so that the individual yield surfaces do not mix together in the smoothing "
45 : "process to produce a result where no stress is admissible (for example, "
46 : "mixing together tensile and compressive failure envelopes).");
47 4820 : params.addRequiredParam<Real>("yield_function_tol",
48 : "The return-map process will be deemed to have converged if all "
49 : "yield functions are within yield_function_tol of zero. If this "
50 : "is set very low then precision-loss might be encountered: if the "
51 : "code detects precision loss then it also deems the return-map "
52 : "process has converged.");
53 4820 : MooseEnum tangent_operator("elastic nonlinear", "nonlinear");
54 4820 : params.addParam<Real>("min_step_size",
55 4820 : 1.0,
56 : "In order to help the Newton-Raphson procedure, the applied strain "
57 : "increment may be applied in sub-increments of size greater than this "
58 : "value. Usually it is better for Moose's nonlinear convergence to "
59 : "increase max_NR_iterations rather than decrease this parameter.");
60 4820 : params.addParam<bool>("warn_about_precision_loss",
61 4820 : false,
62 : "Output a message to the console every "
63 : "time precision-loss is encountered "
64 : "during the Newton-Raphson process");
65 4820 : params.addParam<std::vector<Real>>("admissible_stress",
66 : "A single admissible value of the value of the stress "
67 : "parameters for internal parameters = 0. This is used "
68 : "to initialize the return-mapping algorithm during the first "
69 : "nonlinear iteration. If not given then it is assumed that "
70 : "stress parameters = 0 is admissible.");
71 4820 : MooseEnum smoother_fcn_enum("cos poly1 poly2 poly3", "cos");
72 4820 : params.addParam<MooseEnum>("smoother_function_type",
73 : smoother_fcn_enum,
74 : "Type of smoother function to use. 'cos' means (-a/pi)cos(pi x/2/a), "
75 : "'polyN' means a polynomial of degree 2N+2");
76 4820 : params.addParamNamesToGroup("smoother_function_type", "Advanced");
77 2410 : return params;
78 2410 : }
79 :
80 1809 : MultiParameterPlasticityStressUpdate::MultiParameterPlasticityStressUpdate(
81 1809 : const InputParameters & parameters, unsigned num_sp, unsigned num_yf, unsigned num_intnl)
82 : : StressUpdateBase(parameters),
83 1809 : _num_sp(num_sp),
84 1809 : _definitely_ok_sp(isParamValid("admissible_stress")
85 1809 : ? getParam<std::vector<Real>>("admissible_stress")
86 1809 : : std::vector<Real>(_num_sp, 0.0)),
87 1809 : _Eij(num_sp, std::vector<Real>(num_sp)),
88 1809 : _En(1.0),
89 1809 : _Cij(num_sp, std::vector<Real>(num_sp)),
90 1809 : _num_yf(num_yf),
91 1809 : _num_intnl(num_intnl),
92 3618 : _max_nr_its(getParam<unsigned>("max_NR_iterations")),
93 3618 : _perform_finite_strain_rotations(getParam<bool>("perform_finite_strain_rotations")),
94 3618 : _smoothing_tol(getParam<Real>("smoothing_tol")),
95 3618 : _smoothing_tol2(Utility::pow<2>(getParam<Real>("smoothing_tol"))),
96 3618 : _f_tol(getParam<Real>("yield_function_tol")),
97 3618 : _f_tol2(Utility::pow<2>(getParam<Real>("yield_function_tol"))),
98 3618 : _min_step_size(getParam<Real>("min_step_size")),
99 3618 : _step_one(declareRestartableData<bool>("step_one", true)),
100 3618 : _warn_about_precision_loss(getParam<bool>("warn_about_precision_loss")),
101 :
102 1809 : _plastic_strain(declareProperty<RankTwoTensor>(_base_name + "plastic_strain")),
103 3618 : _plastic_strain_old(getMaterialPropertyOld<RankTwoTensor>(_base_name + "plastic_strain")),
104 1809 : _intnl(declareProperty<std::vector<Real>>(_base_name + "plastic_internal_parameter")),
105 1809 : _intnl_old(
106 1809 : getMaterialPropertyOld<std::vector<Real>>(_base_name + "plastic_internal_parameter")),
107 1809 : _yf(declareProperty<std::vector<Real>>(_base_name + "plastic_yield_function")),
108 1809 : _iter(declareProperty<Real>(_base_name +
109 : "plastic_NR_iterations")), // this is really an unsigned int, but
110 : // for visualisation i convert it to Real
111 1809 : _max_iter_used(declareProperty<Real>(
112 1809 : _base_name + "max_plastic_NR_iterations")), // this is really an unsigned int, but
113 : // for visualisation i convert it to Real
114 3618 : _max_iter_used_old(getMaterialPropertyOld<Real>(_base_name + "max_plastic_NR_iterations")),
115 1809 : _linesearch_needed(
116 1809 : declareProperty<Real>(_base_name + "plastic_linesearch_needed")), // this is really a
117 : // boolean, but for
118 : // visualisation i
119 : // convert it to Real
120 1809 : _trial_sp(num_sp),
121 1809 : _stress_trial(RankTwoTensor()),
122 1809 : _rhs(num_sp + 1),
123 1809 : _dvar_dtrial(num_sp + 1, std::vector<Real>(num_sp, 0.0)),
124 1809 : _ok_sp(num_sp),
125 1809 : _ok_intnl(num_intnl),
126 1809 : _del_stress_params(num_sp),
127 1809 : _current_sp(num_sp),
128 1809 : _current_intnl(num_intnl),
129 1809 : _smoothed_q(yieldAndFlow(yieldAndFlow(num_sp, num_intnl))),
130 1809 : _dsp_scratch(num_sp),
131 1809 : _dsp_trial_scratch(num_sp),
132 1809 : _d2sp_scratch(num_sp),
133 1809 : _yfs_scratch(num_yf),
134 1809 : _sp_params_old_scratch(num_sp),
135 1809 : _delta_nr_params_scratch(num_sp + 1),
136 1809 : _dintnl_scratch(num_intnl, std::vector<Real>(num_sp)),
137 1809 : _jac_scratch((num_sp + 1) * (num_sp + 1)),
138 1809 : _ipiv_scratch(num_sp + 1),
139 1809 : _all_q_scratch(num_yf, yieldAndFlow(num_sp, num_intnl)),
140 :
141 1809 : _smoother_function_type(
142 3618 : parameters.get<MooseEnum>("smoother_function_type").getEnum<SmootherFunctionType>())
143 : {
144 1809 : if (_definitely_ok_sp.size() != _num_sp)
145 0 : mooseError("MultiParameterPlasticityStressUpdate: admissible_stress parameter must consist of ",
146 0 : _num_sp,
147 : " numbers");
148 1809 : }
149 :
150 : void
151 40928 : MultiParameterPlasticityStressUpdate::initQpStatefulProperties()
152 : {
153 40928 : _plastic_strain[_qp].zero();
154 40928 : _intnl[_qp].assign(_num_intnl, 0);
155 40928 : _yf[_qp].assign(_num_yf, 0);
156 40928 : _iter[_qp] = 0.0;
157 40928 : _max_iter_used[_qp] = 0.0;
158 40928 : _linesearch_needed[_qp] = 0.0;
159 40928 : }
160 :
161 : void
162 168 : MultiParameterPlasticityStressUpdate::propagateQpStatefulProperties()
163 : {
164 168 : _plastic_strain[_qp] = _plastic_strain_old[_qp];
165 168 : std::copy(_intnl_old[_qp].begin(), _intnl_old[_qp].end(), _intnl[_qp].begin());
166 168 : _max_iter_used[_qp] = _max_iter_used_old[_qp];
167 168 : }
168 :
169 : void
170 1113907 : MultiParameterPlasticityStressUpdate::updateState(RankTwoTensor & strain_increment,
171 : RankTwoTensor & inelastic_strain_increment,
172 : const RankTwoTensor & rotation_increment,
173 : RankTwoTensor & stress_new,
174 : const RankTwoTensor & stress_old,
175 : const RankFourTensor & elasticity_tensor,
176 : const RankTwoTensor & /*elastic_strain_old*/,
177 : bool compute_full_tangent_operator,
178 : RankFourTensor & tangent_operator)
179 : {
180 : // Size _yf[_qp] appropriately
181 1113907 : _yf[_qp].assign(_num_yf, 0);
182 : // _plastic_strain and _intnl are usually sized appropriately because they are stateful, but this
183 : // Material may be used from a DiracKernel where stateful materials are not allowed. The best we
184 : // can do is:
185 1113907 : if (_intnl[_qp].size() != _num_intnl)
186 0 : initQpStatefulProperties();
187 :
188 1113907 : initializeReturnProcess();
189 :
190 1113907 : if (_t_step >= 2)
191 816664 : _step_one = false;
192 :
193 : // initially assume an elastic deformation
194 1113907 : std::copy(_intnl_old[_qp].begin(), _intnl_old[_qp].end(), _intnl[_qp].begin());
195 :
196 1113907 : _iter[_qp] = 0.0;
197 1113907 : _max_iter_used[_qp] = std::max(_max_iter_used[_qp], _max_iter_used_old[_qp]);
198 1113907 : _linesearch_needed[_qp] = 0.0;
199 :
200 1113907 : computeStressParams(stress_new, _trial_sp);
201 1113907 : yieldFunctionValuesV(_trial_sp, _intnl[_qp], _yf[_qp]);
202 :
203 1113907 : if (yieldF(_yf[_qp]) <= _f_tol)
204 : {
205 624764 : _plastic_strain[_qp] = _plastic_strain_old[_qp];
206 : inelastic_strain_increment.zero();
207 624764 : if (_fe_problem.currentlyComputingJacobian())
208 153888 : tangent_operator = elasticity_tensor;
209 624764 : return;
210 : }
211 :
212 489143 : _stress_trial = stress_new;
213 : /* The trial stress must be inadmissible
214 : * so we need to return to the yield surface. The following
215 : * equations must be satisfied.
216 : *
217 : * 0 = rhs[0] = S[0] - S[0]^trial + ga * E[0, i] * dg/dS[i]
218 : * 0 = rhs[1] = S[1] - S[1]^trial + ga * E[1, i] * dg/dS[i]
219 : * ...
220 : * 0 = rhs[N-1] = S[N-1] - S[N-1]^trial + ga * E[N-1, i] * dg/dS[i]
221 : * 0 = rhs[N] = f(S, intnl)
222 : *
223 : * as well as equations defining intnl parameters as functions of
224 : * stress_params, trial_stress_params and intnl_old
225 : *
226 : * The unknowns are S[0], ..., S[N-1], gaE, and the intnl parameters.
227 : * Here gaE = ga * _En (the _En serves to make gaE similar magnitude to S)
228 : * I find it convenient to solve the first N+1 equations for p, q and gaE,
229 : * while substituting the "intnl parameters" equations into these during the solve process
230 : */
231 :
232 2053268 : for (auto & deriv : _dvar_dtrial)
233 1564125 : deriv.assign(_num_sp, 0.0);
234 :
235 489143 : preReturnMapV(_trial_sp, stress_new, _intnl_old[_qp], _yf[_qp], elasticity_tensor);
236 :
237 489143 : setEffectiveElasticity(elasticity_tensor);
238 :
239 489143 : if (_step_one)
240 : std::copy(_definitely_ok_sp.begin(), _definitely_ok_sp.end(), _ok_sp.begin());
241 : else
242 277020 : computeStressParams(stress_old, _ok_sp);
243 489143 : std::copy(_intnl_old[_qp].begin(), _intnl_old[_qp].end(), _ok_intnl.begin());
244 :
245 : // Return-map problem: must apply the following changes in stress_params,
246 : // and find the returned stress_params and gaE
247 1564125 : for (unsigned i = 0; i < _num_sp; ++i)
248 1074982 : _del_stress_params[i] = _trial_sp[i] - _ok_sp[i];
249 :
250 : Real step_taken = 0.0; // amount of del_stress_params that we've applied and the return-map
251 : // problem has succeeded
252 : Real step_size = 1.0; // potentially can apply del_stress_params in substeps
253 : Real gaE_total = 0.0;
254 :
255 : // In the following sub-stepping procedure it is possible that
256 : // the last step is an elastic step, and therefore _smoothed_q won't
257 : // be computed on the last step, so we have to compute it.
258 : bool smoothed_q_calculated = false;
259 :
260 978286 : while (step_taken < 1.0 && step_size >= _min_step_size)
261 : {
262 489143 : if (1.0 - step_taken < step_size)
263 : // prevent over-shoots of substepping
264 : step_size = 1.0 - step_taken;
265 :
266 : // trial variables in terms of admissible variables
267 1564125 : for (unsigned i = 0; i < _num_sp; ++i)
268 1074982 : _trial_sp[i] = _ok_sp[i] + step_size * _del_stress_params[i];
269 :
270 : // initialize variables that are to be found via Newton-Raphson
271 489143 : _current_sp = _trial_sp;
272 489143 : Real gaE = 0.0;
273 :
274 : // flags indicating failure of Newton-Raphson and line-search
275 : int nr_failure = 0;
276 : int ls_failure = 0;
277 :
278 : // NR iterations taken in this substep
279 : unsigned step_iter = 0;
280 :
281 : // The residual-squared for the line-search
282 489143 : Real res2 = 0.0;
283 :
284 489143 : if (step_size < 1.0 && yieldF(_trial_sp, _ok_intnl) <= _f_tol)
285 : // This is an elastic step
286 : // The "step_size < 1.0" in above condition is for efficiency: we definitely
287 : // know that this is a plastic step if step_size = 1.0
288 : smoothed_q_calculated = false;
289 : else
290 : {
291 : // this is a plastic step
292 :
293 : // initialize current_sp, gaE and current_intnl based on the non-smoothed situation
294 489143 : initializeVarsV(_trial_sp, _ok_intnl, _current_sp, gaE, _current_intnl);
295 : // and find the smoothed yield function, flow potential and derivatives
296 489143 : _smoothed_q = smoothAllQuantities(_current_sp, _current_intnl);
297 : smoothed_q_calculated = true;
298 489143 : calculateRHS(_trial_sp, _current_sp, gaE, _smoothed_q, _rhs);
299 489143 : res2 = calculateRes2(_rhs);
300 :
301 : // Perform a Newton-Raphson with linesearch to get current_sp, gaE, and also _smoothed_q
302 1756452 : while (res2 > _f_tol2 && step_iter < _max_nr_its && nr_failure == 0 && ls_failure == 0)
303 : {
304 : // solve the linear system and store the answer (the "updates") in rhs
305 1267715 : nr_failure = nrStep(_smoothed_q, _trial_sp, _current_sp, _current_intnl, gaE, _rhs);
306 1267715 : if (nr_failure != 0)
307 : break;
308 :
309 : // handle precision loss
310 1267715 : if (precisionLoss(_rhs, _current_sp, gaE))
311 : {
312 406 : if (_warn_about_precision_loss)
313 : {
314 : Moose::err << "MultiParameterPlasticityStressUpdate: precision-loss in element "
315 0 : << _current_elem->id() << " quadpoint=" << _qp << ". Stress_params =";
316 0 : for (unsigned i = 0; i < _num_sp; ++i)
317 0 : Moose::err << " " << _current_sp[i];
318 : Moose::err << " gaE = " << gaE << "\n";
319 : }
320 406 : res2 = 0.0;
321 406 : break;
322 : }
323 :
324 : // apply (parts of) the updates, re-calculate _smoothed_q, and res2
325 1267309 : ls_failure = lineSearch(res2,
326 : _current_sp,
327 : gaE,
328 : _trial_sp,
329 : _smoothed_q,
330 1267309 : _ok_intnl,
331 : _current_intnl,
332 : _rhs,
333 1267309 : _linesearch_needed[_qp]);
334 1267309 : step_iter++;
335 : }
336 : }
337 489143 : if (res2 <= _f_tol2 && step_iter < _max_nr_its && nr_failure == 0 && ls_failure == 0 &&
338 489060 : gaE >= 0.0)
339 : {
340 : // this Newton-Raphson worked fine, or this was an elastic step
341 : std::copy(_current_sp.begin(), _current_sp.end(), _ok_sp.begin());
342 489060 : gaE_total += gaE;
343 489060 : step_taken += step_size;
344 489060 : setIntnlValuesV(_trial_sp, _ok_sp, _ok_intnl, _intnl[_qp]);
345 489060 : std::copy(_intnl[_qp].begin(), _intnl[_qp].end(), _ok_intnl.begin());
346 : // calculate dp/dp_trial, dp/dq_trial, etc, for Jacobian
347 489060 : dVardTrial(!smoothed_q_calculated,
348 : _trial_sp,
349 : _ok_sp,
350 : gaE,
351 : _ok_intnl,
352 489060 : _smoothed_q,
353 : step_size,
354 : compute_full_tangent_operator,
355 489060 : _dvar_dtrial);
356 489060 : if (static_cast<Real>(step_iter) > _iter[_qp])
357 338732 : _iter[_qp] = static_cast<Real>(step_iter);
358 489060 : if (static_cast<Real>(step_iter) > _max_iter_used[_qp])
359 16372 : _max_iter_used[_qp] = static_cast<Real>(step_iter);
360 489060 : step_size *= 1.1;
361 : }
362 : else
363 : {
364 : // Newton-Raphson + line-search process failed
365 83 : step_size *= 0.5;
366 : }
367 : }
368 :
369 489143 : if (step_size < _min_step_size)
370 83 : errorHandler("MultiParameterPlasticityStressUpdate: Minimum step-size violated");
371 :
372 : // success!
373 489060 : finalizeReturnProcess(rotation_increment);
374 489060 : yieldFunctionValuesV(_ok_sp, _intnl[_qp], _yf[_qp]);
375 :
376 489060 : if (!smoothed_q_calculated)
377 0 : _smoothed_q = smoothAllQuantities(_ok_sp, _intnl[_qp]);
378 :
379 489060 : setStressAfterReturnV(
380 489060 : _stress_trial, _ok_sp, gaE_total, _intnl[_qp], _smoothed_q, elasticity_tensor, stress_new);
381 :
382 489060 : setInelasticStrainIncrementAfterReturn(_stress_trial,
383 : gaE_total,
384 : _smoothed_q,
385 : elasticity_tensor,
386 : stress_new,
387 : inelastic_strain_increment);
388 :
389 489060 : strain_increment = strain_increment - inelastic_strain_increment;
390 489060 : _plastic_strain[_qp] = _plastic_strain_old[_qp] + inelastic_strain_increment;
391 :
392 489060 : if (_fe_problem.currentlyComputingJacobian())
393 : // for efficiency, do not compute the tangent operator if not currently computing Jacobian
394 74832 : consistentTangentOperatorV(_stress_trial,
395 : _trial_sp,
396 : stress_new,
397 : _ok_sp,
398 : gaE_total,
399 : _smoothed_q,
400 : elasticity_tensor,
401 : compute_full_tangent_operator,
402 74832 : _dvar_dtrial,
403 : tangent_operator);
404 : }
405 :
406 : MultiParameterPlasticityStressUpdate::yieldAndFlow
407 2026792 : MultiParameterPlasticityStressUpdate::smoothAllQuantities(const std::vector<Real> & stress_params,
408 : const std::vector<Real> & intnl) const
409 : {
410 : mooseAssert(_all_q_scratch.size() == _num_yf,
411 : "_all_q_scratch incorrectly sized in smoothAllQuantities");
412 14586880 : for (auto & q : _all_q_scratch)
413 12560088 : q.reset();
414 2026792 : computeAllQV(stress_params, intnl, _all_q_scratch);
415 :
416 : /* This routine holds the key to my smoothing strategy. It
417 : * may be proved that this smoothing strategy produces a
418 : * yield surface that is both C2 differentiable and convex,
419 : * assuming the individual yield functions are C2 and
420 : * convex too.
421 : * Of course all the derivatives must also be smoothed.
422 : * Also, I assume that d(flow potential)/dstress gets smoothed
423 : * by the Yield Function (which produces a C2 flow potential).
424 : * See the line identified in the loop below.
425 : * Only time will tell whether this is a good strategy, but it
426 : * works well in all tests so far. Convexity is irrelevant
427 : * for the non-associated case, but at least the return-map
428 : * problem should always have a unique solution.
429 : * For two yield functions+flows, labelled 1 and 2, we
430 : * should have
431 : * d(g1 - g2) . d(f1 - f2) >= 0
432 : * If not then the return-map problem for even the
433 : * multi-surface plasticity with no smoothing won't have a
434 : * unique solution. If the multi-surface plasticity has
435 : * a unique solution then the smoothed version defined
436 : * below will too.
437 : */
438 :
439 : // res_f is the index that contains the smoothed yieldAndFlow
440 : std::size_t res_f = 0;
441 :
442 12560088 : for (std::size_t a = 1; a < _all_q_scratch.size(); ++a)
443 : {
444 10533296 : if (_all_q_scratch[res_f].f >= _all_q_scratch[a].f + _smoothing_tol)
445 : // no smoothing is needed: res_f is already indexes the largest yield function
446 8220530 : continue;
447 2312766 : else if (_all_q_scratch[a].f >= _all_q_scratch[res_f].f + _smoothing_tol)
448 : {
449 : // no smoothing is needed, and res_f needs to index to _all_q_scratch[a]
450 : res_f = a;
451 849426 : continue;
452 : }
453 : else
454 : {
455 : // smoothing is required
456 1463340 : const Real f_diff = _all_q_scratch[res_f].f - _all_q_scratch[a].f;
457 1463340 : const Real ism = ismoother(f_diff);
458 1463340 : const Real sm = smoother(f_diff);
459 1463340 : const Real dsm = dsmoother(f_diff);
460 : // we want: _all_q_scratch[res_f].f = 0.5 * _all_q_scratch[res_f].f + _all_q_scratch[a].f +
461 : // _smoothing_tol) + ism, but we have to do the derivatives first
462 5243692 : for (unsigned i = 0; i < _num_sp; ++i)
463 : {
464 13902072 : for (unsigned j = 0; j < _num_sp; ++j)
465 10121720 : _all_q_scratch[res_f].d2g[i][j] =
466 10121720 : 0.5 * (_all_q_scratch[res_f].d2g[i][j] + _all_q_scratch[a].d2g[i][j]) +
467 10121720 : dsm * (_all_q_scratch[res_f].df[j] - _all_q_scratch[a].df[j]) *
468 10121720 : (_all_q_scratch[res_f].dg[i] - _all_q_scratch[a].dg[i]) +
469 10121720 : sm * (_all_q_scratch[res_f].d2g[i][j] - _all_q_scratch[a].d2g[i][j]);
470 11108328 : for (unsigned j = 0; j < _num_intnl; ++j)
471 7327976 : _all_q_scratch[res_f].d2g_di[i][j] =
472 7327976 : 0.5 * (_all_q_scratch[res_f].d2g_di[i][j] + _all_q_scratch[a].d2g_di[i][j]) +
473 7327976 : dsm * (_all_q_scratch[res_f].df_di[j] - _all_q_scratch[a].df_di[j]) *
474 7327976 : (_all_q_scratch[res_f].dg[i] - _all_q_scratch[a].dg[i]) +
475 7327976 : sm * (_all_q_scratch[res_f].d2g_di[i][j] - _all_q_scratch[a].d2g_di[i][j]);
476 : }
477 5243692 : for (unsigned i = 0; i < _num_sp; ++i)
478 : {
479 3780352 : _all_q_scratch[res_f].df[i] =
480 3780352 : 0.5 * (_all_q_scratch[res_f].df[i] + _all_q_scratch[a].df[i]) +
481 3780352 : sm * (_all_q_scratch[res_f].df[i] - _all_q_scratch[a].df[i]);
482 : // whether the following (smoothing g with f's smoother) is a good strategy remains to be
483 : // seen...
484 3780352 : _all_q_scratch[res_f].dg[i] =
485 3780352 : 0.5 * (_all_q_scratch[res_f].dg[i] + _all_q_scratch[a].dg[i]) +
486 3780352 : sm * (_all_q_scratch[res_f].dg[i] - _all_q_scratch[a].dg[i]);
487 : }
488 4312444 : for (unsigned i = 0; i < _num_intnl; ++i)
489 2849104 : _all_q_scratch[res_f].df_di[i] =
490 2849104 : 0.5 * (_all_q_scratch[res_f].df_di[i] + _all_q_scratch[a].df_di[i]) +
491 2849104 : sm * (_all_q_scratch[res_f].df_di[i] - _all_q_scratch[a].df_di[i]);
492 1463340 : _all_q_scratch[res_f].f =
493 1463340 : 0.5 * (_all_q_scratch[res_f].f + _all_q_scratch[a].f + _smoothing_tol) + ism;
494 : }
495 : }
496 2026792 : return _all_q_scratch[res_f];
497 : }
498 :
499 : Real
500 1680308 : MultiParameterPlasticityStressUpdate::ismoother(Real f_diff) const
501 : {
502 1680308 : if (std::abs(f_diff) >= _smoothing_tol)
503 : return 0.0;
504 1680308 : switch (_smoother_function_type)
505 : {
506 1680308 : case SmootherFunctionType::cos:
507 1680308 : return -_smoothing_tol / M_PI * std::cos(0.5 * M_PI * f_diff / _smoothing_tol);
508 0 : case SmootherFunctionType::poly1:
509 0 : return 0.75 / _smoothing_tol *
510 0 : (0.5 * (Utility::pow<2>(f_diff) - _smoothing_tol2) -
511 0 : (_smoothing_tol2 / 12.0) * (Utility::pow<4>(f_diff / _smoothing_tol) - 1.0));
512 0 : case SmootherFunctionType::poly2:
513 0 : return 0.625 / _smoothing_tol *
514 0 : (0.5 * (Utility::pow<2>(f_diff) - _smoothing_tol2) -
515 0 : (_smoothing_tol2 / 30.0) * (Utility::pow<6>(f_diff / _smoothing_tol) - 1.0));
516 0 : case SmootherFunctionType::poly3:
517 0 : return (7.0 / 12.0 / _smoothing_tol) *
518 0 : (0.5 * (Utility::pow<2>(f_diff) - _smoothing_tol2) -
519 0 : (_smoothing_tol2 / 56.0) * (Utility::pow<8>(f_diff / _smoothing_tol) - 1.0));
520 : default:
521 : return 0.0;
522 : }
523 : }
524 :
525 : Real
526 1463340 : MultiParameterPlasticityStressUpdate::smoother(Real f_diff) const
527 : {
528 1463340 : if (std::abs(f_diff) >= _smoothing_tol)
529 : return 0.0;
530 1463340 : switch (_smoother_function_type)
531 : {
532 1463340 : case SmootherFunctionType::cos:
533 1463340 : return 0.5 * std::sin(f_diff * M_PI * 0.5 / _smoothing_tol);
534 0 : case SmootherFunctionType::poly1:
535 0 : return 0.75 / _smoothing_tol *
536 0 : (f_diff - (_smoothing_tol / 3.0) * Utility::pow<3>(f_diff / _smoothing_tol));
537 0 : case SmootherFunctionType::poly2:
538 0 : return 0.625 / _smoothing_tol *
539 0 : (f_diff - (_smoothing_tol / 5.0) * Utility::pow<5>(f_diff / _smoothing_tol));
540 0 : case SmootherFunctionType::poly3:
541 0 : return (7.0 / 12.0 / _smoothing_tol) *
542 0 : (f_diff - (_smoothing_tol / 7.0) * Utility::pow<7>(f_diff / _smoothing_tol));
543 : default:
544 : return 0.0;
545 : }
546 : }
547 :
548 : Real
549 1463340 : MultiParameterPlasticityStressUpdate::dsmoother(Real f_diff) const
550 : {
551 1463340 : if (std::abs(f_diff) >= _smoothing_tol)
552 : return 0.0;
553 1463340 : switch (_smoother_function_type)
554 : {
555 1463340 : case SmootherFunctionType::cos:
556 1463340 : return 0.25 * M_PI / _smoothing_tol * std::cos(f_diff * M_PI * 0.5 / _smoothing_tol);
557 0 : case SmootherFunctionType::poly1:
558 0 : return 0.75 / _smoothing_tol * (1.0 - Utility::pow<2>(f_diff / _smoothing_tol));
559 0 : case SmootherFunctionType::poly2:
560 0 : return 0.625 / _smoothing_tol * (1.0 - Utility::pow<4>(f_diff / _smoothing_tol));
561 0 : case SmootherFunctionType::poly3:
562 0 : return (7.0 / 12.0 / _smoothing_tol) * (1.0 - Utility::pow<6>(f_diff / _smoothing_tol));
563 : default:
564 : return 0.0;
565 : }
566 : }
567 :
568 : int
569 1267309 : MultiParameterPlasticityStressUpdate::lineSearch(Real & res2,
570 : std::vector<Real> & stress_params,
571 : Real & gaE,
572 : const std::vector<Real> & trial_stress_params,
573 : yieldAndFlow & smoothed_q,
574 : const std::vector<Real> & intnl_ok,
575 : std::vector<Real> & intnl,
576 : std::vector<Real> & rhs,
577 : Real & linesearch_needed) const
578 : {
579 1267309 : const Real res2_old = res2;
580 : mooseAssert(_sp_params_old_scratch.size() == _num_sp,
581 : "_sp_params_old_scratch incorrectly sized in lineSearch");
582 1267309 : _sp_params_old_scratch = stress_params;
583 1267309 : const Real gaE_old = gaE;
584 : mooseAssert(_delta_nr_params_scratch.size() == rhs.size(),
585 : "_delta_nr_params_scratch incorrectly sized in lineSearch");
586 1267309 : _delta_nr_params_scratch = rhs;
587 :
588 : Real lam = 1.0; // line-search parameter
589 : const Real lam_min = 1E-10; // minimum value of lam allowed
590 1267309 : const Real slope = -2.0 * res2_old; // "Numerical Recipes" uses -b*A*x, in order to check for
591 : // roundoff, but i hope the nrStep would warn if there were
592 : // problems
593 : Real tmp_lam; // cached value of lam used in quadratic & cubic line search
594 : Real f2 = res2_old; // cached value of f = residual2 used in the cubic in the line search
595 : Real lam2 = lam; // cached value of lam used in the cubic in the line search
596 :
597 : while (true)
598 : {
599 : // update variables using the current line-search parameter
600 5285715 : for (unsigned i = 0; i < _num_sp; ++i)
601 3748066 : stress_params[i] = _sp_params_old_scratch[i] + lam * _delta_nr_params_scratch[i];
602 1537649 : gaE = gaE_old + lam * _delta_nr_params_scratch[_num_sp];
603 :
604 : // and internal parameters
605 1537649 : setIntnlValuesV(trial_stress_params, stress_params, intnl_ok, intnl);
606 :
607 1537649 : smoothed_q = smoothAllQuantities(stress_params, intnl);
608 :
609 : // update rhs for next-time through
610 1537649 : calculateRHS(trial_stress_params, stress_params, gaE, smoothed_q, rhs);
611 1537649 : res2 = calculateRes2(rhs);
612 :
613 : // do the line-search
614 1537649 : if (res2 < res2_old + 1E-4 * lam * slope)
615 : break;
616 270422 : else if (lam < lam_min)
617 : return 1;
618 270340 : else if (lam == 1.0)
619 : {
620 : // model as a quadratic
621 163006 : tmp_lam = -0.5 * slope / (res2 - res2_old - slope);
622 : }
623 : else
624 : {
625 : // model as a cubic
626 107334 : const Real rhs1 = res2 - res2_old - lam * slope;
627 107334 : const Real rhs2 = f2 - res2_old - lam2 * slope;
628 107334 : const Real a = (rhs1 / Utility::pow<2>(lam) - rhs2 / Utility::pow<2>(lam2)) / (lam - lam2);
629 : const Real b =
630 107334 : (-lam2 * rhs1 / Utility::pow<2>(lam) + lam * rhs2 / Utility::pow<2>(lam2)) / (lam - lam2);
631 107334 : if (a == 0.0)
632 0 : tmp_lam = -slope / (2.0 * b);
633 : else
634 : {
635 107334 : const Real disc = Utility::pow<2>(b) - 3.0 * a * slope;
636 107334 : if (disc < 0)
637 0 : tmp_lam = 0.5 * lam;
638 107334 : else if (b <= 0)
639 3772 : tmp_lam = (-b + std::sqrt(disc)) / (3.0 * a);
640 : else
641 103562 : tmp_lam = -slope / (b + std::sqrt(disc));
642 : }
643 107334 : if (tmp_lam > 0.5 * lam)
644 2300 : tmp_lam = 0.5 * lam;
645 : }
646 : lam2 = lam;
647 : f2 = res2;
648 270340 : lam = std::max(tmp_lam, 0.1 * lam);
649 270340 : }
650 :
651 1267227 : if (lam < 1.0)
652 162924 : linesearch_needed = 1.0;
653 : return 0;
654 : }
655 :
656 : int
657 1267715 : MultiParameterPlasticityStressUpdate::nrStep(const yieldAndFlow & smoothed_q,
658 : const std::vector<Real> & trial_stress_params,
659 : const std::vector<Real> & stress_params,
660 : const std::vector<Real> & intnl,
661 : Real gaE,
662 : std::vector<Real> & rhs) const
663 : {
664 : mooseAssert(_dintnl_scratch.size() == _num_intnl, "_dintnl_scratch incorrectly sized in nrStep");
665 1267715 : setIntnlDerivativesV(trial_stress_params, stress_params, intnl, _dintnl_scratch);
666 :
667 : mooseAssert(_jac_scratch.size() == (_num_sp + 1) * (_num_sp + 1),
668 : "_jac_scratch incorrectly sized in nrStep");
669 1267715 : dnRHSdVar(smoothed_q, _dintnl_scratch, stress_params, gaE, _jac_scratch);
670 :
671 : // use LAPACK to solve the linear system
672 1267715 : const PetscBLASInt nrhs = 1;
673 : mooseAssert(_ipiv_scratch.size() == _num_sp + 1, "_ipiv_scratch incorrectly sized in nrStep");
674 : PetscBLASInt info;
675 1267715 : const PetscBLASInt gesv_num_rhs = _num_sp + 1;
676 1267715 : LAPACKgesv_(&gesv_num_rhs,
677 : &nrhs,
678 : &_jac_scratch[0],
679 : &gesv_num_rhs,
680 : &_ipiv_scratch[0],
681 : &rhs[0],
682 : &gesv_num_rhs,
683 : &info);
684 1267715 : return info;
685 : }
686 :
687 : void
688 83 : MultiParameterPlasticityStressUpdate::errorHandler(const std::string & message) const
689 : {
690 166 : throw MooseException(message);
691 : }
692 :
693 : void
694 108472 : MultiParameterPlasticityStressUpdate::initializeReturnProcess()
695 : {
696 108472 : }
697 :
698 : void
699 96696 : MultiParameterPlasticityStressUpdate::finalizeReturnProcess(
700 : const RankTwoTensor & /*rotation_increment*/)
701 : {
702 96696 : }
703 :
704 : void
705 0 : MultiParameterPlasticityStressUpdate::preReturnMapV(
706 : const std::vector<Real> & /*trial_stress_params*/,
707 : const RankTwoTensor & /*stress_trial*/,
708 : const std::vector<Real> & /*intnl_old*/,
709 : const std::vector<Real> & /*yf*/,
710 : const RankFourTensor & /*Eijkl*/)
711 : {
712 0 : }
713 :
714 : Real
715 0 : MultiParameterPlasticityStressUpdate::yieldF(const std::vector<Real> & stress_params,
716 : const std::vector<Real> & intnl) const
717 : {
718 : mooseAssert(_yfs_scratch.size() == _num_yf, "_yfs_scratch incorrectly sized in yieldF");
719 0 : yieldFunctionValuesV(stress_params, intnl, _yfs_scratch);
720 0 : return yieldF(_yfs_scratch);
721 : }
722 :
723 : Real
724 1113907 : MultiParameterPlasticityStressUpdate::yieldF(const std::vector<Real> & yfs) const
725 : {
726 1113907 : Real yf = yfs[0];
727 4226745 : for (std::size_t i = 1; i < yfs.size(); ++i)
728 3112838 : if (yf >= yfs[i] + _smoothing_tol)
729 : // no smoothing is needed, and yf is the biggest yield function
730 2340741 : continue;
731 772097 : else if (yfs[i] >= yf + _smoothing_tol)
732 : // no smoothing is needed, and yfs[i] is the biggest yield function
733 : yf = yfs[i];
734 : else
735 216968 : yf = 0.5 * (yf + yfs[i] + _smoothing_tol) + ismoother(yf - yfs[i]);
736 1113907 : return yf;
737 : }
738 :
739 : void
740 0 : MultiParameterPlasticityStressUpdate::initializeVarsV(const std::vector<Real> & trial_stress_params,
741 : const std::vector<Real> & intnl_old,
742 : std::vector<Real> & stress_params,
743 : Real & gaE,
744 : std::vector<Real> & intnl) const
745 : {
746 0 : gaE = 0.0;
747 : std::copy(trial_stress_params.begin(), trial_stress_params.end(), stress_params.begin());
748 : std::copy(intnl_old.begin(), intnl_old.end(), intnl.begin());
749 0 : }
750 :
751 : void
752 128 : MultiParameterPlasticityStressUpdate::consistentTangentOperatorV(
753 : const RankTwoTensor & stress_trial,
754 : const std::vector<Real> & /*trial_stress_params*/,
755 : const RankTwoTensor & stress,
756 : const std::vector<Real> & /*stress_params*/,
757 : Real gaE,
758 : const yieldAndFlow & smoothed_q,
759 : const RankFourTensor & elasticity_tensor,
760 : bool compute_full_tangent_operator,
761 : const std::vector<std::vector<Real>> & dvar_dtrial,
762 : RankFourTensor & cto)
763 : {
764 128 : cto = elasticity_tensor;
765 128 : if (!compute_full_tangent_operator)
766 0 : return;
767 :
768 128 : const Real ga = gaE / _En;
769 :
770 128 : dstressparam_dstress(stress, _dsp_scratch);
771 128 : dstressparam_dstress(stress_trial, _dsp_trial_scratch);
772 :
773 512 : for (unsigned a = 0; a < _num_sp; ++a)
774 : {
775 1536 : for (unsigned b = 0; b < _num_sp; ++b)
776 : {
777 1152 : const RankTwoTensor t = elasticity_tensor * _dsp_trial_scratch[a];
778 1152 : RankTwoTensor s = _Cij[b][a] * _dsp_scratch[b];
779 4608 : for (unsigned c = 0; c < _num_sp; ++c)
780 3456 : s -= _dsp_scratch[b] * _Cij[b][c] * dvar_dtrial[c][a];
781 1152 : s = elasticity_tensor * s;
782 1152 : cto -= s.outerProduct(t);
783 : }
784 : }
785 :
786 128 : d2stressparam_dstress(stress, _d2sp_scratch);
787 128 : RankFourTensor Tijab;
788 512 : for (unsigned i = 0; i < _num_sp; ++i)
789 768 : Tijab += smoothed_q.dg[i] * _d2sp_scratch[i];
790 128 : Tijab = ga * elasticity_tensor * Tijab;
791 :
792 128 : RankFourTensor inv = RankFourTensor(RankFourTensor::initIdentityFour) + Tijab;
793 : try
794 : {
795 128 : inv = inv.transposeMajor().invSymm();
796 : }
797 0 : catch (const MooseException & e)
798 : {
799 : // Cannot form the inverse, so probably at some degenerate place in stress space.
800 : // Just return with the "best estimate" of the cto.
801 0 : mooseWarning("MultiParameterPlasticityStressUpdate: Cannot invert 1+T in consistent tangent "
802 : "operator computation at quadpoint ",
803 0 : _qp,
804 : " of element ",
805 0 : _current_elem->id());
806 : return;
807 0 : }
808 :
809 128 : cto = (cto.transposeMajor() * inv).transposeMajor();
810 : }
811 :
812 : void
813 96696 : MultiParameterPlasticityStressUpdate::setInelasticStrainIncrementAfterReturn(
814 : const RankTwoTensor & /*stress_trial*/,
815 : Real gaE,
816 : const yieldAndFlow & smoothed_q,
817 : const RankFourTensor & /*elasticity_tensor*/,
818 : const RankTwoTensor & returned_stress,
819 : RankTwoTensor & inelastic_strain_increment)
820 : {
821 96696 : dstressparam_dstress(returned_stress, _dsp_scratch);
822 96696 : inelastic_strain_increment = RankTwoTensor();
823 386784 : for (unsigned i = 0; i < _num_sp; ++i)
824 290088 : inelastic_strain_increment += smoothed_q.dg[i] * _dsp_scratch[i];
825 96696 : inelastic_strain_increment *= gaE / _En;
826 96696 : }
827 :
828 : Real
829 2026792 : MultiParameterPlasticityStressUpdate::calculateRes2(const std::vector<Real> & rhs) const
830 : {
831 : Real res2 = 0.0;
832 8876632 : for (const auto & r : rhs)
833 6849840 : res2 += r * r;
834 2026792 : return res2;
835 : }
836 :
837 : void
838 2026792 : MultiParameterPlasticityStressUpdate::calculateRHS(const std::vector<Real> & trial_stress_params,
839 : const std::vector<Real> & stress_params,
840 : Real gaE,
841 : const yieldAndFlow & smoothed_q,
842 : std::vector<Real> & rhs) const
843 : {
844 2026792 : const Real ga = gaE / _En;
845 6849840 : for (unsigned i = 0; i < _num_sp; ++i)
846 : {
847 4823048 : rhs[i] = stress_params[i] - trial_stress_params[i];
848 16777536 : for (unsigned j = 0; j < _num_sp; ++j)
849 11954488 : rhs[i] += ga * _Eij[i][j] * smoothed_q.dg[j];
850 : }
851 2026792 : rhs[_num_sp] = smoothed_q.f;
852 2026792 : }
853 :
854 : void
855 1342547 : MultiParameterPlasticityStressUpdate::dnRHSdVar(const yieldAndFlow & smoothed_q,
856 : const std::vector<std::vector<Real>> & dintnl,
857 : const std::vector<Real> & /*stress_params*/,
858 : Real gaE,
859 : std::vector<double> & jac) const
860 : {
861 16370216 : for (auto & jac_entry : jac)
862 15027669 : jac_entry = 0.0;
863 :
864 1342547 : const Real ga = gaE / _En;
865 :
866 : unsigned ind = 0;
867 4448319 : for (unsigned var = 0; var < _num_sp; ++var)
868 : {
869 10579350 : for (unsigned rhs = 0; rhs < _num_sp; ++rhs)
870 : {
871 7473578 : if (var == rhs)
872 3105772 : jac[ind] -= 1.0;
873 26206836 : for (unsigned j = 0; j < _num_sp; ++j)
874 : {
875 18733258 : jac[ind] -= ga * _Eij[rhs][j] * smoothed_q.d2g[j][var];
876 55150014 : for (unsigned k = 0; k < _num_intnl; ++k)
877 36416756 : jac[ind] -= ga * _Eij[rhs][j] * smoothed_q.d2g_di[j][k] * dintnl[k][var];
878 : }
879 7473578 : ind++;
880 : }
881 : // now rhs = _num_sp (that is, the yield function)
882 3105772 : jac[ind] -= smoothed_q.df[var];
883 9200676 : for (unsigned k = 0; k < _num_intnl; ++k)
884 6094904 : jac[ind] -= smoothed_q.df_di[k] * dintnl[k][var];
885 3105772 : ind++;
886 : }
887 :
888 : // now var = _num_sp (that is, gaE)
889 4448319 : for (unsigned rhs = 0; rhs < _num_sp; ++rhs)
890 : {
891 10579350 : for (unsigned j = 0; j < _num_sp; ++j)
892 7473578 : jac[ind] -= (1.0 / _En) * _Eij[rhs][j] * smoothed_q.dg[j];
893 3105772 : ind++;
894 : }
895 : // now rhs = _num_sp (that is, the yield function)
896 1342547 : jac[ind] = 0.0;
897 1342547 : }
898 :
899 : void
900 489060 : MultiParameterPlasticityStressUpdate::dVardTrial(bool elastic_only,
901 : const std::vector<Real> & trial_stress_params,
902 : const std::vector<Real> & stress_params,
903 : Real gaE,
904 : const std::vector<Real> & intnl,
905 : const yieldAndFlow & smoothed_q,
906 : Real step_size,
907 : bool compute_full_tangent_operator,
908 : std::vector<std::vector<Real>> & dvar_dtrial) const
909 : {
910 489060 : if (!_fe_problem.currentlyComputingJacobian())
911 414228 : return;
912 :
913 74832 : if (!compute_full_tangent_operator)
914 : return;
915 :
916 74832 : if (elastic_only)
917 : {
918 : // no change to gaE, and all off-diag stuff remains unchanged from previous step
919 0 : for (unsigned v = 0; v < _num_sp; ++v)
920 0 : dvar_dtrial[v][v] += step_size;
921 : return;
922 : }
923 :
924 74832 : const Real ga = gaE / _En;
925 :
926 : // the following uses heap allocations in the belief that the compute time spent on alloc/dealloc
927 : // will be dwarfed by other aspects of implicit time stepping
928 :
929 74832 : std::vector<std::vector<Real>> dintnl(_num_intnl, std::vector<Real>(_num_sp));
930 74832 : setIntnlDerivativesV(trial_stress_params, stress_params, intnl, dintnl);
931 :
932 : // rhs is described elsewhere, the following are changes in rhs wrt the trial_stress_param
933 : // values
934 : // In the following we use d(intnl)/d(trial variable) = - d(intnl)/d(variable)
935 74832 : std::vector<Real> rhs_cto((_num_sp + 1) * _num_sp);
936 :
937 : unsigned ind = 0;
938 225216 : for (unsigned a = 0; a < _num_sp; ++a)
939 : {
940 : // change in RHS[b] wrt changes in stress_param_trial[a]
941 453312 : for (unsigned b = 0; b < _num_sp; ++b)
942 : {
943 302928 : if (a == b)
944 150384 : rhs_cto[ind] -= 1.0;
945 915264 : for (unsigned j = 0; j < _num_sp; ++j)
946 1830096 : for (unsigned k = 0; k < _num_intnl; ++k)
947 1217760 : rhs_cto[ind] -= ga * _Eij[b][j] * smoothed_q.d2g_di[j][k] * dintnl[k][a];
948 302928 : ind++;
949 : }
950 : // now b = _num_sp (that is, the yield function)
951 450384 : for (unsigned k = 0; k < _num_intnl; ++k)
952 300000 : rhs_cto[ind] -= smoothed_q.df_di[k] * dintnl[k][a];
953 150384 : ind++;
954 : }
955 :
956 : // jac = d(-rhs)/d(var)
957 74832 : std::vector<double> jac((_num_sp + 1) * (_num_sp + 1));
958 74832 : dnRHSdVar(smoothed_q, dintnl, stress_params, gaE, jac);
959 :
960 74832 : std::vector<PetscBLASInt> ipiv(_num_sp + 1);
961 : PetscBLASInt info;
962 74832 : const PetscBLASInt gesv_num_rhs = _num_sp + 1;
963 74832 : const PetscBLASInt gesv_num_pq = _num_sp;
964 74832 : LAPACKgesv_(&gesv_num_rhs,
965 : &gesv_num_pq,
966 : &jac[0],
967 : &gesv_num_rhs,
968 : &ipiv[0],
969 : &rhs_cto[0],
970 : &gesv_num_rhs,
971 : &info);
972 74832 : if (info != 0)
973 0 : errorHandler("MultiParameterPlasticityStressUpdate: PETSC LAPACK gsev routine returned with "
974 0 : "error code " +
975 0 : Moose::stringify(info));
976 :
977 : ind = 0;
978 74832 : std::vector<std::vector<Real>> dvarn_dtrialn(_num_sp + 1, std::vector<Real>(_num_sp, 0.0));
979 225216 : for (unsigned spt = 0; spt < _num_sp; ++spt) // loop over trial stress-param variables
980 : {
981 453312 : for (unsigned v = 0; v < _num_sp; ++v) // loop over variables in NR procedure
982 : {
983 302928 : dvarn_dtrialn[v][spt] = rhs_cto[ind];
984 302928 : ind++;
985 : }
986 : // the final NR variable is gaE
987 150384 : dvarn_dtrialn[_num_sp][spt] = rhs_cto[ind];
988 150384 : ind++;
989 : }
990 :
991 74832 : const std::vector<std::vector<Real>> dvar_dtrial_old = dvar_dtrial;
992 :
993 225216 : for (unsigned v = 0; v < _num_sp; ++v) // loop over variables in NR procedure
994 : {
995 453312 : for (unsigned spt = 0; spt < _num_sp; ++spt) // loop over trial stress-param variables
996 : {
997 302928 : dvar_dtrial[v][spt] = step_size * dvarn_dtrialn[v][spt];
998 915264 : for (unsigned a = 0; a < _num_sp; ++a)
999 612336 : dvar_dtrial[v][spt] += dvarn_dtrialn[v][a] * dvar_dtrial_old[a][spt];
1000 : }
1001 : }
1002 : // for gaE the formulae are a little different
1003 : const unsigned v = _num_sp;
1004 225216 : for (unsigned spt = 0; spt < _num_sp; ++spt)
1005 : {
1006 150384 : dvar_dtrial[v][spt] += step_size * dvarn_dtrialn[v][spt]; // note +=
1007 453312 : for (unsigned a = 0; a < _num_sp; ++a)
1008 302928 : dvar_dtrial[v][spt] += dvarn_dtrialn[v][a] * dvar_dtrial_old[a][spt];
1009 : }
1010 74832 : }
1011 :
1012 : bool
1013 1267715 : MultiParameterPlasticityStressUpdate::precisionLoss(const std::vector<Real> & solution,
1014 : const std::vector<Real> & stress_params,
1015 : Real gaE) const
1016 : {
1017 1267715 : if (std::abs(solution[_num_sp]) > 1E-13 * std::abs(gaE))
1018 : return false;
1019 3592 : for (unsigned i = 0; i < _num_sp; ++i)
1020 3186 : if (std::abs(solution[i]) > 1E-13 * std::abs(stress_params[i]))
1021 : return false;
1022 : return true;
1023 : }
|
