www.mooseframework.org
FiniteStrainCPSlipRateRes.C
Go to the documentation of this file.
1 /****************************************************************/
2 /* MOOSE - Multiphysics Object Oriented Simulation Environment */
3 /* */
4 /* All contents are licensed under LGPL V2.1 */
5 /* See LICENSE for full restrictions */
6 /****************************************************************/
7 
9 #include "libmesh/utility.h"
10 
11 template <>
12 InputParameters
14 {
15  InputParameters params = validParams<FiniteStrainCrystalPlasticity>();
16  return params;
17 }
18 
19 FiniteStrainCPSlipRateRes::FiniteStrainCPSlipRateRes(const InputParameters & parameters)
20  : FiniteStrainCrystalPlasticity(parameters),
21  _resid(_nss),
22  _slip_rate(_nss),
23  _dsliprate_dgss(_nss),
24  _jacob(_nss, _nss),
25  _dsliprate_dsliprate(_nss, _nss)
26 {
27 }
28 
29 void
31 {
33  solveStress();
34  if (_err_tol)
35  return;
37 }
38 
39 void
41 {
43  _slip_rate.zero();
44 }
45 
46 void
48 {
49  Real rnorm, rnorm0, rnorm_prev;
50  unsigned int iter = 0;
51 
52 #ifdef DEBUG
53  std::vector<Real> rnormst(_maxiter + 1), slipratest(_maxiter + 1); // Use for Debugging
54 #endif
55 
57  if (_err_tol)
58  return;
59  rnorm = calcResidNorm();
60  rnorm0 = rnorm;
61 
62 #ifdef DEBUG
63  rnormst[iter] = rnorm;
64  Real slipratemax = 0.0;
65  for (unsigned int i = 0; i < _nss; ++i)
66  if (std::abs(_slip_rate(i)) > slipratemax)
67  slipratemax = std::abs(_slip_rate(i));
68  slipratest[iter] = slipratemax;
69 #endif
70 
71  while (rnorm > _rtol * rnorm0 && rnorm0 > _abs_tol && iter < _maxiter)
72  {
73  calcUpdate();
74 
75  DenseVector<Real> update = _resid;
76 
77  _slip_rate -= update;
78 
80  if (_err_tol)
81  return;
82  rnorm_prev = rnorm;
83  rnorm = calcResidNorm();
84 
85  if (_use_line_search && rnorm > rnorm_prev && !lineSearchUpdateSlipRate(rnorm_prev, update))
86  {
87 #ifdef DEBUG
88  mooseWarning("FiniteStrainCrystalPLasticity: Failed with line search");
89 #endif
90  _err_tol = true;
91  return;
92  }
93 
95 
96  if (_use_line_search)
97  rnorm = calcResidNorm();
98  iter++;
99 
100 #ifdef DEBUG
101  slipratemax = 0.0;
102  for (unsigned int i = 0; i < _nss; ++i)
103  if (std::abs(_slip_rate(i)) > slipratemax)
104  slipratemax = std::abs(_slip_rate(i));
105  rnormst[iter] = rnorm;
106  slipratest[iter] = slipratemax;
107 #endif
108  }
109 
110  if (iter == _maxiter)
111  {
112 #ifdef DEBUG
113  mooseWarning("FiniteStrainCPSlipRateRes: NR exceeds maximum iteration ", iter, " ", rnorm);
114 #endif
115  _err_tol = true;
116  return;
117  }
118 }
119 
120 void
122 {
124  if (_err_tol)
125  return;
127 }
128 
129 void
131 {
132  RankTwoTensor eqv_slip_incr, ce, ee;
133  const RankTwoTensor iden(RankTwoTensor::initIdentity);
134 
136  _slip_incr *= _dt;
137 
138  for (unsigned int i = 0; i < _nss; ++i)
139  eqv_slip_incr += _s0[i] * _slip_incr(i);
140 
141  eqv_slip_incr = iden - eqv_slip_incr;
142 
143  _fp_inv = _fp_old_inv * eqv_slip_incr;
144  _fe = _dfgrd_tmp * _fp_inv;
145 
146  ce = _fe.transpose() * _fe;
147  ee = ce - iden;
148  ee *= 0.5;
149 
150  _pk2_tmp = _elasticity_tensor[_qp] * ee;
151 
152  for (unsigned int i = 0; i < _nss; ++i)
153  _tau(i) = _pk2_tmp.doubleContraction(_s0[i]);
154 
157 
158  if (_err_tol)
159  return;
160 
161  for (unsigned int i = 0; i < _nss; ++i)
162  _resid(i) = _slip_rate(i) - _slip_incr(i) / _dt;
163 }
164 
165 void
167 {
168  //_dsliprate_dsliprate not reinitialized to zero, hence order is important
171 
172  for (unsigned int i = 0; i < _nss; ++i)
173  for (unsigned int j = 0; j < _nss; ++j)
174  {
175  _jacob(i, j) = 0.0;
176  if (i == j)
177  _jacob(i, j) += 1.0;
178  _jacob(i, j) -= _dsliprate_dsliprate(i, j);
179  }
180 }
181 
182 void
184 {
185  RankFourTensor dfedfpinv, deedfe, dfpinvdpk2;
186  std::vector<RankTwoTensor> dtaudpk2(_nss), dfpinvdsliprate(_nss);
187 
188  for (unsigned int i = 0; i < _nss; ++i)
189  {
190  dtaudpk2[i] = _s0[i];
191  dfpinvdsliprate[i] = -_fp_old_inv * _s0[i] * _dt;
192  }
193 
194  for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
195  for (unsigned int j = 0; j < LIBMESH_DIM; ++j)
196  for (unsigned int k = 0; k < LIBMESH_DIM; ++k)
197  dfedfpinv(i, j, k, j) = _dfgrd_tmp(i, k);
198 
199  for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
200  for (unsigned int j = 0; j < LIBMESH_DIM; ++j)
201  for (unsigned int k = 0; k < LIBMESH_DIM; ++k)
202  {
203  deedfe(i, j, k, i) = deedfe(i, j, k, i) + _fe(k, j) * 0.5;
204  deedfe(i, j, k, j) = deedfe(i, j, k, j) + _fe(k, i) * 0.5;
205  }
206 
207  RankFourTensor dpk2dfpinv;
208 
209  dpk2dfpinv = _elasticity_tensor[_qp] * deedfe * dfedfpinv;
210 
211  for (unsigned int i = 0; i < _nss; ++i)
212  for (unsigned int j = 0; j < _nss; ++j)
213  _dsliprate_dsliprate(i, j) =
214  _dslipdtau(i) * dtaudpk2[i].doubleContraction(dpk2dfpinv * dfpinvdsliprate[j]);
215 }
216 
217 void
219 {
220  for (unsigned int i = 0; i < _nss; ++i)
221  for (unsigned int j = 0; j < _nss; ++j)
223 }
224 
225 void
227 {
229 
230  if (_err_tol)
231  return;
232 
233  _dslipdtau *= 1.0 / _dt;
234 
235  for (unsigned int i = 0; i < _nss; ++i)
236  _dsliprate_dgss(i) = -_a0(i) / _xm(i) *
237  std::pow(std::abs(_tau(i) / _gss_tmp[i]), 1.0 / _xm(i) - 1.0) * _tau(i) /
238  std::pow(_gss_tmp[i], 2.0);
239 }
240 
241 void
243 {
244  DenseMatrix<Real> A = _jacob;
245  DenseVector<Real> r(_nss);
246  DenseVector<Real> x(_nss);
247 
248  r = _resid;
249 
250  A.lu_solve(r, x);
251 
252  _resid = x;
253 }
254 
255 Real
257 {
258  Real rnorm = 0.0;
259  for (unsigned int i = 0; i < _nss; ++i)
260  rnorm += Utility::pow<2>(_resid(i));
261  rnorm = std::sqrt(rnorm) / _nss;
262 
263  return rnorm;
264 }
265 
266 bool
268  const DenseVector<Real> & update)
269 {
270  if (_lsrch_method == "CUT_HALF")
271  {
272  Real rnorm;
273  Real step = 1.0;
274  do
275  {
276  for (unsigned int i = 0; i < update.size(); ++i)
277  _slip_rate(i) += step * update(i);
278 
279  step /= 2.0;
280 
281  for (unsigned int i = 0; i < update.size(); ++i)
282  _slip_rate(i) -= step * update(i);
283 
285  if (_err_tol)
286  return false;
287  rnorm = calcResidNorm();
288  } while (rnorm > rnorm_prev && step > _min_lsrch_step);
289 
290  if (rnorm > rnorm_prev && step <= _min_lsrch_step)
291  return false;
292 
293  return true;
294  }
295  else if (_lsrch_method == "BISECTION")
296  {
297  unsigned int count = 0;
298  Real step_a = 0.0;
299  Real step_b = 1.0;
300  Real step = 1.0;
301  Real s_m = 1000.0;
302  Real rnorm = 1000.0;
303 
304  Real s_b = calcResidDotProdUpdate(update);
305  Real rnorm1 = calcResidNorm();
306 
307  for (unsigned int i = 0; i < update.size(); ++i)
308  _slip_rate(i) += update(i);
309 
311  Real s_a = calcResidDotProdUpdate(update);
312  Real rnorm0 = calcResidNorm();
313 
314  for (unsigned int i = 0; i < update.size(); ++i)
315  _slip_rate(i) -= update(i);
316 
317  if ((rnorm1 / rnorm0) < _lsrch_tol || s_a * s_b > 0)
318  {
320  return true;
321  }
322 
323  while ((rnorm / rnorm0) > _lsrch_tol && count < _lsrch_max_iter)
324  {
325 
326  for (unsigned int i = 0; i < update.size(); ++i)
327  _slip_rate(i) += step * update(i);
328 
329  step = 0.5 * (step_a + step_b);
330 
331  for (unsigned int i = 0; i < update.size(); ++i)
332  _slip_rate(i) -= step * update(i);
333 
335  s_m = calcResidDotProdUpdate(update);
336  rnorm = calcResidNorm();
337 
338  if (s_m * s_a < 0.0)
339  {
340  step_b = step;
341  s_b = s_m;
342  }
343  if (s_m * s_b < 0.0)
344  {
345  step_a = step;
346  s_a = s_m;
347  }
348  count++;
349  }
350 
351  if ((rnorm / rnorm0) < _lsrch_tol && count < _lsrch_max_iter)
352  return true;
353 
354  return false;
355  }
356  else
357  {
358  mooseError("Line search meothod is not provided.");
359  return false;
360  }
361 }
362 
363 Real
364 FiniteStrainCPSlipRateRes::calcResidDotProdUpdate(const DenseVector<Real> & update)
365 {
366  Real dotprod = 0.0;
367  for (unsigned int i = 0; i < _nss; ++i)
368  dotprod += _resid(i) * update(i);
369  return dotprod;
370 }
bool _use_line_search
Flag to activate line serach.
virtual void solveStress()
This function solves for stress, updates plastic deformation gradient.
const MaterialProperty< RankFourTensor > & _elasticity_tensor
virtual void getSlipIncrements()
This function updates the slip system resistances.
Real calcResidDotProdUpdate(const DenseVector< Real > &)
This function calculates the dot product of residual and update.
virtual void update_slip_system_resistance()
This function updates the slip system resistances.
virtual void preSolveStress()
This function sets variable for internal variable solve.
Real _rtol
Stress residual equation relative tolerance.
unsigned int _maxiter
Maximum number of iterations for stress update.
const unsigned int _nss
Number of slip system resistance.
void calcUpdate()
This function calculates and updates the residual of slip rate.
FiniteStrainCPSlipRateRes(const InputParameters &parameters)
DenseMatrix< Real > _dsliprate_dsliprate
unsigned int _lsrch_max_iter
Line search bisection method maximum iteration number.
virtual Real calcResidNorm()
This function calculates the residual norm.
virtual void calcResidJacobSlipRate()
This function calculates residual and jacobian of slip rate.
virtual void preSolveStress()
This function set variables for stress solve.
virtual void calcJacobianSlipRate()
This function calculates jacobian of slip rate.
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)
InputParameters validParams< FiniteStrainCrystalPlasticity >()
virtual void calcDtauDsliprate()
This function calculates partial derivative of resolved shear stress with respect to split rate...
virtual void postSolveStress()
This function update stress and plastic deformation gradient after solve.
Real _lsrch_tol
Line search bisection method tolerance.
virtual void calcDgssDsliprate()
This function calculates partial derivative of slip system resistances with respect to split rate...
virtual void calcResidualSlipRate()
This function calculates residual of slip rate.
bool lineSearchUpdateSlipRate(const Real, const DenseVector< Real > &)
This function performs the line search update.
InputParameters validParams< FiniteStrainCPSlipRateRes >()
virtual void solveStatevar()
This function solves internal variables.
Real _abs_tol
Stress residual equation absolute tolerance.
virtual void getSlipIncrements()
This function updates the slip increments.
Real _min_lsrch_step
Minimum line search step size.