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 "FiniteStrainCPSlipRateRes.h"
11 : #include "libmesh/utility.h"
12 :
13 : registerMooseObjectDeprecated("TensorMechanicsApp", FiniteStrainCPSlipRateRes, "11/15/2024 12:00");
14 :
15 : InputParameters
16 36 : FiniteStrainCPSlipRateRes::validParams()
17 : {
18 36 : InputParameters params = FiniteStrainCrystalPlasticity::validParams();
19 36 : params.addClassDescription("Deprecated class: please use CrystalPlasticityKalidindiUpdate and "
20 : "ComputeMultipleCrystalPlasticityStress instead.");
21 36 : return params;
22 0 : }
23 :
24 27 : FiniteStrainCPSlipRateRes::FiniteStrainCPSlipRateRes(const InputParameters & parameters)
25 : : FiniteStrainCrystalPlasticity(parameters),
26 27 : _resid(_nss),
27 27 : _slip_rate(_nss),
28 27 : _dsliprate_dgss(_nss),
29 27 : _jacob(_nss, _nss),
30 54 : _dsliprate_dsliprate(_nss, _nss)
31 : {
32 27 : }
33 :
34 : void
35 42456 : FiniteStrainCPSlipRateRes::solveStatevar()
36 : {
37 42456 : preSolveStress();
38 42456 : solveStress();
39 42456 : if (_err_tol)
40 : return;
41 32656 : postSolveStress();
42 : }
43 :
44 : void
45 42456 : FiniteStrainCPSlipRateRes::preSolveStress()
46 : {
47 42456 : FiniteStrainCrystalPlasticity::preSolveStress();
48 : _slip_rate.zero();
49 42456 : }
50 :
51 : void
52 42456 : FiniteStrainCPSlipRateRes::solveStress()
53 : {
54 : Real rnorm, rnorm0, rnorm_prev;
55 : unsigned int iter = 0;
56 :
57 : #ifdef DEBUG
58 : std::vector<Real> rnormst(_maxiter + 1), slipratest(_maxiter + 1); // Use for Debugging
59 : #endif
60 :
61 42456 : calcResidJacobSlipRate();
62 42456 : if (_err_tol)
63 : return;
64 32656 : rnorm = calcResidNorm();
65 : rnorm0 = rnorm;
66 :
67 : #ifdef DEBUG
68 : rnormst[iter] = rnorm;
69 : Real slipratemax = 0.0;
70 : for (unsigned int i = 0; i < _nss; ++i)
71 : if (std::abs(_slip_rate(i)) > slipratemax)
72 : slipratemax = std::abs(_slip_rate(i));
73 : slipratest[iter] = slipratemax;
74 : #endif
75 :
76 310384 : while (rnorm > _rtol * rnorm0 && rnorm0 > _abs_tol && iter < _maxiter)
77 : {
78 277728 : calcUpdate();
79 :
80 : DenseVector<Real> update = _resid;
81 :
82 277728 : _slip_rate -= update;
83 :
84 277728 : calcResidualSlipRate();
85 277728 : if (_err_tol)
86 : return;
87 : rnorm_prev = rnorm;
88 277728 : rnorm = calcResidNorm();
89 :
90 277728 : if (_use_line_search && rnorm > rnorm_prev && !lineSearchUpdateSlipRate(rnorm_prev, update))
91 : {
92 : #ifdef DEBUG
93 : mooseWarning("FiniteStrainCrystalPLasticity: Failed with line search");
94 : #endif
95 0 : _err_tol = true;
96 0 : return;
97 : }
98 :
99 277728 : calcJacobianSlipRate();
100 :
101 277728 : if (_use_line_search)
102 139184 : rnorm = calcResidNorm();
103 277728 : iter++;
104 :
105 : #ifdef DEBUG
106 : slipratemax = 0.0;
107 : for (unsigned int i = 0; i < _nss; ++i)
108 : if (std::abs(_slip_rate(i)) > slipratemax)
109 : slipratemax = std::abs(_slip_rate(i));
110 : rnormst[iter] = rnorm;
111 : slipratest[iter] = slipratemax;
112 : #endif
113 : }
114 :
115 32656 : if (iter == _maxiter)
116 : {
117 : #ifdef DEBUG
118 : mooseWarning("FiniteStrainCPSlipRateRes: NR exceeds maximum iteration ", iter, " ", rnorm);
119 : #endif
120 0 : _err_tol = true;
121 0 : return;
122 : }
123 : }
124 :
125 : void
126 42456 : FiniteStrainCPSlipRateRes::calcResidJacobSlipRate()
127 : {
128 42456 : calcResidualSlipRate();
129 42456 : if (_err_tol)
130 : return;
131 32656 : calcJacobianSlipRate();
132 : }
133 :
134 : void
135 320184 : FiniteStrainCPSlipRateRes::calcResidualSlipRate()
136 : {
137 320184 : RankTwoTensor eqv_slip_incr, ce, ee;
138 320184 : const RankTwoTensor iden(RankTwoTensor::initIdentity);
139 :
140 : _slip_incr = _slip_rate;
141 320184 : _slip_incr *= _dt;
142 :
143 4162392 : for (unsigned int i = 0; i < _nss; ++i)
144 3842208 : eqv_slip_incr += _s0[i] * _slip_incr(i);
145 :
146 320184 : eqv_slip_incr = iden - eqv_slip_incr;
147 :
148 320184 : _fp_inv = _fp_old_inv * eqv_slip_incr;
149 320184 : _fe = _dfgrd_tmp * _fp_inv;
150 :
151 320184 : ce = _fe.transpose() * _fe;
152 320184 : ee = ce - iden;
153 320184 : ee *= 0.5;
154 :
155 320184 : _pk2_tmp = _elasticity_tensor[_qp] * ee;
156 :
157 4162392 : for (unsigned int i = 0; i < _nss; ++i)
158 3842208 : _tau(i) = _pk2_tmp.doubleContraction(_s0[i]);
159 :
160 320184 : update_slip_system_resistance();
161 320184 : getSlipIncrements();
162 :
163 320184 : if (_err_tol)
164 : return;
165 :
166 4034992 : for (unsigned int i = 0; i < _nss; ++i)
167 3724608 : _resid(i) = _slip_rate(i) - _slip_incr(i) / _dt;
168 : }
169 :
170 : void
171 310384 : FiniteStrainCPSlipRateRes::calcJacobianSlipRate()
172 : {
173 : //_dsliprate_dsliprate not reinitialized to zero, hence order is important
174 310384 : calcDtauDsliprate();
175 310384 : calcDgssDsliprate();
176 :
177 4034992 : for (unsigned int i = 0; i < _nss; ++i)
178 48419904 : for (unsigned int j = 0; j < _nss; ++j)
179 : {
180 44695296 : _jacob(i, j) = 0.0;
181 44695296 : if (i == j)
182 3724608 : _jacob(i, j) += 1.0;
183 44695296 : _jacob(i, j) -= _dsliprate_dsliprate(i, j);
184 : }
185 310384 : }
186 :
187 : void
188 310384 : FiniteStrainCPSlipRateRes::calcDtauDsliprate()
189 : {
190 310384 : RankFourTensor dfedfpinv, deedfe, dfpinvdpk2;
191 310384 : std::vector<RankTwoTensor> dtaudpk2(_nss), dfpinvdsliprate(_nss);
192 :
193 4034992 : for (unsigned int i = 0; i < _nss; ++i)
194 : {
195 3724608 : dtaudpk2[i] = _s0[i];
196 3724608 : dfpinvdsliprate[i] = -_fp_old_inv * _s0[i] * _dt;
197 : }
198 :
199 1241536 : for (const auto i : make_range(Moose::dim))
200 3724608 : for (const auto j : make_range(Moose::dim))
201 11173824 : for (const auto k : make_range(Moose::dim))
202 8380368 : dfedfpinv(i, j, k, j) = _dfgrd_tmp(i, k);
203 :
204 1241536 : for (const auto i : make_range(Moose::dim))
205 3724608 : for (const auto j : make_range(Moose::dim))
206 11173824 : for (const auto k : make_range(Moose::dim))
207 : {
208 8380368 : deedfe(i, j, k, i) = deedfe(i, j, k, i) + _fe(k, j) * 0.5;
209 8380368 : deedfe(i, j, k, j) = deedfe(i, j, k, j) + _fe(k, i) * 0.5;
210 : }
211 :
212 310384 : RankFourTensor dpk2dfpinv;
213 :
214 310384 : dpk2dfpinv = _elasticity_tensor[_qp] * deedfe * dfedfpinv;
215 :
216 4034992 : for (unsigned int i = 0; i < _nss; ++i)
217 48419904 : for (unsigned int j = 0; j < _nss; ++j)
218 44695296 : _dsliprate_dsliprate(i, j) =
219 44695296 : _dslipdtau(i) * dtaudpk2[i].doubleContraction(dpk2dfpinv * dfpinvdsliprate[j]);
220 310384 : }
221 :
222 : void
223 310384 : FiniteStrainCPSlipRateRes::calcDgssDsliprate()
224 : {
225 4034992 : for (unsigned int i = 0; i < _nss; ++i)
226 48419904 : for (unsigned int j = 0; j < _nss; ++j)
227 44695296 : _dsliprate_dsliprate(i, j) += _dsliprate_dgss(i) * _dgss_dsliprate(i, j);
228 310384 : }
229 :
230 : void
231 320184 : FiniteStrainCPSlipRateRes::getSlipIncrements()
232 : {
233 320184 : FiniteStrainCrystalPlasticity::getSlipIncrements();
234 :
235 320184 : if (_err_tol)
236 : return;
237 :
238 310384 : _dslipdtau *= 1.0 / _dt;
239 :
240 4034992 : for (unsigned int i = 0; i < _nss; ++i)
241 3724608 : _dsliprate_dgss(i) = -_a0(i) / _xm(i) *
242 3724608 : std::pow(std::abs(_tau(i) / _gss_tmp[i]), 1.0 / _xm(i) - 1.0) * _tau(i) /
243 3724608 : std::pow(_gss_tmp[i], 2.0);
244 : }
245 :
246 : void
247 277728 : FiniteStrainCPSlipRateRes::calcUpdate()
248 : {
249 277728 : DenseMatrix<Real> A = _jacob;
250 277728 : DenseVector<Real> r(_nss);
251 277728 : DenseVector<Real> x(_nss);
252 :
253 : r = _resid;
254 :
255 277728 : A.lu_solve(r, x);
256 :
257 : _resid = x;
258 277728 : }
259 :
260 : Real
261 449568 : FiniteStrainCPSlipRateRes::calcResidNorm()
262 : {
263 : Real rnorm = 0.0;
264 5844384 : for (unsigned int i = 0; i < _nss; ++i)
265 5394816 : rnorm += Utility::pow<2>(_resid(i));
266 449568 : rnorm = std::sqrt(rnorm) / _nss;
267 :
268 449568 : return rnorm;
269 : }
270 :
271 : bool
272 0 : FiniteStrainCPSlipRateRes::lineSearchUpdateSlipRate(const Real rnorm_prev,
273 : const DenseVector<Real> & update)
274 : {
275 0 : if (_lsrch_method == "CUT_HALF")
276 : {
277 : Real rnorm;
278 : Real step = 1.0;
279 : do
280 : {
281 0 : for (unsigned int i = 0; i < update.size(); ++i)
282 0 : _slip_rate(i) += step * update(i);
283 :
284 0 : step /= 2.0;
285 :
286 0 : for (unsigned int i = 0; i < update.size(); ++i)
287 0 : _slip_rate(i) -= step * update(i);
288 :
289 0 : calcResidualSlipRate();
290 0 : if (_err_tol)
291 : return false;
292 0 : rnorm = calcResidNorm();
293 0 : } while (rnorm > rnorm_prev && step > _min_lsrch_step);
294 :
295 0 : if (rnorm > rnorm_prev && step <= _min_lsrch_step)
296 : return false;
297 :
298 : return true;
299 : }
300 0 : else if (_lsrch_method == "BISECTION")
301 : {
302 : unsigned int count = 0;
303 : Real step_a = 0.0;
304 : Real step_b = 1.0;
305 : Real step = 1.0;
306 : Real s_m = 1000.0;
307 : Real rnorm = 1000.0;
308 :
309 0 : Real s_b = calcResidDotProdUpdate(update);
310 0 : Real rnorm1 = calcResidNorm();
311 :
312 0 : for (unsigned int i = 0; i < update.size(); ++i)
313 0 : _slip_rate(i) += update(i);
314 :
315 0 : calcResidualSlipRate();
316 0 : Real s_a = calcResidDotProdUpdate(update);
317 0 : Real rnorm0 = calcResidNorm();
318 :
319 0 : for (unsigned int i = 0; i < update.size(); ++i)
320 0 : _slip_rate(i) -= update(i);
321 :
322 0 : if ((rnorm1 / rnorm0) < _lsrch_tol || s_a * s_b > 0)
323 : {
324 0 : calcResidualSlipRate();
325 0 : return true;
326 : }
327 :
328 0 : while ((rnorm / rnorm0) > _lsrch_tol && count < _lsrch_max_iter)
329 : {
330 :
331 0 : for (unsigned int i = 0; i < update.size(); ++i)
332 0 : _slip_rate(i) += step * update(i);
333 :
334 0 : step = 0.5 * (step_a + step_b);
335 :
336 0 : for (unsigned int i = 0; i < update.size(); ++i)
337 0 : _slip_rate(i) -= step * update(i);
338 :
339 0 : calcResidualSlipRate();
340 0 : s_m = calcResidDotProdUpdate(update);
341 0 : rnorm = calcResidNorm();
342 :
343 0 : if (s_m * s_a < 0.0)
344 : {
345 : step_b = step;
346 : s_b = s_m;
347 : }
348 0 : if (s_m * s_b < 0.0)
349 : {
350 : step_a = step;
351 : s_a = s_m;
352 : }
353 0 : count++;
354 : }
355 :
356 0 : if ((rnorm / rnorm0) < _lsrch_tol && count < _lsrch_max_iter)
357 : return true;
358 :
359 : return false;
360 : }
361 : else
362 : {
363 0 : mooseError("Line search meothod is not provided.");
364 : return false;
365 : }
366 : }
367 :
368 : Real
369 0 : FiniteStrainCPSlipRateRes::calcResidDotProdUpdate(const DenseVector<Real> & update)
370 : {
371 : Real dotprod = 0.0;
372 0 : for (unsigned int i = 0; i < _nss; ++i)
373 0 : dotprod += _resid(i) * update(i);
374 0 : return dotprod;
375 : }
|
