17 #include <petscblaslapack.h> 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");
34 _svd_tol(_params.
get<
Real>(
"linear_dependent")),
46 std::vector<Real> & s)
48 PetscBLASInt bm = r.size();
51 s.resize(std::min(bm, bn));
62 std::vector<double>
a(bm * 6);
65 for (
int col = 0; col < 3; ++col)
66 for (
int row = 0; row < bm; ++row)
67 a[ind++] = r[row](0, col);
68 for (
int col = 3; col < 5; ++col)
69 for (
int row = 0; row < bm; ++row)
70 a[ind++] = r[row](1, col - 2);
71 for (
int row = 0; row < bm; ++row)
72 a[ind++] = r[row](2, 2);
76 PetscBLASInt sizeu = 1;
77 std::vector<double> u(sizeu);
78 PetscBLASInt sizevt = 1;
79 std::vector<double> vt(sizevt);
81 PetscBLASInt sizework =
83 std::vector<double> work(sizework);
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)
115 unsigned num_active = r.size();
120 std::vector<double> s;
123 mooseError(
"In finding the SVD in the return-map algorithm, the PETSC LAPACK gesvd routine " 124 "returned with error code ",
129 unsigned int num_lin_dep = 0;
131 unsigned i = s.size();
138 if (num_lin_dep == 0 && num_active <= 6)
147 std::vector<RankTwoTensor> df_dstress;
150 typedef std::pair<Real, unsigned> pair_for_sorting;
151 std::vector<pair_for_sorting> dist(num_active);
152 for (
unsigned i = 0; i < num_active; ++i)
154 dist[i].first =
f[i] / df_dstress[i].L2norm();
157 std::sort(dist.begin(), dist.end());
160 bool equals_detected =
false;
161 for (
unsigned i = 0; i < num_active - 1; ++i)
162 if (std::abs(dist[i].first - dist[i + 1].first) <
_min_f_tol)
164 equals_detected =
true;
168 std::sort(dist.begin(), dist.end());
170 std::vector<bool> scheduled_for_deactivation;
171 scheduled_for_deactivation.assign(num_active,
false);
180 current_yf = dist[num_active - 1].second;
182 std::vector<RankTwoTensor> r_tmp = {r[current_yf]};
184 unsigned num_kept_active = 1;
185 for (
unsigned yf_to_try = 2; yf_to_try <= num_active; ++yf_to_try)
187 current_yf = dist[num_active - yf_to_try].second;
189 scheduled_for_deactivation[current_yf] =
true;
190 else if (num_kept_active >= 6)
192 scheduled_for_deactivation[current_yf] =
true;
195 r_tmp.push_back(r[current_yf]);
198 mooseError(
"In finding the SVD in the return-map algorithm, the PETSC LAPACK gesvd routine " 199 "returned with error code ",
201 if (s[s.size() - 1] <
_svd_tol * s[0])
203 scheduled_for_deactivation[current_yf] =
true;
209 if (num_lin_dep == 0 && num_active <= 6)
216 unsigned int old_active_number = 0;
217 for (
unsigned surface = 0; surface <
_num_surfaces; ++surface)
220 if (scheduled_for_deactivation[old_active_number])
221 deactivated_due_to_ld[surface] =
true;
228 const std::vector<Real> & intnl_old,
229 const std::vector<Real> & intnl,
230 const std::vector<Real> & pm,
232 std::vector<Real> &
f,
233 std::vector<RankTwoTensor> & r,
235 std::vector<Real> & ic,
236 const std::vector<bool> & active)
241 "Size of intnl_old is " << intnl_old.size()
242 <<
" which is incorrect in calculateConstraints");
244 "Size of intnl is " << intnl.size() <<
" which is incorrect in calculateConstraints");
246 "Size of pm is " << pm.size() <<
" which is incorrect in calculateConstraints");
248 "Size of active is " << active.size()
249 <<
" which is incorrect in calculateConstraints");
258 for (
unsigned surface = 0; surface <
_num_surfaces; ++surface)
260 epp += pm[surface] * r[ind++];
268 std::vector<unsigned int> active_surfaces;
269 std::vector<unsigned int>::iterator active_surface;
273 if (active_surfaces.size() > 0)
276 ic.push_back(intnl[
model] - intnl_old[
model]);
277 for (active_surface = active_surfaces.begin(); active_surface != active_surfaces.end();
279 ic[ic.size() - 1] += pm[*active_surface] * h[ind++];
288 const std::vector<Real> & intnl_old,
289 const std::vector<Real> & intnl,
290 const std::vector<Real> & pm,
292 std::vector<Real> & rhs,
293 const std::vector<bool> & active,
295 std::vector<bool> & deactivated_due_to_ld)
300 "Size of intnl_old is " << intnl_old.size() <<
" which is incorrect in calculateRHS");
302 "Size of intnl is " << intnl.size() <<
" which is incorrect in calculateRHS");
304 "Size of pm is " << pm.size() <<
" which is incorrect in calculateRHS");
306 "Size of active is " << active.size() <<
" which is incorrect in calculateRHS");
310 std::vector<Real> ic;
312 std::vector<RankTwoTensor> r;
321 for (
unsigned surface = 0; surface <
_num_surfaces; ++surface)
322 active_not_deact[surface] = (active[surface] && !deactivated_due_to_ld[surface]);
324 unsigned num_active_f = 0;
325 for (
unsigned surface = 0; surface <
_num_surfaces; ++surface)
326 if (active_not_deact[surface])
329 unsigned num_active_ic = 0;
334 unsigned int dim = 3;
335 unsigned int system_size = 6 + num_active_f + num_active_ic;
339 rhs.resize(system_size);
342 for (
unsigned i = 0; i <
dim; ++i)
343 for (
unsigned j = 0;
j <= i; ++
j)
344 rhs[ind++] = -epp(i,
j);
345 unsigned active_surface = 0;
346 for (
unsigned surface = 0; surface <
_num_surfaces; ++surface)
349 if (!deactivated_due_to_ld[surface])
350 rhs[ind++] = -
f[active_surface];
353 unsigned active_model = 0;
358 rhs[ind++] = -ic[active_model];
362 mooseAssert(ind == system_size,
"Incorrect filling of the rhs in calculateRHS");
367 const std::vector<Real> & intnl,
368 const std::vector<Real> & pm,
370 const std::vector<bool> & active,
371 const std::vector<bool> & deactivated_due_to_ld,
372 std::vector<std::vector<Real>> & jac)
377 "Size of intnl is " << intnl.size() <<
" which is incorrect in calculateJacobian");
379 "Size of pm is " << pm.size() <<
" which is incorrect in calculateJacobian");
381 "Size of active is " << active.size() <<
" which is incorrect in calculateJacobian");
383 "Size of deactivated_due_to_ld is " << deactivated_due_to_ld.size()
384 <<
" which is incorrect in calculateJacobian");
387 unsigned active_surface_ind = 0;
390 for (
unsigned surface = 0; surface <
_num_surfaces; ++surface)
391 active_surface[surface] = (active[surface] && !deactivated_due_to_ld[surface]);
392 unsigned num_active_surface = 0;
393 for (
unsigned surface = 0; surface <
_num_surfaces; ++surface)
394 if (active_surface[surface])
395 num_active_surface++;
397 std::vector<bool> active_model(
402 unsigned num_active_model = 0;
404 if (active_model[
model])
408 std::vector<unsigned int> active_model_index(
_num_models);
410 if (active_model[
model])
411 active_model_index[
model] = ind++;
413 active_model_index[
model] =
416 std::vector<RankTwoTensor> df_dstress;
419 std::vector<Real> df_dintnl;
422 std::vector<RankTwoTensor> r;
425 std::vector<RankFourTensor> dr_dstress;
428 std::vector<RankTwoTensor> dr_dintnl;
434 std::vector<RankTwoTensor> dh_dstress;
437 std::vector<Real> dh_dintnl;
443 for (
unsigned surface = 0; surface <
_num_surfaces; ++surface)
445 depp_dstress += pm[surface] * dr_dstress[ind++];
446 depp_dstress += E_inv;
449 std::vector<RankTwoTensor> depp_dpm;
450 depp_dpm.resize(num_active_surface);
452 active_surface_ind = 0;
453 for (
unsigned surface = 0; surface <
_num_surfaces; ++surface)
457 if (active_surface[surface])
459 depp_dpm[active_surface_ind++] = r[ind];
465 std::vector<RankTwoTensor> depp_dintnl;
468 for (
unsigned surface = 0; surface <
_num_surfaces; ++surface)
473 if (active_model[model_num])
475 depp_dintnl[active_model_index[model_num]] += pm[surface] * dr_dintnl[ind];
485 std::vector<RankTwoTensor> dic_dstress;
488 for (
unsigned surface = 0; surface <
_num_surfaces; ++surface)
493 if (active_model[model_num])
495 dic_dstress[active_model_index[model_num]] += pm[surface] * dh_dstress[ind];
500 std::vector<std::vector<Real>> dic_dpm;
501 dic_dpm.resize(num_active_model);
503 active_surface_ind = 0;
505 dic_dpm[
model].assign(num_active_surface, 0);
506 for (
unsigned surface = 0; surface <
_num_surfaces; ++surface)
510 if (active_surface[surface])
515 dic_dpm[active_model_index[model_num]][active_surface_ind] = h[ind];
516 active_surface_ind++;
522 std::vector<std::vector<Real>> dic_dintnl;
523 dic_dintnl.resize(num_active_model);
526 dic_dintnl[
model].assign(num_active_model, 0);
530 for (
unsigned surface = 0; surface <
_num_surfaces; ++surface)
535 if (active_model[model_num])
537 dic_dintnl[active_model_index[model_num]][active_model_index[model_num]] +=
538 pm[surface] * dh_dintnl[ind];
543 unsigned int dim = 3;
544 unsigned int system_size =
545 6 + num_active_surface + num_active_model;
546 jac.resize(system_size);
547 for (
unsigned i = 0; i < system_size; ++i)
548 jac[i].assign(system_size, 0);
550 unsigned int row_num = 0;
551 unsigned int col_num = 0;
552 for (
unsigned i = 0; i <
dim; ++i)
553 for (
unsigned j = 0;
j <= i; ++
j)
555 for (
unsigned k = 0;
k <
dim; ++
k)
556 for (
unsigned l = 0; l <=
k; ++l)
557 jac[col_num][row_num++] =
558 depp_dstress(i,
j,
k, l) +
559 (
k != l ? depp_dstress(i,
j, l,
k)
561 for (
unsigned surface = 0; surface < num_active_surface; ++surface)
562 jac[col_num][row_num++] = depp_dpm[surface](i,
j);
563 for (
unsigned a = 0;
a < num_active_model; ++
a)
564 jac[col_num][row_num++] = depp_dintnl[
a](i,
j);
570 for (
unsigned surface = 0; surface <
_num_surfaces; ++surface)
571 if (active_surface[surface])
573 for (
unsigned k = 0;
k <
dim; ++
k)
574 for (
unsigned l = 0; l <=
k; ++l)
575 jac[col_num][row_num++] =
576 df_dstress[ind](
k, l) +
577 (
k != l ? df_dstress[ind](l,
k)
579 for (
unsigned beta = 0; beta < num_active_surface; ++beta)
580 jac[col_num][row_num++] = 0;
582 if (active_model[
model])
585 jac[col_num][row_num++] = df_dintnl[ind];
587 jac[col_num][row_num++] = 0;
594 for (
unsigned a = 0;
a < num_active_model; ++
a)
596 for (
unsigned k = 0;
k <
dim; ++
k)
597 for (
unsigned l = 0; l <=
k; ++l)
598 jac[col_num][row_num++] =
599 dic_dstress[
a](
k, l) +
600 (
k != l ? dic_dstress[
a](l,
k)
603 jac[col_num][row_num++] = dic_dpm[
a][
alpha];
604 for (
unsigned b = 0;
b < num_active_model; ++
b)
605 jac[col_num][row_num++] = dic_dintnl[
a][
b];
610 mooseAssert(col_num == system_size,
"Incorrect filling of cols in Jacobian");
615 const std::vector<Real> & intnl_old,
616 const std::vector<Real> & intnl,
617 const std::vector<Real> & pm,
621 std::vector<Real> & dpm,
622 std::vector<Real> & dintnl,
623 const std::vector<bool> & active,
624 std::vector<bool> & deactivated_due_to_ld)
627 std::vector<Real> rhs;
628 calculateRHS(stress, intnl_old, intnl, pm, delta_dp, rhs, active,
true, deactivated_due_to_ld);
630 std::vector<std::vector<Real>> jac;
634 PetscBLASInt system_size = rhs.size();
636 std::vector<double>
a(system_size * system_size);
639 for (
int col = 0; col < system_size; ++col)
640 for (
int row = 0; row < system_size; ++row)
641 a[ind++] = jac[row][col];
643 PetscBLASInt nrhs = 1;
644 std::vector<PetscBLASInt> ipiv(system_size);
646 LAPACKgesv_(&system_size, &nrhs, &
a[0], &system_size, &ipiv[0], &rhs[0], &system_size, &
info);
649 mooseError(
"In solving the linear system in a Newton-Raphson process, the PETSC LAPACK gsev " 650 "routine returned with error code ",
655 for (
unsigned surface = 0; surface <
_num_surfaces; ++surface)
656 active_not_deact[surface] = (active[surface] && !deactivated_due_to_ld[surface]);
658 unsigned int dim = 3;
661 for (
unsigned i = 0; i <
dim; ++i)
662 for (
unsigned j = 0;
j <= i; ++
j)
663 dstress(i,
j) = dstress(
j, i) = rhs[ind++];
665 for (
unsigned surface = 0; surface <
_num_surfaces; ++surface)
666 if (active_not_deact[surface])
667 dpm[surface] = rhs[ind++];
671 dintnl[
model] = rhs[ind++];
673 mooseAssert(static_cast<int>(ind) == system_size,
674 "Incorrect extracting of changes from NR solution in nrStep");
virtual void calculateRHS(const RankTwoTensor &stress, const std::vector< Real > &intnl_old, const std::vector< Real > &intnl, const std::vector< Real > &pm, const RankTwoTensor &delta_dp, std::vector< Real > &rhs, const std::vector< bool > &active, bool eliminate_ld, std::vector< bool > &deactivated_due_to_ld)
Calculate the RHS which is rhs = -(epp(0,0), epp(1,0), epp(1,1), epp(2,0), epp(2,1), epp(2,2), f[0], f[1], ..., f[num_f], ic[0], ic[1], ..., ic[num_ic])
virtual int singularValuesOfR(const std::vector< RankTwoTensor > &r, std::vector< Real > &s)
Performs a singular-value decomposition of r and returns the singular values.
virtual void dyieldFunction_dstress(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< RankTwoTensor > &df_dstress)
The derivative of the active yield function(s) with respect to stress.
Real _svd_tol
Tolerance on the minimum ratio of singular values before flow-directions are deemed linearly dependen...
unsigned int modelNumber(unsigned int surface)
returns the model number, given the surface number
void mooseError(Args &&... args)
virtual void dhardPotential_dintnl(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< Real > &dh_dintnl)
The derivative of the active hardening potentials with respect to the active internal parameters...
virtual void calculateJacobian(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< Real > &pm, const RankFourTensor &E_inv, const std::vector< bool > &active, const std::vector< bool > &deactivated_due_to_ld, std::vector< std::vector< Real >> &jac)
d(rhs)/d(dof)
void seed(std::size_t i, unsigned int seed)
virtual void eliminateLinearDependence(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< Real > &f, const std::vector< RankTwoTensor > &r, const std::vector< bool > &active, std::vector< bool > &deactivated_due_to_ld)
Performs a number of singular-value decompositions to check for linear-dependence of the active direc...
virtual void yieldFunction(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< Real > &f)
The active yield function(s)
static InputParameters validParams()
virtual void dhardPotential_dstress(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< RankTwoTensor > &dh_dstress)
The derivative of the active hardening potentials with respect to stress By assumption in the Userobj...
unsigned int _num_surfaces
Number of surfaces within the plastic models.
MultiPlasticityRawComponentAssembler holds and computes yield functions, flow directions, etc, for use in FiniteStrainMultiPlasticity.
MultiPlasticityLinearSystem(const MooseObject *moose_object)
Real f(Real x)
Test function for Brents method.
void assign(const TypeTensor< T2 > &)
virtual void dflowPotential_dintnl(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< RankTwoTensor > &dr_dintnl)
The derivative of the active flow potentials with respect to the active internal parameters The UserO...
bool anyActiveSurfaces(int model, const std::vector< bool > &active)
returns true if any internal surfaces of the given model are active according to 'active' ...
virtual void flowPotential(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< RankTwoTensor > &r)
The active flow potential(s) - one for each yield function.
const PertinentGeochemicalSystem model(database, {"H2O", "H+", "HCO3-", "O2(aq)", "Ca++", ">(s)FeOH", "radius_neg1", "radius_neg1.5"}, {"Calcite"}, {}, {"Calcite_asdf"}, {"CH4(aq)"}, {">(s)FeOCa+"}, "O2(aq)", "e-")
static InputParameters validParams()
unsigned int _num_models
Number of plastic models for this material.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
static const std::string alpha
virtual void dflowPotential_dstress(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< RankFourTensor > &dr_dstress)
The derivative of the active flow potential(s) with respect to stress.
virtual void hardPotential(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< Real > &h)
The active hardening potentials (one for each internal parameter and for each yield function) by assu...
static const std::complex< double > j(0, 1)
Complex number "j" (also known as "i")
virtual void dyieldFunction_dintnl(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< Real > &df_dintnl)
The derivative of active yield function(s) with respect to their internal parameters (the user object...
virtual void nrStep(const RankTwoTensor &stress, const std::vector< Real > &intnl_old, const std::vector< Real > &intnl, const std::vector< Real > &pm, const RankFourTensor &E_inv, const RankTwoTensor &delta_dp, RankTwoTensor &dstress, std::vector< Real > &dpm, std::vector< Real > &dintnl, const std::vector< bool > &active, std::vector< bool > &deactivated_due_to_ld)
Performs one Newton-Raphson step.
static const std::string k
virtual void calculateConstraints(const RankTwoTensor &stress, const std::vector< Real > &intnl_old, const std::vector< Real > &intnl, const std::vector< Real > &pm, const RankTwoTensor &delta_dp, std::vector< Real > &f, std::vector< RankTwoTensor > &r, RankTwoTensor &epp, std::vector< Real > &ic, const std::vector< bool > &active)
The constraints.
const Elem & get(const ElemType type_in)
std::vector< const SolidMechanicsPlasticModel * > _f
User objects that define the yield functions, flow potentials, etc.
void activeSurfaces(int model, const std::vector< bool > &active, std::vector< unsigned int > &active_surfaces)
Returns the external surface number(s) of the active surfaces of the given model This may be of size=...
Real _min_f_tol
Minimum value of the _f_tol parameters for the Yield Function User Objects.