Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://www.mooseframework.org
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #include "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 2087 : MultiPlasticityLinearSystem::validParams()
21 : {
22 2087 : InputParameters params = MultiPlasticityRawComponentAssembler::validParams();
23 6261 : params.addRangeCheckedParam<Real>("linear_dependent",
24 4174 : 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 2087 : return params;
30 0 : }
31 :
32 1561 : MultiPlasticityLinearSystem::MultiPlasticityLinearSystem(const MooseObject * moose_object)
33 : : MultiPlasticityRawComponentAssembler(moose_object),
34 3122 : _svd_tol(_params.get<Real>("linear_dependent")),
35 1561 : _min_f_tol(-1.0)
36 : {
37 3971 : for (unsigned model = 0; model < _num_models; ++model)
38 2410 : if (_min_f_tol == -1.0 || _min_f_tol > _f[model]->_f_tol)
39 1525 : _min_f_tol = _f[model]->_f_tol;
40 :
41 : MooseRandom::seed(0);
42 1561 : }
43 :
44 : int
45 102932 : MultiPlasticityLinearSystem::singularValuesOfR(const std::vector<RankTwoTensor> & r,
46 : std::vector<Real> & s)
47 : {
48 102932 : PetscBLASInt bm = r.size();
49 102932 : PetscBLASInt bn = 6;
50 :
51 102932 : 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 102932 : std::vector<double> a(bm * 6);
63 : // Fill in the a "matrix" by going down columns
64 : unsigned ind = 0;
65 411728 : for (int col = 0; col < 3; ++col)
66 1166874 : for (int row = 0; row < bm; ++row)
67 858078 : a[ind++] = r[row](0, col);
68 308796 : for (int col = 3; col < 5; ++col)
69 777916 : for (int row = 0; row < bm; ++row)
70 572052 : a[ind++] = r[row](1, col - 2);
71 388958 : for (int row = 0; row < bm; ++row)
72 286026 : 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 102932 : PetscBLASInt sizeu = 1;
77 102932 : std::vector<double> u(sizeu);
78 102932 : PetscBLASInt sizevt = 1;
79 102932 : std::vector<double> vt(sizevt);
80 :
81 102932 : PetscBLASInt sizework =
82 102932 : 16 * (bm + 6); // this is above the lowerbound specified in the LAPACK doco
83 102932 : std::vector<double> work(sizework);
84 :
85 : PetscBLASInt info;
86 :
87 102932 : 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 205864 : return info;
103 : }
104 :
105 : void
106 281794 : 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 281794 : deactivated_due_to_ld.resize(_num_surfaces, false);
114 :
115 281794 : unsigned num_active = r.size();
116 :
117 281794 : if (num_active <= 1)
118 271064 : return;
119 :
120 : std::vector<double> s;
121 58306 : int info = singularValuesOfR(r, s);
122 58306 : 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 58306 : unsigned i = s.size();
132 87486 : while (i-- > 0)
133 87486 : if (s[i] < _svd_tol * s[0])
134 29180 : num_lin_dep++;
135 : else
136 : break;
137 :
138 58306 : 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 10730 : dyieldFunction_dstress(stress, intnl, active, df_dstress);
149 :
150 : typedef std::pair<Real, unsigned> pair_for_sorting;
151 10730 : std::vector<pair_for_sorting> dist(num_active);
152 66120 : for (unsigned i = 0; i < num_active; ++i)
153 : {
154 55390 : dist[i].first = f[i] / df_dstress[i].L2norm();
155 55390 : dist[i].second = i;
156 : }
157 10730 : 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 55390 : for (unsigned i = 0; i < num_active - 1; ++i)
162 44660 : if (std::abs(dist[i].first - dist[i + 1].first) < _min_f_tol)
163 : {
164 : equals_detected = true;
165 5158 : dist[i].first += _min_f_tol * (MooseRandom::rand() - 0.5);
166 : }
167 10730 : if (equals_detected)
168 3728 : 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 10730 : current_yf = dist[num_active - 1].second;
181 : // the one with largest dist
182 10730 : std::vector<RankTwoTensor> r_tmp = {r[current_yf]};
183 :
184 : unsigned num_kept_active = 1;
185 46768 : for (unsigned yf_to_try = 2; yf_to_try <= num_active; ++yf_to_try)
186 : {
187 44626 : current_yf = dist[num_active - yf_to_try].second;
188 44626 : if (num_active == 2) // shortcut to we don't have to singularValuesOfR
189 : scheduled_for_deactivation[current_yf] = true;
190 44626 : 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 44626 : r_tmp.push_back(r[current_yf]);
196 44626 : info = singularValuesOfR(r_tmp, s);
197 44626 : 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 44626 : if (s[s.size() - 1] < _svd_tol * s[0])
202 : {
203 : scheduled_for_deactivation[current_yf] = true;
204 : r_tmp.pop_back();
205 32698 : num_lin_dep--;
206 : }
207 : else
208 11928 : num_kept_active++;
209 44626 : 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 84858 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
218 74128 : if (active[surface])
219 : {
220 55390 : if (scheduled_for_deactivation[old_active_number])
221 : deactivated_due_to_ld[surface] = true;
222 55390 : old_active_number++;
223 : }
224 10730 : }
225 :
226 : void
227 697414 : 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 697414 : yieldFunction(stress, intnl, active, f);
253 :
254 : // flow directions and "epp"
255 697414 : flowPotential(stress, intnl, active, r);
256 697414 : epp = RankTwoTensor();
257 : unsigned ind = 0;
258 1993630 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
259 1296216 : if (active[surface])
260 959152 : epp += pm[surface] * r[ind++]; // note, even the deactivated_due_to_ld must get added in
261 697414 : epp -= delta_dp;
262 :
263 : // internal constraints
264 : std::vector<Real> h;
265 697414 : hardPotential(stress, intnl, active, h);
266 697414 : ic.resize(0);
267 : ind = 0;
268 : std::vector<unsigned int> active_surfaces;
269 : std::vector<unsigned int>::iterator active_surface;
270 1845678 : for (unsigned model = 0; model < _num_models; ++model)
271 : {
272 1148264 : activeSurfaces(model, active, active_surfaces);
273 1148264 : if (active_surfaces.size() > 0)
274 : {
275 : // some surfaces are active in this model, so must form an internal constraint
276 944112 : ic.push_back(intnl[model] - intnl_old[model]);
277 1903264 : for (active_surface = active_surfaces.begin(); active_surface != active_surfaces.end();
278 : ++active_surface)
279 959152 : 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 697414 : }
285 :
286 : void
287 281794 : 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 281794 : RankTwoTensor epp; // the plastic-strain constraint ("direction constraint")
310 : std::vector<Real> ic; // the "internal constraints"
311 :
312 : std::vector<RankTwoTensor> r;
313 281794 : calculateConstraints(stress, intnl_old, intnl, pm, delta_dp, f, r, epp, ic, active);
314 :
315 281794 : if (eliminate_ld)
316 281794 : eliminateLinearDependence(stress, intnl, f, r, active, deactivated_due_to_ld);
317 : else
318 0 : deactivated_due_to_ld.assign(_num_surfaces, false);
319 :
320 281794 : std::vector<bool> active_not_deact(_num_surfaces);
321 820586 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
322 721044 : active_not_deact[surface] = (active[surface] && !deactivated_due_to_ld[surface]);
323 :
324 : unsigned num_active_f = 0;
325 820586 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
326 538792 : if (active_not_deact[surface])
327 356540 : num_active_f++;
328 :
329 : unsigned num_active_ic = 0;
330 746610 : for (unsigned model = 0; model < _num_models; ++model)
331 464816 : if (anyActiveSurfaces(model, active_not_deact))
332 350252 : num_active_ic++;
333 :
334 : unsigned int dim = 3;
335 281794 : 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 281794 : rhs.resize(system_size);
340 :
341 : unsigned ind = 0;
342 1127176 : for (unsigned i = 0; i < dim; ++i)
343 2536146 : for (unsigned j = 0; j <= i; ++j)
344 1690764 : rhs[ind++] = -epp(i, j);
345 : unsigned active_surface = 0;
346 820586 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
347 538792 : if (active[surface])
348 : {
349 389238 : if (!deactivated_due_to_ld[surface])
350 356540 : rhs[ind++] = -f[active_surface];
351 389238 : active_surface++;
352 : }
353 : unsigned active_model = 0;
354 746610 : for (unsigned model = 0; model < _num_models; ++model)
355 464816 : if (anyActiveSurfaces(model, active))
356 : {
357 381718 : if (anyActiveSurfaces(model, active_not_deact))
358 350252 : rhs[ind++] = -ic[active_model];
359 381718 : active_model++;
360 : }
361 :
362 : mooseAssert(ind == system_size, "Incorrect filling of the rhs in calculateRHS");
363 281794 : }
364 :
365 : void
366 281794 : 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 281794 : std::vector<bool> active_surface(_num_surfaces); // active and not deactivated_due_to_ld
390 820586 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
391 721044 : active_surface[surface] = (active[surface] && !deactivated_due_to_ld[surface]);
392 : unsigned num_active_surface = 0;
393 820586 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
394 538792 : if (active_surface[surface])
395 356540 : num_active_surface++;
396 :
397 : std::vector<bool> active_model(
398 281794 : _num_models); // whether a model has surfaces that are active and not deactivated_due_to_ld
399 746610 : for (unsigned model = 0; model < _num_models; ++model)
400 464816 : active_model[model] = anyActiveSurfaces(model, active_surface);
401 :
402 : unsigned num_active_model = 0;
403 746610 : for (unsigned model = 0; model < _num_models; ++model)
404 464816 : if (active_model[model])
405 350252 : num_active_model++;
406 :
407 : ind = 0;
408 281794 : std::vector<unsigned int> active_model_index(_num_models);
409 746610 : for (unsigned model = 0; model < _num_models; ++model)
410 464816 : if (active_model[model])
411 350252 : active_model_index[model] = ind++;
412 : else
413 114564 : active_model_index[model] =
414 114564 : _num_models + 1; // just a dummy, that will probably cause a crash if something goes wrong
415 :
416 : std::vector<RankTwoTensor> df_dstress;
417 281794 : dyieldFunction_dstress(stress, intnl, active_surface, df_dstress);
418 :
419 : std::vector<Real> df_dintnl;
420 281794 : dyieldFunction_dintnl(stress, intnl, active_surface, df_dintnl);
421 :
422 : std::vector<RankTwoTensor> r;
423 281794 : flowPotential(stress, intnl, active, r);
424 :
425 : std::vector<RankFourTensor> dr_dstress;
426 281794 : dflowPotential_dstress(stress, intnl, active, dr_dstress);
427 :
428 : std::vector<RankTwoTensor> dr_dintnl;
429 281794 : dflowPotential_dintnl(stress, intnl, active, dr_dintnl);
430 :
431 : std::vector<Real> h;
432 281794 : hardPotential(stress, intnl, active, h);
433 :
434 : std::vector<RankTwoTensor> dh_dstress;
435 281794 : dhardPotential_dstress(stress, intnl, active, dh_dstress);
436 :
437 : std::vector<Real> dh_dintnl;
438 281794 : dhardPotential_dintnl(stress, intnl, active, dh_dintnl);
439 :
440 : // d(epp)/dstress = sum_{active alpha} pm[alpha]*dr_dstress
441 281794 : RankFourTensor depp_dstress;
442 : ind = 0;
443 820586 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
444 538792 : if (active[surface]) // includes deactivated_due_to_ld
445 778476 : depp_dstress += pm[surface] * dr_dstress[ind++];
446 281794 : depp_dstress += E_inv;
447 :
448 : // d(epp)/dpm_{active_surface_index} = r_{active_surface_index}
449 : std::vector<RankTwoTensor> depp_dpm;
450 281794 : depp_dpm.resize(num_active_surface);
451 : ind = 0;
452 : active_surface_ind = 0;
453 820586 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
454 : {
455 538792 : if (active[surface])
456 : {
457 389238 : if (active_surface[surface]) // do not include the deactived_due_to_ld, since their pm are not
458 : // dofs in the NR
459 356540 : depp_dpm[active_surface_ind++] = r[ind];
460 389238 : ind++;
461 : }
462 : }
463 :
464 : // d(epp)/dintnl_{active_model_index} = sum(pm[asdf]*dr_dintnl[fdsa])
465 : std::vector<RankTwoTensor> depp_dintnl;
466 281794 : depp_dintnl.assign(num_active_model, RankTwoTensor());
467 : ind = 0;
468 820586 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
469 : {
470 538792 : if (active[surface])
471 : {
472 389238 : unsigned int model_num = modelNumber(surface);
473 389238 : if (active_model[model_num]) // only include models with surfaces which are still active after
474 : // deactivated_due_to_ld
475 357772 : depp_dintnl[active_model_index[model_num]] += pm[surface] * dr_dintnl[ind];
476 389238 : 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 281794 : dic_dstress.assign(num_active_model, RankTwoTensor());
487 : ind = 0;
488 820586 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
489 : {
490 538792 : if (active[surface])
491 : {
492 389238 : unsigned int model_num = modelNumber(surface);
493 389238 : 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 357772 : dic_dstress[active_model_index[model_num]] += pm[surface] * dh_dstress[ind];
496 389238 : ind++;
497 : }
498 : }
499 :
500 : std::vector<std::vector<Real>> dic_dpm;
501 281794 : dic_dpm.resize(num_active_model);
502 : ind = 0;
503 : active_surface_ind = 0;
504 632046 : for (unsigned model = 0; model < num_active_model; ++model)
505 350252 : dic_dpm[model].assign(num_active_surface, 0);
506 820586 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
507 : {
508 538792 : if (active[surface])
509 : {
510 389238 : if (active_surface[surface]) // only take derivs wrt active-but-not-deactivated_due_to_ld pm
511 : {
512 356540 : 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 356540 : dic_dpm[active_model_index[model_num]][active_surface_ind] = h[ind];
516 356540 : active_surface_ind++;
517 : }
518 389238 : ind++;
519 : }
520 : }
521 :
522 : std::vector<std::vector<Real>> dic_dintnl;
523 281794 : dic_dintnl.resize(num_active_model);
524 632046 : for (unsigned model = 0; model < num_active_model; ++model)
525 : {
526 350252 : dic_dintnl[model].assign(num_active_model, 0);
527 350252 : dic_dintnl[model][model] = 1; // deriv wrt internal parameter
528 : }
529 : ind = 0;
530 820586 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
531 : {
532 538792 : if (active[surface])
533 : {
534 389238 : unsigned int model_num = modelNumber(surface);
535 389238 : if (active_model[model_num]) // only the models that contain surfaces that are still active
536 : // after deactivation_due_to_ld
537 357772 : dic_dintnl[active_model_index[model_num]][active_model_index[model_num]] +=
538 357772 : pm[surface] * dh_dintnl[ind];
539 389238 : ind++;
540 : }
541 : }
542 :
543 : unsigned int dim = 3;
544 281794 : unsigned int system_size =
545 281794 : 6 + num_active_surface + num_active_model; // "6" comes from symmeterizing epp
546 281794 : jac.resize(system_size);
547 2679350 : for (unsigned i = 0; i < system_size; ++i)
548 2397556 : jac[i].assign(system_size, 0);
549 :
550 : unsigned int row_num = 0;
551 : unsigned int col_num = 0;
552 1127176 : for (unsigned i = 0; i < dim; ++i)
553 2536146 : for (unsigned j = 0; j <= i; ++j)
554 : {
555 6763056 : for (unsigned k = 0; k < dim; ++k)
556 15216876 : for (unsigned l = 0; l <= k; ++l)
557 10144584 : jac[col_num][row_num++] =
558 10144584 : depp_dstress(i, j, k, l) +
559 10144584 : (k != l ? depp_dstress(i, j, l, k)
560 : : 0); // extra part is needed because i assume dstress(i, j) = dstress(j, i)
561 3830004 : for (unsigned surface = 0; surface < num_active_surface; ++surface)
562 2139240 : jac[col_num][row_num++] = depp_dpm[surface](i, j);
563 3792276 : for (unsigned a = 0; a < num_active_model; ++a)
564 2101512 : jac[col_num][row_num++] = depp_dintnl[a](i, j);
565 : row_num = 0;
566 1690764 : col_num++;
567 : }
568 :
569 : ind = 0;
570 820586 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
571 538792 : if (active_surface[surface])
572 : {
573 1426160 : for (unsigned k = 0; k < dim; ++k)
574 3208860 : for (unsigned l = 0; l <= k; ++l)
575 2139240 : jac[col_num][row_num++] =
576 3208860 : df_dstress[ind](k, l) +
577 2139240 : (k != l ? df_dstress[ind](l, k)
578 : : 0); // extra part is needed because i assume dstress(i, j) = dstress(j, i)
579 897212 : for (unsigned beta = 0; beta < num_active_surface; ++beta)
580 540672 : jac[col_num][row_num++] = 0; // df_dpm
581 1088572 : for (unsigned model = 0; model < _num_models; ++model)
582 732032 : if (active_model[model]) // only use df_dintnl for models in active_model
583 : {
584 524536 : if (modelNumber(surface) == model)
585 356540 : jac[col_num][row_num++] = df_dintnl[ind];
586 : else
587 167996 : jac[col_num][row_num++] = 0;
588 : }
589 356540 : ind++;
590 : row_num = 0;
591 356540 : col_num++;
592 : }
593 :
594 632046 : for (unsigned a = 0; a < num_active_model; ++a)
595 : {
596 1401008 : for (unsigned k = 0; k < dim; ++k)
597 3152268 : for (unsigned l = 0; l <= k; ++l)
598 2101512 : jac[col_num][row_num++] =
599 3152268 : dic_dstress[a](k, l) +
600 2101512 : (k != l ? dic_dstress[a](l, k)
601 : : 0); // extra part is needed because i assume dstress(i, j) = dstress(j, i)
602 874788 : for (unsigned alpha = 0; alpha < num_active_surface; ++alpha)
603 524536 : jac[col_num][row_num++] = dic_dpm[a][alpha];
604 865292 : for (unsigned b = 0; b < num_active_model; ++b)
605 515040 : jac[col_num][row_num++] = dic_dintnl[a][b];
606 : row_num = 0;
607 350252 : col_num++;
608 : }
609 :
610 : mooseAssert(col_num == system_size, "Incorrect filling of cols in Jacobian");
611 563588 : }
612 :
613 : void
614 281794 : 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 281794 : calculateRHS(stress, intnl_old, intnl, pm, delta_dp, rhs, active, true, deactivated_due_to_ld);
629 :
630 : std::vector<std::vector<Real>> jac;
631 281794 : calculateJacobian(stress, intnl, pm, E_inv, active, deactivated_due_to_ld, jac);
632 :
633 : // prepare for LAPACKgesv_ routine provided by PETSc
634 281794 : PetscBLASInt system_size = rhs.size();
635 :
636 281794 : std::vector<double> a(system_size * system_size);
637 : // Fill in the a "matrix" by going down columns
638 : unsigned ind = 0;
639 2679350 : for (int col = 0; col < system_size; ++col)
640 23128428 : for (int row = 0; row < system_size; ++row)
641 20730872 : a[ind++] = jac[row][col];
642 :
643 281794 : PetscBLASInt nrhs = 1;
644 281794 : std::vector<PetscBLASInt> ipiv(system_size);
645 : PetscBLASInt info;
646 281794 : LAPACKgesv_(&system_size, &nrhs, &a[0], &system_size, &ipiv[0], &rhs[0], &system_size, &info);
647 :
648 281794 : 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 281794 : std::vector<bool> active_not_deact(_num_surfaces);
655 820586 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
656 721044 : active_not_deact[surface] = (active[surface] && !deactivated_due_to_ld[surface]);
657 :
658 : unsigned int dim = 3;
659 : ind = 0;
660 :
661 1127176 : for (unsigned i = 0; i < dim; ++i)
662 2536146 : for (unsigned j = 0; j <= i; ++j)
663 1690764 : dstress(i, j) = dstress(j, i) = rhs[ind++];
664 281794 : dpm.assign(_num_surfaces, 0);
665 820586 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
666 538792 : if (active_not_deact[surface])
667 356540 : dpm[surface] = rhs[ind++];
668 281794 : dintnl.assign(_num_models, 0);
669 746610 : for (unsigned model = 0; model < _num_models; ++model)
670 464816 : if (anyActiveSurfaces(model, active_not_deact))
671 350252 : dintnl[model] = rhs[ind++];
672 :
673 : mooseAssert(static_cast<int>(ind) == system_size,
674 : "Incorrect extracting of changes from NR solution in nrStep");
675 563588 : }
|
