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