11 #include "petscblaslapack.h" 12 #include "libmesh/utility.h" 25 params.
addClassDescription(
"Deprecated class: please use CrystalPlasticityKalidindiUpdate and " 26 "ComputeMultipleCrystalPlasticityStress instead. Crystal Plasticity " 27 "base class: FCC system with power law flow rule implemented");
29 params.
addParam<std::vector<Real>>(
"gprops", {},
"Initial values of slip system resistances");
30 params.
addParam<std::vector<Real>>(
"hprops", {},
"Hardening properties");
31 params.
addParam<std::vector<Real>>(
"flowprops", {},
"Parameters used in slip rate equations");
33 "Name of the file containing the slip system");
35 "slip_sys_res_prop_file_name",
37 "Name of the file containing the initial values of slip system resistances");
39 "slip_sys_flow_prop_file_name",
41 "Name of the file containing the values of slip rate equation parameters");
43 "slip_sys_hard_prop_file_name",
45 "Name of the file containing the values of hardness evolution parameters");
46 params.
addParam<
Real>(
"rtol", 1e-6,
"Constitutive stress residue relative tolerance");
47 params.
addParam<
Real>(
"abs_tol", 1e-6,
"Constitutive stress residue absolute tolerance");
48 params.
addParam<
Real>(
"gtol", 1e2,
"Constitutive slip system resistance residual tolerance");
49 params.
addParam<
Real>(
"slip_incr_tol", 2e-2,
"Maximum allowable slip in an increment");
50 params.
addParam<
unsigned int>(
"maxiter", 100,
"Maximum number of iterations for stress update");
52 "maxitergss", 100,
"Maximum number of iterations for slip system resistance update");
54 "num_slip_sys_flowrate_props",
56 "Number of flow rate properties for a slip system");
57 params.
addParam<UserObjectName>(
"read_prop_user_object",
58 "The ElementReadPropertyFile " 59 "GeneralUserObject to read element " 60 "specific property values from file");
61 MooseEnum tan_mod_options(
"exact none",
"none");
64 "Type of tangent moduli for preconditioner: default elastic");
65 MooseEnum intvar_read_options(
"slip_sys_file slip_sys_res_file none",
"none");
69 "Read from options for initial value of internal variables: Default from .i file");
70 params.
addParam<
unsigned int>(
"num_slip_sys_props",
72 "Number of slip system specific properties provided in the file " 73 "containing slip system normals and directions");
75 "gen_random_stress_flag",
77 "Flag to generate random stress to perform time cutback on constitutive failure");
78 params.
addParam<
bool>(
"input_random_scaling_var",
80 "Flag to input scaling variable: _Cijkl(0,0,0,0) when false");
83 "Random scaling variable: Large value can cause non-positive definiteness");
87 "Random integer used to generate random stress when constitutive failure occurs");
89 "maximum_substep_iteration", 1,
"Maximum number of substep iteration");
90 params.
addParam<
bool>(
"use_line_search",
false,
"Use line search in constitutive update");
91 params.
addParam<
Real>(
"min_line_search_step_size", 0.01,
"Minimum line search step size");
92 params.
addParam<
Real>(
"line_search_tol", 0.5,
"Line search bisection method tolerance");
94 "line_search_maxiter", 20,
"Line search bisection method maximum number of iteration");
95 MooseEnum line_search_method(
"CUT_HALF BISECTION",
"CUT_HALF");
97 "line_search_method", line_search_method,
"The method used in line search");
104 _nss(getParam<
int>(
"nss")),
105 _gprops(getParam<
std::vector<
Real>>(
"gprops")),
106 _hprops(getParam<
std::vector<
Real>>(
"hprops")),
107 _flowprops(getParam<
std::vector<
Real>>(
"flowprops")),
108 _slip_sys_file_name(getParam<FileName>(
"slip_sys_file_name")),
109 _slip_sys_res_prop_file_name(getParam<FileName>(
"slip_sys_res_prop_file_name")),
110 _slip_sys_flow_prop_file_name(getParam<FileName>(
"slip_sys_flow_prop_file_name")),
111 _slip_sys_hard_prop_file_name(getParam<FileName>(
"slip_sys_hard_prop_file_name")),
112 _rtol(getParam<
Real>(
"rtol")),
113 _abs_tol(getParam<
Real>(
"abs_tol")),
114 _gtol(getParam<
Real>(
"gtol")),
115 _slip_incr_tol(getParam<
Real>(
"slip_incr_tol")),
116 _maxiter(getParam<unsigned
int>(
"maxiter")),
117 _maxiterg(getParam<unsigned
int>(
"maxitergss")),
118 _num_slip_sys_flowrate_props(getParam<unsigned
int>(
"num_slip_sys_flowrate_props")),
119 _tan_mod_type(getParam<
MooseEnum>(
"tan_mod_type")),
120 _intvar_read_type(getParam<
MooseEnum>(
"intvar_read_type")),
121 _num_slip_sys_props(getParam<unsigned
int>(
"num_slip_sys_props")),
122 _gen_rndm_stress_flag(getParam<bool>(
"gen_random_stress_flag")),
123 _input_rndm_scale_var(getParam<bool>(
"input_random_scaling_var")),
124 _rndm_scale_var(getParam<
Real>(
"random_scaling_var")),
125 _rndm_seed(getParam<unsigned
int>(
"random_seed")),
126 _max_substep_iter(getParam<unsigned
int>(
"maximum_substep_iteration")),
127 _use_line_search(getParam<bool>(
"use_line_search")),
128 _min_lsrch_step(getParam<
Real>(
"min_line_search_step_size")),
129 _lsrch_tol(getParam<
Real>(
"line_search_tol")),
130 _lsrch_max_iter(getParam<unsigned
int>(
"line_search_maxiter")),
131 _lsrch_method(getParam<
MooseEnum>(
"line_search_method")),
141 _gss(declareProperty<
std::vector<
Real>>(
"gss")),
142 _gss_old(getMaterialPropertyOld<
std::vector<
Real>>(
144 _acc_slip(declareProperty<
Real>(
"acc_slip")),
146 getMaterialPropertyOld<
Real>(
"acc_slip")),
149 _deformation_gradient(getMaterialProperty<
RankTwoTensor>(
"deformation_gradient")),
150 _deformation_gradient_old(getMaterialPropertyOld<
RankTwoTensor>(
"deformation_gradient")),
151 _elasticity_tensor_name(_base_name +
"elasticity_tensor"),
152 _elasticity_tensor(getMaterialPropertyByName<
RankFourTensor>(_elasticity_tensor_name)),
154 _mo(_nss * LIBMESH_DIM),
155 _no(_nss * LIBMESH_DIM),
162 _dgss_dsliprate(_nss, _nss)
182 mooseError(
"Crystal Plasticity Error: Specify number of internal variable's initial values to " 183 "be read from slip system file");
238 for (
unsigned int i = 0; i <
_nss; ++i)
253 for (
unsigned int i = 0; i <
_nss; ++i)
255 mooseError(
"Error FiniteStrainCrystalPlasticity: Premature end of slip_sys_res_prop file");
265 mooseError(
"FiniteStrainCrystalPLasticity: Error in reading slip system resistance properties: " 266 "Specify input in .i file or in slip_sys_res_prop_file or in slip_sys_file");
270 unsigned int num_data_grp = 3;
273 for (
unsigned int i = 0; i <
_gprops.size() / num_data_grp; ++i)
278 vs =
_gprops[i * num_data_grp];
279 ve =
_gprops[i * num_data_grp + 1];
281 if (vs <= 0 || ve <= 0)
282 mooseError(
"FiniteStrainCrystalPLasticity: Indices in gss property read must be positive " 288 if (vs != floor(vs) || ve != floor(ve))
289 mooseError(
"FiniteStrainCrystalPLasticity: Error in reading slip system resistances: Values " 290 "specifying start and end number of slip system groups should be integer");
292 is =
static_cast<unsigned int>(vs);
293 ie =
static_cast<unsigned int>(ve);
296 mooseError(
"FiniteStrainCrystalPLasticity: Start index is = ",
298 " should be greater than end index ie = ",
300 " in slip system resistance property read");
302 for (
unsigned int j =
is;
j <= ie; ++
j)
306 for (
unsigned int i = 0; i <
_nss; ++i)
308 mooseError(
"FiniteStrainCrystalPLasticity: Value of resistance for slip system ",
325 std::vector<Real> vec;
328 for (
unsigned int i = 0; i <
_nss; ++i)
331 if (!(file >> vec[
j]))
333 "Error FiniteStrainCrystalPlasticity: Premature end of slip_sys_flow_rate_param file");
347 mooseError(
"FiniteStrainCrystalPLasticity: Error in reading flow rate properties: Specify " 348 "input in .i file or a slip_sys_flow_prop_file_name");
357 for (
unsigned int i = 0; i <
_flowprops.size() / num_data_grp; ++i)
365 if (vs <= 0 || ve <= 0)
366 mooseError(
"FiniteStrainCrystalPLasticity: Indices in flow rate parameter read must be " 367 "positive integers: is = ",
372 if (vs != floor(vs) || ve != floor(ve))
373 mooseError(
"FiniteStrainCrystalPLasticity: Error in reading flow props: Values specifying " 374 "start and end number of slip system groups should be integer");
376 is =
static_cast<unsigned int>(vs);
377 ie =
static_cast<unsigned int>(ve);
380 mooseError(
"FiniteStrainCrystalPLasticity: Start index is = ",
382 " should be greater than end index ie = ",
384 " in flow rate parameter read");
386 for (
unsigned int j =
is;
j <= ie; ++
j)
393 for (
unsigned int i = 0; i <
_nss; ++i)
395 if (!(
_a0(i) > 0.0 &&
_xm(i) > 0.0))
398 "FiniteStrainCrystalPlasticity: Non-positive flow rate parameters ",
_a0(i),
",",
_xm(i));
415 mooseError(
"FiniteStrainCrystalPLasticity: Error in reading hardness properties: Specify input " 416 "in .i file or a slip_sys_hard_prop_file_name");
428 Real vec[LIBMESH_DIM];
429 std::ifstream fileslipsys;
435 for (
unsigned int i = 0; i <
_nss; ++i)
439 if (!(fileslipsys >> vec[
j]))
440 mooseError(
"Crystal Plasticity Error: Premature end of file reading slip system file \n");
444 mag = Utility::pow<2>(vec[0]) + Utility::pow<2>(vec[1]) + Utility::pow<2>(vec[2]);
447 for (
unsigned j = 0;
j < LIBMESH_DIM; ++
j)
448 _no(i * LIBMESH_DIM +
j) = vec[
j] / mag;
452 if (!(fileslipsys >> vec[
j]))
453 mooseError(
"Crystal Plasticity Error: Premature end of file reading slip system file \n");
456 mag = Utility::pow<2>(vec[0]) + Utility::pow<2>(vec[1]) + Utility::pow<2>(vec[2]);
460 _mo(i * LIBMESH_DIM +
j) = vec[
j] / mag;
464 mag +=
_mo(i * LIBMESH_DIM +
j) *
_no(i * LIBMESH_DIM +
j);
468 "Crystal Plasicity Error: Slip direction and normal not orthonormal, System number = ",
475 mooseError(
"Crystal Plasticity Error: Premature end of file reading slip system file - " 476 "check in slip system file read input options/values\n");
495 unsigned int substep_iter = 1;
496 unsigned int num_substep = 1;
513 _dt = dt_original / num_substep;
515 for (
unsigned int istep = 0; istep < num_substep; ++istep)
522 if (istep == num_substep - 1)
543 mooseWarning(
"FiniteStrainCrystalPlasticity: Failure with substepping");
601 mooseError(
"FiniteStrainCrystalPlasticity: Constitutive failure");
645 std::vector<Real> gss_prev(
_nss);
663 for (
unsigned i = 0; i <
_nss; ++i)
676 mooseWarning(
"FiniteStrainCrystalPLasticity: Hardness Integration error gmax", gmax,
"\n");
733 unsigned int iter = 0;
736 Real rnorm, rnorm0, rnorm_prev;
743 "FiniteStrainCrystalPLasticity: Slip increment exceeds tolerance - Element number ",
766 "FiniteStrainCrystalPLasticity: Slip increment exceeds tolerance - Element number ",
780 mooseWarning(
"FiniteStrainCrystalPLasticity: Failed with line search");
795 mooseWarning(
"FiniteStrainCrystalPLasticity: Stress Integration error rmax = ", rnorm);
838 DenseVector<Real> hb(
_nss);
844 for (
unsigned int i = 0; i <
_nss; ++i)
850 for (
unsigned int i = 0; i <
_nss; ++i)
855 for (
unsigned int i = 0; i <
_nss; ++i)
862 for (
unsigned int j = 0;
j <
_nss; ++
j)
864 unsigned int iplane, jplane;
868 if (iplane == jplane)
898 ce_pk2 = ce_pk2 /
_fe.
det();
901 for (
unsigned int i = 0; i <
_nss; ++i)
909 eqv_slip_incr.
zero();
910 for (
unsigned int i = 0; i <
_nss; ++i)
913 eqv_slip_incr = iden - eqv_slip_incr;
931 std::vector<RankTwoTensor> dtaudpk2(
_nss), dfpinvdslip(
_nss);
933 for (
unsigned int i = 0; i <
_nss; ++i)
935 dtaudpk2[i] =
_s0[i];
948 deedfe(i,
j,
k, i) = deedfe(i,
j,
k, i) +
_fe(
k,
j) * 0.5;
949 deedfe(i,
j,
k,
j) = deedfe(i,
j,
k,
j) +
_fe(
k, i) * 0.5;
952 for (
unsigned int i = 0; i <
_nss; ++i)
953 dfpinvdpk2 += (dfpinvdslip[i] *
_dslipdtau(i)).outerProduct(dtaudpk2[i]);
963 for (
unsigned int i = 0; i <
_nss; ++i)
966 std::copysign(1.0,
_tau(i)) *
_dt;
977 for (
unsigned int i = 0; i <
_nss; ++i)
996 PetscScalar cmat[LIBMESH_DIM][LIBMESH_DIM], work[10];
997 PetscReal w[LIBMESH_DIM];
998 PetscBLASInt nd = LIBMESH_DIM, lwork = 10,
info;
1000 c =
a.transpose() *
a;
1004 cmat[i][
j] =
c(i,
j);
1006 LAPACKsyev_(
"V",
"U", &nd, &cmat[0][0], &nd, w, work, &lwork, &
info);
1009 mooseError(
"FiniteStrainCrystalPLasticity: DSYEV function call in getMatRot function failed");
1018 evec(i,
j) = cmat[i][
j];
1051 DenseVector<Real> mo(LIBMESH_DIM *
_nss), no(LIBMESH_DIM *
_nss);
1054 for (
unsigned int i = 0; i <
_nss; ++i)
1058 mo(i * LIBMESH_DIM +
j) = 0.0;
1060 mo(i * LIBMESH_DIM +
j) =
1066 no(i * LIBMESH_DIM +
j) = 0.0;
1068 no(i * LIBMESH_DIM +
j) =
1074 for (
unsigned int i = 0; i <
_nss; ++i)
1077 _s0[i](
j,
k) = mo(i * LIBMESH_DIM +
j) * no(i * LIBMESH_DIM +
k);
1093 deedfe(i,
j,
k, i) = deedfe(i,
j,
k, i) +
_fe(
k,
j) * 0.5;
1094 deedfe(i,
j,
k,
j) = deedfe(i,
j,
k,
j) +
_fe(
k, i) * 0.5;
1097 usingTensorIndices(i_, j_, k_, l_);
1107 tan_mod(i,
j, i, l) = tan_mod(i,
j, i, l) + pk2fet(l,
j);
1108 tan_mod(i,
j,
j, l) = tan_mod(i,
j,
j, l) + fepk2(i, l);
1111 tan_mod += dsigdpk2dfe;
1152 unsigned int count = 0;
1157 Real rnorm = 1000.0;
1169 if ((rnorm1 / rnorm0) <
_lsrch_tol || s_a * s_b > 0)
1178 step = 0.5 * (step_b + step_a);
1184 if (s_m * s_a < 0.0)
1189 if (s_m * s_b < 0.0)
1204 mooseError(
"Line search meothod is not provided.");
const MaterialProperty< RankTwoTensor > & _pk2_old
void internalVariableUpdateNRiteration()
This function updates internal variables after each NewTon Raphson iteration (_fp_inv) ...
MaterialProperty< RankFourTensor > & _Jacobian_mult
derivative of stress w.r.t. strain (_dstress_dstrain)
MaterialProperty< RankTwoTensor > & _lag_e
unsigned int _rndm_seed
Seed value.
RankTwoTensorTempl< Real > inverse() const
bool _read_from_slip_sys_file
const MaterialProperty< RankTwoTensor > & _deformation_gradient
std::vector< Real > _gss_tmp_old
ComputeStressBase is the base class for stress tensors computed from MOOSE's strain calculators...
bool _input_rndm_scale_var
Input option for scaling variable to generate random stress when convergence fails.
static void initRandom(unsigned int)
virtual void getHardnessParams()
This function assign flow rate parameters from .i file - see test.
bool _use_line_search
Flag to activate line serach.
const MaterialProperty< std::vector< Real > > & _gss_old
static RankFourTensorTempl< Real > IdentityFour()
const MaterialProperty< RankFourTensor > & _elasticity_tensor
Elasticity tensor material property.
virtual void calcJacobian(RankFourTensor &)
This function calculate jacobian.
static InputParameters validParams()
virtual void computeQpStress()
This function updates the stress at a quadrature point.
MaterialProperty< RankTwoTensor > & _pk2
virtual void initQpStatefulProperties()
This function initializes the stateful properties such as stress, plastic deformation gradient...
RankTwoTensor _fp_prev_inv
virtual void solveStress()
This function solves for stress, updates plastic deformation gradient.
std::string _slip_sys_flow_prop_file_name
File should contain values of the flow rate equation parameters.
virtual void solveStatevar()
This function solves internal variables.
FiniteStrainCrystalPlasticity(const InputParameters ¶meters)
std::vector< Real > _flowprops
virtual void update_slip_system_resistance()
This function updates the slip system resistances.
static constexpr std::size_t dim
unsigned int _num_slip_sys_flowrate_props
Number of slip system flow rate parameters.
virtual void calcResidual(RankTwoTensor &)
This function calculate stress residual.
const MaterialProperty< RankTwoTensor > & _fp_old
MaterialProperty< Real > & _acc_slip
RankTwoTensor getMatRot(const RankTwoTensor &a)
This function perform RU decomposition to obtain the rotation tensor.
ADRealEigenVector< T, D, asd > sqrt(const ADRealEigenVector< T, D, asd > &)
virtual void readFileHardnessParams()
This function read hardness parameters from file.
virtual void updateGss()
This function updates the slip system resistances.
Real _rtol
Stress residual equation relative tolerance.
void mooseWarning(Args &&... args) const
std::string _slip_sys_hard_prop_file_name
The hardening parameters in this class are read from .i file. The user can override to read from file...
unsigned int _maxiter
Maximum number of iterations for stress update.
virtual void resize(const std::size_t size) override final
virtual void assignSlipSysRes()
This function assign initial values of slip system resistances/internal variables read from getSlipSy...
Real _dfgrd_scale_factor
Scales the substepping increment to obtain deformation gradient at a substep iteration.
const unsigned int _nss
Number of slip system resistance.
RankTwoTensor _pk2_tmp_old
ADRealEigenVector< T, D, asd > abs(const ADRealEigenVector< T, D, asd > &)
virtual RankFourTensor elastoPlasticTangentModuli()
This function calculate the exact tangent moduli for preconditioner.
RankTwoTensor _delta_dfgrd
Flag to check whether convergence is achieved.
virtual RankFourTensor elasticTangentModuli()
This function calculate the elastic tangent moduli for preconditioner.
MaterialProperty< std::vector< Real > > & _gss
void addIa(const Real &a)
bool checkFileReadable(const std::string &filename, bool check_line_endings=false, bool throw_on_unreadable=true, bool check_for_git_lfs_pointer=true)
static RankTwoTensorTempl< Real > genRandomSymmTensor(Real stddev, Real mean)
virtual RankFourTensor calcTangentModuli()
This function calculate the tangent moduli for preconditioner.
bool line_search_update(const Real rnorm_prev, const RankTwoTensor)
This function performs the line search update.
const MaterialProperty< RankTwoTensor > & _deformation_gradient_old
Real doubleContraction(const RankTwoTensorTempl< Real > &a) const
const MaterialProperty< RankTwoTensor > & _crysrot
void calc_schmid_tensor()
This function calculate the Schmid tensor.
DenseVector< Real > _slip_sys_props
PetscErrorCode PetscInt const PetscInt IS * is
MaterialProperty< RankTwoTensor > & _fp
virtual void postSolveStatevar()
This function update internal variable after solve.
unsigned int _lsrch_max_iter
Line search bisection method maximum iteration number.
virtual void getFlowRateParams()
This function assign flow rate parameters - see test.
static InputParameters validParams()
Real _gtol
Internal variable update equation tolerance.
DenseVector< Real > _dslipdtau
std::string _slip_sys_res_prop_file_name
File should contain initial values of the slip system resistances.
virtual void postSolveQp()
This function update stress and internal variable after solve.
virtual void preSolveStress()
This function set variables for stress solve.
bool _gen_rndm_stress_flag
RankTwoTensorTempl< Real > transpose() const
std::vector< Real > _gss_tmp
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual void postSolveStress()
This function update stress and plastic deformation gradient after solve.
virtual void calc_resid_jacob(RankTwoTensor &, RankFourTensor &)
This function calls the residual and jacobian functions used in the stress update algorithm...
MooseEnum _tan_mod_type
Type of tangent moduli calculation.
RankFourTensorTempl< Real > times(const RankTwoTensorTempl< Real > &b) const
Real _lsrch_tol
Line search bisection method tolerance.
virtual void computeQpElasticityTensor()
This function updates the elasticity tensor at a quadrature point.
IntRange< T > make_range(T beg, T end)
std::vector< RankTwoTensor > _s0
DenseMatrix< Real > _dgss_dsliprate
virtual void getSlipSystems()
This function reads slip system from file - see test.
RankTwoTensor _fp_old_inv
void mooseError(Args &&... args) const
unsigned int _num_slip_sys_props
Number of slip system specific properties provided in the file containing slip system normals and dir...
std::vector< Real > _hprops
virtual void readFileFlowRateParams()
This function read flow rate parameters from file - see test.
unsigned int _max_substep_iter
Maximum number of substep iterations.
static const std::complex< double > j(0, 1)
Complex number "j" (also known as "i")
MaterialProperty< RankTwoTensor > & _stress
Stress material property.
Real _rndm_scale_var
Scaling value.
bool _first_step_iter
Flags to reset variables and reinitialize variables.
virtual void initSlipSysProps()
This function initializes slip system resistances.
DenseVector< Real > _slip_incr
FiniteStrainCrystalPlasticity uses the multiplicative decomposition of deformation gradient and solve...
registerMooseObjectDeprecated("SolidMechanicsApp", FiniteStrainCrystalPlasticity, "11/15/2024 12:00")
virtual void getInitSlipSysRes()
This function assign slip system resistances - see test.
virtual void initAdditionalProps()
This function initializes additional parameters.
virtual void preSolveStatevar()
This function set variables for internal variable solve.
RankTwoTensor _dfgrd_tmp_old
MaterialProperty< RankTwoTensor > & _update_rot
Real _abs_tol
Stress residual equation absolute tolerance.
std::vector< Real > _gprops
virtual void readFileInitSlipSysRes()
This function read slip system resistances from file - see test.
unsigned int _maxiterg
Maximum number of iterations for internal variable update.
RankFourTensorTempl< T > invSymm() const
virtual void getSlipIncrements()
This function updates the slip increments.
MooseUnits pow(const MooseUnits &, int)
const MaterialProperty< Real > & _acc_slip_old
static const std::string k
RankTwoTensor get_current_rotation(const RankTwoTensor &a)
This function perform RU decomposition to obtain the rotation tensor.
void ErrorVector unsigned int
virtual void solveQp()
This function solves stress and internal variables.
const Elem *const & _current_elem
MooseEnum _intvar_read_type
Read from options for initial values of internal variables.
Real _min_lsrch_step
Minimum line search step size.
virtual void preSolveQp()
This function set variables for stress and internal variable solve.
std::string _slip_sys_file_name
File should contain slip plane normal and direction. See test.
Real _slip_incr_tol
Slip increment tolerance.