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