Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
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("SolidMechanicsApp", FiniteStrainCPSlipRateRes, "11/15/2024 12:00");
14 :
15 : InputParameters
16 48 : FiniteStrainCPSlipRateRes::validParams()
17 : {
18 48 : InputParameters params = FiniteStrainCrystalPlasticity::validParams();
19 48 : params.addClassDescription("Deprecated class: please use CrystalPlasticityKalidindiUpdate and "
20 : "ComputeMultipleCrystalPlasticityStress instead.");
21 48 : return params;
22 0 : }
23 :
24 36 : FiniteStrainCPSlipRateRes::FiniteStrainCPSlipRateRes(const InputParameters & parameters)
25 : : FiniteStrainCrystalPlasticity(parameters),
26 36 : _resid(_nss),
27 36 : _slip_rate(_nss),
28 36 : _dsliprate_dgss(_nss),
29 36 : _jacob(_nss, _nss),
30 72 : _dsliprate_dsliprate(_nss, _nss)
31 : {
32 36 : }
33 :
34 : void
35 61560 : FiniteStrainCPSlipRateRes::solveStatevar()
36 : {
37 61560 : preSolveStress();
38 61560 : solveStress();
39 61560 : if (_err_tol)
40 : return;
41 47312 : postSolveStress();
42 : }
43 :
44 : void
45 61560 : FiniteStrainCPSlipRateRes::preSolveStress()
46 : {
47 61560 : FiniteStrainCrystalPlasticity::preSolveStress();
48 : _slip_rate.zero();
49 61560 : }
50 :
51 : void
52 61560 : 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 61560 : calcResidJacobSlipRate();
62 61560 : if (_err_tol)
63 : return;
64 47312 : 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 453088 : while (rnorm > _rtol * rnorm0 && rnorm0 > _abs_tol && iter < _maxiter)
77 : {
78 405776 : calcUpdate();
79 :
80 : DenseVector<Real> update = _resid;
81 :
82 405776 : _slip_rate -= update;
83 :
84 405776 : calcResidualSlipRate();
85 405776 : if (_err_tol)
86 : return;
87 : rnorm_prev = rnorm;
88 405776 : rnorm = calcResidNorm();
89 :
90 405776 : 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 405776 : calcJacobianSlipRate();
100 :
101 405776 : if (_use_line_search)
102 195592 : rnorm = calcResidNorm();
103 405776 : 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 47312 : 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 61560 : FiniteStrainCPSlipRateRes::calcResidJacobSlipRate()
127 : {
128 61560 : calcResidualSlipRate();
129 61560 : if (_err_tol)
130 : return;
131 47312 : calcJacobianSlipRate();
132 : }
133 :
134 : void
135 467336 : FiniteStrainCPSlipRateRes::calcResidualSlipRate()
136 : {
137 467336 : RankTwoTensor eqv_slip_incr, ce, ee;
138 467336 : const RankTwoTensor iden(RankTwoTensor::initIdentity);
139 :
140 : _slip_incr = _slip_rate;
141 467336 : _slip_incr *= _dt;
142 :
143 6075368 : for (unsigned int i = 0; i < _nss; ++i)
144 5608032 : eqv_slip_incr += _s0[i] * _slip_incr(i);
145 :
146 467336 : eqv_slip_incr = iden - eqv_slip_incr;
147 :
148 467336 : _fp_inv = _fp_old_inv * eqv_slip_incr;
149 467336 : _fe = _dfgrd_tmp * _fp_inv;
150 :
151 467336 : ce = _fe.transpose() * _fe;
152 467336 : ee = ce - iden;
153 467336 : ee *= 0.5;
154 :
155 467336 : _pk2_tmp = _elasticity_tensor[_qp] * ee;
156 :
157 6075368 : for (unsigned int i = 0; i < _nss; ++i)
158 5608032 : _tau(i) = _pk2_tmp.doubleContraction(_s0[i]);
159 :
160 467336 : update_slip_system_resistance();
161 467336 : getSlipIncrements();
162 :
163 467336 : if (_err_tol)
164 : return;
165 :
166 5890144 : for (unsigned int i = 0; i < _nss; ++i)
167 5437056 : _resid(i) = _slip_rate(i) - _slip_incr(i) / _dt;
168 : }
169 :
170 : void
171 453088 : FiniteStrainCPSlipRateRes::calcJacobianSlipRate()
172 : {
173 : //_dsliprate_dsliprate not reinitialized to zero, hence order is important
174 453088 : calcDtauDsliprate();
175 453088 : calcDgssDsliprate();
176 :
177 5890144 : for (unsigned int i = 0; i < _nss; ++i)
178 70681728 : for (unsigned int j = 0; j < _nss; ++j)
179 : {
180 65244672 : _jacob(i, j) = 0.0;
181 65244672 : if (i == j)
182 5437056 : _jacob(i, j) += 1.0;
183 65244672 : _jacob(i, j) -= _dsliprate_dsliprate(i, j);
184 : }
185 453088 : }
186 :
187 : void
188 453088 : FiniteStrainCPSlipRateRes::calcDtauDsliprate()
189 : {
190 453088 : RankFourTensor dfedfpinv, deedfe, dfpinvdpk2;
191 453088 : std::vector<RankTwoTensor> dtaudpk2(_nss), dfpinvdsliprate(_nss);
192 :
193 5890144 : for (unsigned int i = 0; i < _nss; ++i)
194 : {
195 5437056 : dtaudpk2[i] = _s0[i];
196 5437056 : dfpinvdsliprate[i] = -_fp_old_inv * _s0[i] * _dt;
197 : }
198 :
199 1812352 : for (const auto i : make_range(Moose::dim))
200 5437056 : for (const auto j : make_range(Moose::dim))
201 16311168 : for (const auto k : make_range(Moose::dim))
202 12233376 : dfedfpinv(i, j, k, j) = _dfgrd_tmp(i, k);
203 :
204 1812352 : for (const auto i : make_range(Moose::dim))
205 5437056 : for (const auto j : make_range(Moose::dim))
206 16311168 : for (const auto k : make_range(Moose::dim))
207 : {
208 12233376 : deedfe(i, j, k, i) = deedfe(i, j, k, i) + _fe(k, j) * 0.5;
209 12233376 : deedfe(i, j, k, j) = deedfe(i, j, k, j) + _fe(k, i) * 0.5;
210 : }
211 :
212 453088 : RankFourTensor dpk2dfpinv;
213 :
214 453088 : dpk2dfpinv = _elasticity_tensor[_qp] * deedfe * dfedfpinv;
215 :
216 5890144 : for (unsigned int i = 0; i < _nss; ++i)
217 70681728 : for (unsigned int j = 0; j < _nss; ++j)
218 65244672 : _dsliprate_dsliprate(i, j) =
219 65244672 : _dslipdtau(i) * dtaudpk2[i].doubleContraction(dpk2dfpinv * dfpinvdsliprate[j]);
220 453088 : }
221 :
222 : void
223 453088 : FiniteStrainCPSlipRateRes::calcDgssDsliprate()
224 : {
225 5890144 : for (unsigned int i = 0; i < _nss; ++i)
226 70681728 : for (unsigned int j = 0; j < _nss; ++j)
227 65244672 : _dsliprate_dsliprate(i, j) += _dsliprate_dgss(i) * _dgss_dsliprate(i, j);
228 453088 : }
229 :
230 : void
231 467336 : FiniteStrainCPSlipRateRes::getSlipIncrements()
232 : {
233 467336 : FiniteStrainCrystalPlasticity::getSlipIncrements();
234 :
235 467336 : if (_err_tol)
236 : return;
237 :
238 453088 : _dslipdtau *= 1.0 / _dt;
239 :
240 5890144 : for (unsigned int i = 0; i < _nss; ++i)
241 5437056 : _dsliprate_dgss(i) = -_a0(i) / _xm(i) *
242 5437056 : std::pow(std::abs(_tau(i) / _gss_tmp[i]), 1.0 / _xm(i) - 1.0) * _tau(i) /
243 5437056 : std::pow(_gss_tmp[i], 2.0);
244 : }
245 :
246 : void
247 405776 : FiniteStrainCPSlipRateRes::calcUpdate()
248 : {
249 405776 : DenseMatrix<Real> A = _jacob;
250 405776 : DenseVector<Real> r(_nss);
251 405776 : DenseVector<Real> x(_nss);
252 :
253 : r = _resid;
254 :
255 405776 : A.lu_solve(r, x);
256 :
257 : _resid = x;
258 405776 : }
259 :
260 : Real
261 648680 : FiniteStrainCPSlipRateRes::calcResidNorm()
262 : {
263 : Real rnorm = 0.0;
264 8432840 : for (unsigned int i = 0; i < _nss; ++i)
265 : rnorm += Utility::pow<2>(_resid(i));
266 648680 : rnorm = std::sqrt(rnorm) / _nss;
267 :
268 648680 : 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 : }
|
