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 "MultiPlasticityLinearSystem.h"
11 : #include "RankFourTensor.h"
12 :
13 : // Following is for perturbing distances in eliminating linearly-dependent directions
14 : #include "MooseRandom.h"
15 :
16 : // Following is used to access PETSc's LAPACK routines
17 : #include <petscblaslapack.h>
18 :
19 : InputParameters
20 4174 : MultiPlasticityLinearSystem::validParams()
21 : {
22 4174 : InputParameters params = MultiPlasticityRawComponentAssembler::validParams();
23 12522 : params.addRangeCheckedParam<Real>("linear_dependent",
24 8348 : 1E-4,
25 : "linear_dependent>=0 & linear_dependent<1",
26 : "Flow directions are considered linearly dependent if the "
27 : "smallest singular value is less than linear_dependent times "
28 : "the largest singular value");
29 4174 : return params;
30 0 : }
31 :
32 3122 : MultiPlasticityLinearSystem::MultiPlasticityLinearSystem(const MooseObject * moose_object)
33 : : MultiPlasticityRawComponentAssembler(moose_object),
34 6244 : _svd_tol(_params.get<Real>("linear_dependent")),
35 3122 : _min_f_tol(-1.0)
36 : {
37 7942 : for (unsigned model = 0; model < _num_models; ++model)
38 4820 : if (_min_f_tol == -1.0 || _min_f_tol > _f[model]->_f_tol)
39 3050 : _min_f_tol = _f[model]->_f_tol;
40 :
41 : MooseRandom::seed(0);
42 3122 : }
43 :
44 : int
45 185952 : MultiPlasticityLinearSystem::singularValuesOfR(const std::vector<RankTwoTensor> & r,
46 : std::vector<Real> & s)
47 : {
48 185952 : PetscBLASInt bm = r.size();
49 185952 : PetscBLASInt bn = 6;
50 :
51 185952 : s.resize(std::min(bm, bn));
52 :
53 : // prepare for gesvd or gesdd routine provided by PETSc
54 : // Want to find the singular values of matrix
55 : // ( r[0](0,0) r[0](0,1) r[0](0,2) r[0](1,1) r[0](1,2) r[0](2,2) )
56 : // ( r[1](0,0) r[1](0,1) r[1](0,2) r[1](1,1) r[1](1,2) r[1](2,2) )
57 : // a = ( r[2](0,0) r[2](0,1) r[2](0,2) r[2](1,1) r[2](1,2) r[2](2,2) )
58 : // ( r[3](0,0) r[3](0,1) r[3](0,2) r[3](1,1) r[3](1,2) r[3](2,2) )
59 : // ( r[4](0,0) r[4](0,1) r[4](0,2) r[4](1,1) r[4](1,2) r[4](2,2) )
60 : // bm = 5
61 :
62 185952 : std::vector<double> a(bm * 6);
63 : // Fill in the a "matrix" by going down columns
64 : unsigned ind = 0;
65 743808 : for (int col = 0; col < 3; ++col)
66 2145456 : for (int row = 0; row < bm; ++row)
67 1587600 : a[ind++] = r[row](0, col);
68 557856 : for (int col = 3; col < 5; ++col)
69 1430304 : for (int row = 0; row < bm; ++row)
70 1058400 : a[ind++] = r[row](1, col - 2);
71 715152 : for (int row = 0; row < bm; ++row)
72 529200 : a[ind++] = r[row](2, 2);
73 :
74 : // u and vt are dummy variables because they won't
75 : // get referenced due to the "N" and "N" choices
76 185952 : PetscBLASInt sizeu = 1;
77 185952 : std::vector<double> u(sizeu);
78 185952 : PetscBLASInt sizevt = 1;
79 185952 : std::vector<double> vt(sizevt);
80 :
81 185952 : PetscBLASInt sizework =
82 185952 : 16 * (bm + 6); // this is above the lowerbound specified in the LAPACK doco
83 185952 : std::vector<double> work(sizework);
84 :
85 : PetscBLASInt info;
86 :
87 185952 : LAPACKgesvd_("N",
88 : "N",
89 : &bm,
90 : &bn,
91 : &a[0],
92 : &bm,
93 : &s[0],
94 : &u[0],
95 : &sizeu,
96 : &vt[0],
97 : &sizevt,
98 : &work[0],
99 : &sizework,
100 : &info);
101 :
102 371904 : return info;
103 : }
104 :
105 : void
106 544304 : MultiPlasticityLinearSystem::eliminateLinearDependence(const RankTwoTensor & stress,
107 : const std::vector<Real> & intnl,
108 : const std::vector<Real> & f,
109 : const std::vector<RankTwoTensor> & r,
110 : const std::vector<bool> & active,
111 : std::vector<bool> & deactivated_due_to_ld)
112 : {
113 544304 : deactivated_due_to_ld.resize(_num_surfaces, false);
114 :
115 544304 : unsigned num_active = r.size();
116 :
117 544304 : if (num_active <= 1)
118 523312 : return;
119 :
120 : std::vector<double> s;
121 98832 : int info = singularValuesOfR(r, s);
122 98832 : if (info != 0)
123 0 : mooseError("In finding the SVD in the return-map algorithm, the PETSC LAPACK gesvd routine "
124 : "returned with error code ",
125 : info);
126 :
127 : // num_lin_dep are the number of linearly dependent
128 : // "r vectors", if num_active <= 6
129 : unsigned int num_lin_dep = 0;
130 :
131 98832 : unsigned i = s.size();
132 155708 : while (i-- > 0)
133 155708 : if (s[i] < _svd_tol * s[0])
134 56876 : num_lin_dep++;
135 : else
136 : break;
137 :
138 98832 : if (num_lin_dep == 0 && num_active <= 6)
139 : return;
140 :
141 : // From here on, some flow directions are linearly dependent
142 :
143 : // Find the signed "distance" of the current (stress, internal) configuration
144 : // from the yield surfaces. This distance will not be precise, but
145 : // i want to preferentially deactivate yield surfaces that are close
146 : // to the current stress point.
147 : std::vector<RankTwoTensor> df_dstress;
148 20992 : dyieldFunction_dstress(stress, intnl, active, df_dstress);
149 :
150 : typedef std::pair<Real, unsigned> pair_for_sorting;
151 20992 : std::vector<pair_for_sorting> dist(num_active);
152 129172 : for (unsigned i = 0; i < num_active; ++i)
153 : {
154 108180 : dist[i].first = f[i] / df_dstress[i].L2norm();
155 108180 : dist[i].second = i;
156 : }
157 20992 : std::sort(dist.begin(), dist.end()); // sorted in ascending order
158 :
159 : // There is a potential problem when we have equal f[i], for it can give oscillations
160 : bool equals_detected = false;
161 108180 : for (unsigned i = 0; i < num_active - 1; ++i)
162 87188 : if (std::abs(dist[i].first - dist[i + 1].first) < _min_f_tol)
163 : {
164 : equals_detected = true;
165 10032 : dist[i].first += _min_f_tol * (MooseRandom::rand() - 0.5);
166 : }
167 20992 : if (equals_detected)
168 7248 : std::sort(dist.begin(), dist.end()); // sorted in ascending order
169 :
170 : std::vector<bool> scheduled_for_deactivation;
171 : scheduled_for_deactivation.assign(num_active, false);
172 :
173 : // In the following loop we go through all the flow directions, from
174 : // the one with the largest dist, to the one with the smallest dist,
175 : // adding them one-by-one into r_tmp. Upon each addition we check
176 : // for linear-dependence. if LD is found, we schedule the most
177 : // recently added flow direction for deactivation, and pop it
178 : // back off r_tmp
179 : unsigned current_yf;
180 20992 : current_yf = dist[num_active - 1].second;
181 : // the one with largest dist
182 20992 : std::vector<RankTwoTensor> r_tmp = {r[current_yf]};
183 :
184 : unsigned num_kept_active = 1;
185 91288 : for (unsigned yf_to_try = 2; yf_to_try <= num_active; ++yf_to_try)
186 : {
187 87120 : current_yf = dist[num_active - yf_to_try].second;
188 87120 : if (num_active == 2) // shortcut to we don't have to singularValuesOfR
189 : scheduled_for_deactivation[current_yf] = true;
190 87120 : else if (num_kept_active >= 6) // shortcut to we don't have to singularValuesOfR: there can
191 : // never be > 6 linearly-independent r vectors
192 : scheduled_for_deactivation[current_yf] = true;
193 : else
194 : {
195 87120 : r_tmp.push_back(r[current_yf]);
196 87120 : info = singularValuesOfR(r_tmp, s);
197 87120 : if (info != 0)
198 0 : mooseError("In finding the SVD in the return-map algorithm, the PETSC LAPACK gesvd routine "
199 : "returned with error code ",
200 : info);
201 87120 : if (s[s.size() - 1] < _svd_tol * s[0])
202 : {
203 : scheduled_for_deactivation[current_yf] = true;
204 : r_tmp.pop_back();
205 63732 : num_lin_dep--;
206 : }
207 : else
208 23388 : num_kept_active++;
209 87120 : if (num_lin_dep == 0 && num_active <= 6)
210 : // have taken out all the vectors that were linearly dependent
211 : // so no point continuing
212 : break;
213 : }
214 : }
215 :
216 : unsigned int old_active_number = 0;
217 166112 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
218 145120 : if (active[surface])
219 : {
220 108180 : if (scheduled_for_deactivation[old_active_number])
221 : deactivated_due_to_ld[surface] = true;
222 108180 : old_active_number++;
223 : }
224 20992 : }
225 :
226 : void
227 1308900 : MultiPlasticityLinearSystem::calculateConstraints(const RankTwoTensor & stress,
228 : const std::vector<Real> & intnl_old,
229 : const std::vector<Real> & intnl,
230 : const std::vector<Real> & pm,
231 : const RankTwoTensor & delta_dp,
232 : std::vector<Real> & f,
233 : std::vector<RankTwoTensor> & r,
234 : RankTwoTensor & epp,
235 : std::vector<Real> & ic,
236 : const std::vector<bool> & active)
237 : {
238 : // see comments at the start of .h file
239 :
240 : mooseAssert(intnl_old.size() == _num_models,
241 : "Size of intnl_old is " << intnl_old.size()
242 : << " which is incorrect in calculateConstraints");
243 : mooseAssert(intnl.size() == _num_models,
244 : "Size of intnl is " << intnl.size() << " which is incorrect in calculateConstraints");
245 : mooseAssert(pm.size() == _num_surfaces,
246 : "Size of pm is " << pm.size() << " which is incorrect in calculateConstraints");
247 : mooseAssert(active.size() == _num_surfaces,
248 : "Size of active is " << active.size()
249 : << " which is incorrect in calculateConstraints");
250 :
251 : // yield functions
252 1308900 : yieldFunction(stress, intnl, active, f);
253 :
254 : // flow directions and "epp"
255 1308900 : flowPotential(stress, intnl, active, r);
256 1308900 : epp = RankTwoTensor();
257 : unsigned ind = 0;
258 3724532 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
259 2415632 : if (active[surface])
260 1746096 : epp += pm[surface] * r[ind++]; // note, even the deactivated_due_to_ld must get added in
261 1308900 : epp -= delta_dp;
262 :
263 : // internal constraints
264 : std::vector<Real> h;
265 1308900 : hardPotential(stress, intnl, active, h);
266 1308900 : ic.resize(0);
267 : ind = 0;
268 : std::vector<unsigned int> active_surfaces;
269 : std::vector<unsigned int>::iterator active_surface;
270 3429524 : for (unsigned model = 0; model < _num_models; ++model)
271 : {
272 2120624 : activeSurfaces(model, active, active_surfaces);
273 2120624 : if (active_surfaces.size() > 0)
274 : {
275 : // some surfaces are active in this model, so must form an internal constraint
276 1716144 : ic.push_back(intnl[model] - intnl_old[model]);
277 3462240 : for (active_surface = active_surfaces.begin(); active_surface != active_surfaces.end();
278 : ++active_surface)
279 1746096 : ic[ic.size() - 1] += pm[*active_surface] * h[ind++]; // we know the correct one is h[ind]
280 : // since it was constructed in the same
281 : // manner
282 : }
283 : }
284 1308900 : }
285 :
286 : void
287 544304 : MultiPlasticityLinearSystem::calculateRHS(const RankTwoTensor & stress,
288 : const std::vector<Real> & intnl_old,
289 : const std::vector<Real> & intnl,
290 : const std::vector<Real> & pm,
291 : const RankTwoTensor & delta_dp,
292 : std::vector<Real> & rhs,
293 : const std::vector<bool> & active,
294 : bool eliminate_ld,
295 : std::vector<bool> & deactivated_due_to_ld)
296 : {
297 : // see comments at the start of .h file
298 :
299 : mooseAssert(intnl_old.size() == _num_models,
300 : "Size of intnl_old is " << intnl_old.size() << " which is incorrect in calculateRHS");
301 : mooseAssert(intnl.size() == _num_models,
302 : "Size of intnl is " << intnl.size() << " which is incorrect in calculateRHS");
303 : mooseAssert(pm.size() == _num_surfaces,
304 : "Size of pm is " << pm.size() << " which is incorrect in calculateRHS");
305 : mooseAssert(active.size() == _num_surfaces,
306 : "Size of active is " << active.size() << " which is incorrect in calculateRHS");
307 :
308 : std::vector<Real> f; // the yield functions
309 544304 : RankTwoTensor epp; // the plastic-strain constraint ("direction constraint")
310 : std::vector<Real> ic; // the "internal constraints"
311 :
312 : std::vector<RankTwoTensor> r;
313 544304 : calculateConstraints(stress, intnl_old, intnl, pm, delta_dp, f, r, epp, ic, active);
314 :
315 544304 : if (eliminate_ld)
316 544304 : eliminateLinearDependence(stress, intnl, f, r, active, deactivated_due_to_ld);
317 : else
318 0 : deactivated_due_to_ld.assign(_num_surfaces, false);
319 :
320 544304 : std::vector<bool> active_not_deact(_num_surfaces);
321 1580928 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
322 1397232 : active_not_deact[surface] = (active[surface] && !deactivated_due_to_ld[surface]);
323 :
324 : unsigned num_active_f = 0;
325 1580928 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
326 1036624 : if (active_not_deact[surface])
327 676016 : num_active_f++;
328 :
329 : unsigned num_active_ic = 0;
330 1433424 : for (unsigned model = 0; model < _num_models; ++model)
331 889120 : if (anyActiveSurfaces(model, active_not_deact))
332 663504 : num_active_ic++;
333 :
334 : unsigned int dim = 3;
335 544304 : unsigned int system_size = 6 + num_active_f + num_active_ic; // "6" comes from symmeterizing epp,
336 : // num_active_f comes from "f",
337 : // num_active_f comes from "ic"
338 :
339 544304 : rhs.resize(system_size);
340 :
341 : unsigned ind = 0;
342 2177216 : for (unsigned i = 0; i < dim; ++i)
343 4898736 : for (unsigned j = 0; j <= i; ++j)
344 3265824 : rhs[ind++] = -epp(i, j);
345 : unsigned active_surface = 0;
346 1580928 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
347 1036624 : if (active[surface])
348 : {
349 739748 : if (!deactivated_due_to_ld[surface])
350 676016 : rhs[ind++] = -f[active_surface];
351 739748 : active_surface++;
352 : }
353 : unsigned active_model = 0;
354 1433424 : for (unsigned model = 0; model < _num_models; ++model)
355 889120 : if (anyActiveSurfaces(model, active))
356 : {
357 724772 : if (anyActiveSurfaces(model, active_not_deact))
358 663504 : rhs[ind++] = -ic[active_model];
359 724772 : active_model++;
360 : }
361 :
362 : mooseAssert(ind == system_size, "Incorrect filling of the rhs in calculateRHS");
363 544304 : }
364 :
365 : void
366 544304 : MultiPlasticityLinearSystem::calculateJacobian(const RankTwoTensor & stress,
367 : const std::vector<Real> & intnl,
368 : const std::vector<Real> & pm,
369 : const RankFourTensor & E_inv,
370 : const std::vector<bool> & active,
371 : const std::vector<bool> & deactivated_due_to_ld,
372 : std::vector<std::vector<Real>> & jac)
373 : {
374 : // see comments at the start of .h file
375 :
376 : mooseAssert(intnl.size() == _num_models,
377 : "Size of intnl is " << intnl.size() << " which is incorrect in calculateJacobian");
378 : mooseAssert(pm.size() == _num_surfaces,
379 : "Size of pm is " << pm.size() << " which is incorrect in calculateJacobian");
380 : mooseAssert(active.size() == _num_surfaces,
381 : "Size of active is " << active.size() << " which is incorrect in calculateJacobian");
382 : mooseAssert(deactivated_due_to_ld.size() == _num_surfaces,
383 : "Size of deactivated_due_to_ld is " << deactivated_due_to_ld.size()
384 : << " which is incorrect in calculateJacobian");
385 :
386 : unsigned ind = 0;
387 : unsigned active_surface_ind = 0;
388 :
389 544304 : std::vector<bool> active_surface(_num_surfaces); // active and not deactivated_due_to_ld
390 1580928 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
391 1397232 : active_surface[surface] = (active[surface] && !deactivated_due_to_ld[surface]);
392 : unsigned num_active_surface = 0;
393 1580928 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
394 1036624 : if (active_surface[surface])
395 676016 : num_active_surface++;
396 :
397 : std::vector<bool> active_model(
398 544304 : _num_models); // whether a model has surfaces that are active and not deactivated_due_to_ld
399 1433424 : for (unsigned model = 0; model < _num_models; ++model)
400 889120 : active_model[model] = anyActiveSurfaces(model, active_surface);
401 :
402 : unsigned num_active_model = 0;
403 1433424 : for (unsigned model = 0; model < _num_models; ++model)
404 889120 : if (active_model[model])
405 663504 : num_active_model++;
406 :
407 : ind = 0;
408 544304 : std::vector<unsigned int> active_model_index(_num_models);
409 1433424 : for (unsigned model = 0; model < _num_models; ++model)
410 889120 : if (active_model[model])
411 663504 : active_model_index[model] = ind++;
412 : else
413 225616 : active_model_index[model] =
414 225616 : _num_models + 1; // just a dummy, that will probably cause a crash if something goes wrong
415 :
416 : std::vector<RankTwoTensor> df_dstress;
417 544304 : dyieldFunction_dstress(stress, intnl, active_surface, df_dstress);
418 :
419 : std::vector<Real> df_dintnl;
420 544304 : dyieldFunction_dintnl(stress, intnl, active_surface, df_dintnl);
421 :
422 : std::vector<RankTwoTensor> r;
423 544304 : flowPotential(stress, intnl, active, r);
424 :
425 : std::vector<RankFourTensor> dr_dstress;
426 544304 : dflowPotential_dstress(stress, intnl, active, dr_dstress);
427 :
428 : std::vector<RankTwoTensor> dr_dintnl;
429 544304 : dflowPotential_dintnl(stress, intnl, active, dr_dintnl);
430 :
431 : std::vector<Real> h;
432 544304 : hardPotential(stress, intnl, active, h);
433 :
434 : std::vector<RankTwoTensor> dh_dstress;
435 544304 : dhardPotential_dstress(stress, intnl, active, dh_dstress);
436 :
437 : std::vector<Real> dh_dintnl;
438 544304 : dhardPotential_dintnl(stress, intnl, active, dh_dintnl);
439 :
440 : // d(epp)/dstress = sum_{active alpha} pm[alpha]*dr_dstress
441 544304 : RankFourTensor depp_dstress;
442 : ind = 0;
443 1580928 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
444 1036624 : if (active[surface]) // includes deactivated_due_to_ld
445 1479496 : depp_dstress += pm[surface] * dr_dstress[ind++];
446 544304 : depp_dstress += E_inv;
447 :
448 : // d(epp)/dpm_{active_surface_index} = r_{active_surface_index}
449 : std::vector<RankTwoTensor> depp_dpm;
450 544304 : depp_dpm.resize(num_active_surface);
451 : ind = 0;
452 : active_surface_ind = 0;
453 1580928 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
454 : {
455 1036624 : if (active[surface])
456 : {
457 739748 : if (active_surface[surface]) // do not include the deactived_due_to_ld, since their pm are not
458 : // dofs in the NR
459 676016 : depp_dpm[active_surface_ind++] = r[ind];
460 739748 : ind++;
461 : }
462 : }
463 :
464 : // d(epp)/dintnl_{active_model_index} = sum(pm[asdf]*dr_dintnl[fdsa])
465 : std::vector<RankTwoTensor> depp_dintnl;
466 544304 : depp_dintnl.assign(num_active_model, RankTwoTensor());
467 : ind = 0;
468 1580928 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
469 : {
470 1036624 : if (active[surface])
471 : {
472 739748 : unsigned int model_num = modelNumber(surface);
473 739748 : if (active_model[model_num]) // only include models with surfaces which are still active after
474 : // deactivated_due_to_ld
475 678480 : depp_dintnl[active_model_index[model_num]] += pm[surface] * dr_dintnl[ind];
476 739748 : ind++;
477 : }
478 : }
479 :
480 : // df_dstress has been calculated above
481 : // df_dpm is always zero
482 : // df_dintnl has been calculated above, but only the active_surface+active_model stuff needs to be
483 : // included in Jacobian: see below
484 :
485 : std::vector<RankTwoTensor> dic_dstress;
486 544304 : dic_dstress.assign(num_active_model, RankTwoTensor());
487 : ind = 0;
488 1580928 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
489 : {
490 1036624 : if (active[surface])
491 : {
492 739748 : unsigned int model_num = modelNumber(surface);
493 739748 : if (active_model[model_num]) // only include ic for models with active_surface (ie, if model
494 : // only contains deactivated_due_to_ld don't include it)
495 678480 : dic_dstress[active_model_index[model_num]] += pm[surface] * dh_dstress[ind];
496 739748 : ind++;
497 : }
498 : }
499 :
500 : std::vector<std::vector<Real>> dic_dpm;
501 544304 : dic_dpm.resize(num_active_model);
502 : ind = 0;
503 : active_surface_ind = 0;
504 1207808 : for (unsigned model = 0; model < num_active_model; ++model)
505 663504 : dic_dpm[model].assign(num_active_surface, 0);
506 1580928 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
507 : {
508 1036624 : if (active[surface])
509 : {
510 739748 : if (active_surface[surface]) // only take derivs wrt active-but-not-deactivated_due_to_ld pm
511 : {
512 676016 : unsigned int model_num = modelNumber(surface);
513 : // if (active_model[model_num]) // do not need this check as if the surface has
514 : // active_surface, the model must be deemed active!
515 676016 : dic_dpm[active_model_index[model_num]][active_surface_ind] = h[ind];
516 676016 : active_surface_ind++;
517 : }
518 739748 : ind++;
519 : }
520 : }
521 :
522 : std::vector<std::vector<Real>> dic_dintnl;
523 544304 : dic_dintnl.resize(num_active_model);
524 1207808 : for (unsigned model = 0; model < num_active_model; ++model)
525 : {
526 663504 : dic_dintnl[model].assign(num_active_model, 0);
527 663504 : dic_dintnl[model][model] = 1; // deriv wrt internal parameter
528 : }
529 : ind = 0;
530 1580928 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
531 : {
532 1036624 : if (active[surface])
533 : {
534 739748 : unsigned int model_num = modelNumber(surface);
535 739748 : if (active_model[model_num]) // only the models that contain surfaces that are still active
536 : // after deactivation_due_to_ld
537 678480 : dic_dintnl[active_model_index[model_num]][active_model_index[model_num]] +=
538 678480 : pm[surface] * dh_dintnl[ind];
539 739748 : ind++;
540 : }
541 : }
542 :
543 : unsigned int dim = 3;
544 544304 : unsigned int system_size =
545 544304 : 6 + num_active_surface + num_active_model; // "6" comes from symmeterizing epp
546 544304 : jac.resize(system_size);
547 5149648 : for (unsigned i = 0; i < system_size; ++i)
548 4605344 : jac[i].assign(system_size, 0);
549 :
550 : unsigned int row_num = 0;
551 : unsigned int col_num = 0;
552 2177216 : for (unsigned i = 0; i < dim; ++i)
553 4898736 : for (unsigned j = 0; j <= i; ++j)
554 : {
555 13063296 : for (unsigned k = 0; k < dim; ++k)
556 29392416 : for (unsigned l = 0; l <= k; ++l)
557 19594944 : jac[col_num][row_num++] =
558 19594944 : depp_dstress(i, j, k, l) +
559 19594944 : (k != l ? depp_dstress(i, j, l, k)
560 : : 0); // extra part is needed because i assume dstress(i, j) = dstress(j, i)
561 7321920 : for (unsigned surface = 0; surface < num_active_surface; ++surface)
562 4056096 : jac[col_num][row_num++] = depp_dpm[surface](i, j);
563 7246848 : for (unsigned a = 0; a < num_active_model; ++a)
564 3981024 : jac[col_num][row_num++] = depp_dintnl[a](i, j);
565 : row_num = 0;
566 3265824 : col_num++;
567 : }
568 :
569 : ind = 0;
570 1580928 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
571 1036624 : if (active_surface[surface])
572 : {
573 2704064 : for (unsigned k = 0; k < dim; ++k)
574 6084144 : for (unsigned l = 0; l <= k; ++l)
575 4056096 : jac[col_num][row_num++] =
576 6084144 : df_dstress[ind](k, l) +
577 4056096 : (k != l ? df_dstress[ind](l, k)
578 : : 0); // extra part is needed because i assume dstress(i, j) = dstress(j, i)
579 1684736 : for (unsigned beta = 0; beta < num_active_surface; ++beta)
580 1008720 : jac[col_num][row_num++] = 0; // df_dpm
581 2061648 : for (unsigned model = 0; model < _num_models; ++model)
582 1385632 : if (active_model[model]) // only use df_dintnl for models in active_model
583 : {
584 976576 : if (modelNumber(surface) == model)
585 676016 : jac[col_num][row_num++] = df_dintnl[ind];
586 : else
587 300560 : jac[col_num][row_num++] = 0;
588 : }
589 676016 : ind++;
590 : row_num = 0;
591 676016 : col_num++;
592 : }
593 :
594 1207808 : for (unsigned a = 0; a < num_active_model; ++a)
595 : {
596 2654016 : for (unsigned k = 0; k < dim; ++k)
597 5971536 : for (unsigned l = 0; l <= k; ++l)
598 3981024 : jac[col_num][row_num++] =
599 5971536 : dic_dstress[a](k, l) +
600 3981024 : (k != l ? dic_dstress[a](l, k)
601 : : 0); // extra part is needed because i assume dstress(i, j) = dstress(j, i)
602 1640080 : for (unsigned alpha = 0; alpha < num_active_surface; ++alpha)
603 976576 : jac[col_num][row_num++] = dic_dpm[a][alpha];
604 1621152 : for (unsigned b = 0; b < num_active_model; ++b)
605 957648 : jac[col_num][row_num++] = dic_dintnl[a][b];
606 : row_num = 0;
607 663504 : col_num++;
608 : }
609 :
610 : mooseAssert(col_num == system_size, "Incorrect filling of cols in Jacobian");
611 1088608 : }
612 :
613 : void
614 544304 : MultiPlasticityLinearSystem::nrStep(const RankTwoTensor & stress,
615 : const std::vector<Real> & intnl_old,
616 : const std::vector<Real> & intnl,
617 : const std::vector<Real> & pm,
618 : const RankFourTensor & E_inv,
619 : const RankTwoTensor & delta_dp,
620 : RankTwoTensor & dstress,
621 : std::vector<Real> & dpm,
622 : std::vector<Real> & dintnl,
623 : const std::vector<bool> & active,
624 : std::vector<bool> & deactivated_due_to_ld)
625 : {
626 : // Calculate RHS and Jacobian
627 : std::vector<Real> rhs;
628 544304 : calculateRHS(stress, intnl_old, intnl, pm, delta_dp, rhs, active, true, deactivated_due_to_ld);
629 :
630 : std::vector<std::vector<Real>> jac;
631 544304 : calculateJacobian(stress, intnl, pm, E_inv, active, deactivated_due_to_ld, jac);
632 :
633 : // prepare for LAPACKgesv_ routine provided by PETSc
634 544304 : PetscBLASInt system_size = rhs.size();
635 :
636 544304 : std::vector<double> a(system_size * system_size);
637 : // Fill in the a "matrix" by going down columns
638 : unsigned ind = 0;
639 5149648 : for (int col = 0; col < system_size; ++col)
640 44194048 : for (int row = 0; row < system_size; ++row)
641 39588704 : a[ind++] = jac[row][col];
642 :
643 544304 : PetscBLASInt nrhs = 1;
644 544304 : std::vector<PetscBLASInt> ipiv(system_size);
645 : PetscBLASInt info;
646 544304 : LAPACKgesv_(&system_size, &nrhs, &a[0], &system_size, &ipiv[0], &rhs[0], &system_size, &info);
647 :
648 544304 : if (info != 0)
649 0 : mooseError("In solving the linear system in a Newton-Raphson process, the PETSC LAPACK gsev "
650 : "routine returned with error code ",
651 : info);
652 :
653 : // Extract the results back to dstress, dpm and dintnl
654 544304 : std::vector<bool> active_not_deact(_num_surfaces);
655 1580928 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
656 1397232 : active_not_deact[surface] = (active[surface] && !deactivated_due_to_ld[surface]);
657 :
658 : unsigned int dim = 3;
659 : ind = 0;
660 :
661 2177216 : for (unsigned i = 0; i < dim; ++i)
662 4898736 : for (unsigned j = 0; j <= i; ++j)
663 3265824 : dstress(i, j) = dstress(j, i) = rhs[ind++];
664 544304 : dpm.assign(_num_surfaces, 0);
665 1580928 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
666 1036624 : if (active_not_deact[surface])
667 676016 : dpm[surface] = rhs[ind++];
668 544304 : dintnl.assign(_num_models, 0);
669 1433424 : for (unsigned model = 0; model < _num_models; ++model)
670 889120 : if (anyActiveSurfaces(model, active_not_deact))
671 663504 : dintnl[model] = rhs[ind++];
672 :
673 : mooseAssert(static_cast<int>(ind) == system_size,
674 : "Incorrect extracting of changes from NR solution in nrStep");
675 1088608 : }
|
