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 "LagrangianObjectiveRate.h"
11 : #include "MooseEnum.h"
12 : #include "MooseError.h"
13 :
14 : #include <tuple>
15 :
16 : // ============================================================================
17 : // File-local helpers (formerly static members of the rate class hierarchy)
18 : // ============================================================================
19 : namespace
20 : {
21 : /// Build the J tensor that defines the linear advection `sigma_{n+1} = J(dQ)^{-1} sigma`.
22 : RankFourTensor
23 30248705 : updateTensor(const RankTwoTensor & dQ)
24 : {
25 : const auto I = RankTwoTensor::Identity();
26 : usingTensorIndices(i, j, k, l);
27 60497410 : return (1.0 + dQ.trace()) * I.times<i, k, j, l>(I) - dQ.times<i, k, j, l>(I) -
28 60497410 : I.times<i, k, j, l>(dQ);
29 : }
30 :
31 : /// Apply the linear advection to a stress, returning `(advected_stress, Jinv)`.
32 : std::tuple<RankTwoTensor, RankFourTensor>
33 1046966 : advectStress(const RankTwoTensor & S0, const RankTwoTensor & dQ)
34 : {
35 1046966 : const RankFourTensor J = updateTensor(dQ);
36 1046966 : const RankFourTensor Jinv = J.inverse();
37 1046966 : const RankTwoTensor S = Jinv * S0;
38 1046966 : return {S, Jinv};
39 : }
40 :
41 : /// Residual-only advection: solve `J(dQ) : S = S0` for the single stress RHS instead of forming
42 : /// the full fourth-order inverse. `Jinv` is only needed as an operator on the Jacobian, so on
43 : /// residual-only sweeps a one-RHS 9x9 solve replaces the 9-column inverse (~2x fewer flops on the
44 : /// most frequently evaluated path). Correct for any dQ: `J` maps the symmetric subspace to itself,
45 : /// so the solution equals `Jinv * S0` and stays symmetric even when dQ is not.
46 : RankTwoTensor
47 29201739 : advectStressSolve(const RankTwoTensor & S0, const RankTwoTensor & dQ)
48 : {
49 29201739 : const RankFourTensor J = updateTensor(dQ);
50 : // Row/col index (a, b) -> (3a + b); A_{(ij)(kl)} = J_ijkl matches RankFourTensor::operator*.
51 : Eigen::Matrix<Real, 9, 9> A;
52 : Eigen::Matrix<Real, 9, 1> b;
53 116806956 : for (const auto i : make_range(3u))
54 350420868 : for (const auto j : make_range(3u))
55 : {
56 262815651 : b(3 * i + j) = S0(i, j);
57 1051262604 : for (const auto k : make_range(3u))
58 3153787812 : for (const auto l : make_range(3u))
59 2365340859 : A(3 * i + j, 3 * k + l) = J(i, j, k, l);
60 : }
61 29201739 : const Eigen::Matrix<Real, 9, 1> x = A.partialPivLu().solve(b);
62 29201739 : RankTwoTensor S;
63 116806956 : for (const auto i : make_range(3u))
64 350420868 : for (const auto j : make_range(3u))
65 262815651 : S(i, j) = x(3 * i + j);
66 29201739 : return S;
67 : }
68 :
69 : /// Derivative of the linear advection action with respect to the kinematic tensor.
70 : RankFourTensor
71 1046966 : stressAdvectionDerivative(const RankTwoTensor & S)
72 : {
73 : const auto I = RankTwoTensor::Identity();
74 : usingTensorIndices(i, j, k, l);
75 1046966 : return S.times<i, j, k, l>(I) - I.times<i, k, l, j>(S) - S.times<i, l, j, k>(I);
76 : }
77 :
78 : /// Consistent tangent for the linear template.
79 : RankFourTensor
80 1046966 : cauchyJacobian(const RankFourTensor & Jinv,
81 : const RankFourTensor & small_jacobian,
82 : const RankFourTensor & U)
83 : {
84 1046966 : return Jinv * (small_jacobian - U);
85 : }
86 :
87 : /// Rodrigues exponential of a skew 3x3 tensor W. When `dR_dW` is non-null it is filled with the
88 : /// unconstrained chain-rule derivative; otherwise the R4 algebra is skipped. The upstream
89 : /// `d(dW)/dF` projects out non-skew perturbations.
90 : RankTwoTensor
91 3333960 : rotationFromVorticity(const RankTwoTensor & W, RankFourTensor * dR_dW)
92 : {
93 : // For a skew W in 3D, theta = sqrt((W : W) / 2). R = exp(W) via Rodrigues:
94 : // R = I + f(theta) W + g(theta) W^2, f = sin theta / theta, g = (1 - cos theta) / theta^2.
95 : // For theta -> 0, fall back to the Taylor expansion: R ~= I + W + W^2/2.
96 3333960 : const Real theta2 = 0.5 * W.doubleContraction(W);
97 3333960 : const Real theta = std::sqrt(theta2);
98 : const auto I2 = RankTwoTensor::Identity();
99 3333960 : const RankTwoTensor W2 = W * W;
100 :
101 : Real f, g, df_dth, dg_dth;
102 : const Real small_theta = 1.0e-7;
103 3333960 : if (theta < small_theta)
104 : {
105 29153 : f = 1.0 - theta2 / 6.0;
106 29153 : g = 0.5 - theta2 / 24.0;
107 29153 : df_dth = -theta / 3.0;
108 29153 : dg_dth = -theta / 12.0;
109 : }
110 : else
111 : {
112 3304807 : const Real s = std::sin(theta), c = std::cos(theta);
113 3304807 : f = s / theta;
114 3304807 : g = (1.0 - c) / theta2;
115 3304807 : df_dth = (theta * c - s) / theta2;
116 3304807 : dg_dth = (theta * s - 2.0 * (1.0 - c)) / (theta * theta2);
117 : }
118 3333960 : const RankTwoTensor R = I2 + f * W + g * W2;
119 :
120 3333960 : if (!dR_dW)
121 : return R;
122 :
123 : // dR_ij/dW_mn = (df/dtheta * dtheta/dW_mn) W_ij + f * delta_im delta_jn
124 : // + (dg/dtheta * dtheta/dW_mn) (W^2)_ij + g * (delta_im W_nj + W_im delta_jn)
125 : // dtheta/dW_mn = W_mn / (2theta) (from d(theta^2)/dW = W).
126 : usingTensorIndices(i_, j_, m_, n_);
127 : const RankFourTensor d_W2_dW =
128 79885 : I2.template times<i_, m_, n_, j_>(W) + W.template times<i_, m_, j_, n_>(I2);
129 79885 : *dR_dW = f * RankFourTensor::IdentityFour() + g * d_W2_dW;
130 79885 : if (theta >= small_theta)
131 : {
132 71937 : const Real inv_2theta = 1.0 / (2.0 * theta);
133 : // (df/dtheta * W_mn / (2theta)) * W_ij -> (df/dtheta * inv_2theta) * W (x) W (output indices
134 : // ij,mn).
135 143874 : *dR_dW += (df_dth * inv_2theta) * W.template times<i_, j_, m_, n_>(W);
136 143874 : *dR_dW += (dg_dth * inv_2theta) * W2.template times<i_, j_, m_, n_>(W);
137 : }
138 : // (For theta < small_theta the dtheta-dependent contributions are O(theta) -> 0; the
139 : // identity + W^2 parts above cover the small-angle limit cleanly.)
140 :
141 : return R;
142 : }
143 : }
144 :
145 : namespace LagrangianObjectiveRates
146 : {
147 :
148 : // ============================================================================
149 : // Truesdell
150 : // ============================================================================
151 : Outputs
152 23928368 : truesdell(const Inputs & in, bool need_jacobian)
153 : {
154 23928368 : Outputs out;
155 23928368 : if (!need_jacobian)
156 : {
157 22914807 : out.cauchy_stress = advectStressSolve(in.cauchy_stress_old + in.dS, in.dL);
158 22914807 : return out;
159 : }
160 :
161 1013561 : auto [S, Jinv] = advectStress(in.cauchy_stress_old + in.dS, in.dL);
162 1013561 : out.cauchy_stress = S;
163 1013561 : out.dcauchy_stress_d_eigenstrain = -Jinv * in.small_jacobian;
164 :
165 1013561 : const RankFourTensor U = stressAdvectionDerivative(S);
166 1013561 : out.cauchy_jacobian = cauchyJacobian(Jinv, in.small_jacobian, U);
167 : return out;
168 : }
169 :
170 : // ============================================================================
171 : // Jaumann
172 : // ============================================================================
173 : Outputs
174 3160988 : jaumann(const Inputs & in, bool need_jacobian)
175 : {
176 3160988 : Outputs out;
177 3160988 : if (!need_jacobian)
178 : {
179 3143979 : out.cauchy_stress = advectStressSolve(in.cauchy_stress_old + in.dS, in.dW);
180 3143979 : return out;
181 : }
182 :
183 17009 : auto [S, Jinv] = advectStress(in.cauchy_stress_old + in.dS, in.dW);
184 17009 : out.cauchy_stress = S;
185 17009 : out.dcauchy_stress_d_eigenstrain = -Jinv * in.small_jacobian;
186 :
187 17009 : const RankFourTensor d_dW_d_dL = in.d_dW_d_F * in.d_dL_d_F.inverse();
188 17009 : const RankFourTensor U = stressAdvectionDerivative(S) * d_dW_d_dL;
189 17009 : out.cauchy_jacobian = cauchyJacobian(Jinv, in.small_jacobian, U);
190 : return out;
191 : }
192 :
193 : // ============================================================================
194 : // Green-Naghdi
195 : // ============================================================================
196 : Outputs
197 3159349 : greenNaghdi(const Inputs & in, bool need_jacobian)
198 : {
199 : usingTensorIndices(i, j, k, l, m);
200 :
201 : const RankTwoTensor I = RankTwoTensor::Identity();
202 3159349 : const RankTwoTensor dR = in.rotation * in.rotation_old.transpose() - I;
203 3159349 : const RankTwoTensor dO = dR * in.inv_df;
204 :
205 3159349 : Outputs out;
206 3159349 : if (!need_jacobian)
207 : {
208 3142953 : out.cauchy_stress = advectStressSolve(in.cauchy_stress_old + in.dS, dO);
209 3142953 : return out;
210 : }
211 :
212 16396 : auto [S, Jinv] = advectStress(in.cauchy_stress_old + in.dS, dO);
213 16396 : out.cauchy_stress = S;
214 16396 : out.dcauchy_stress_d_eigenstrain = -Jinv * in.small_jacobian;
215 :
216 16396 : const RankFourTensor & d_R_d_F = in.d_rotation_d_F;
217 16396 : const RankFourTensor d_F_d_dL = in.d_dL_d_F.inverse();
218 16396 : const RankTwoTensor T = in.rotation_old.transpose() * in.inv_df;
219 :
220 16396 : const RankFourTensor d_invdf_d_F = -in.inv_df.times<i, k, l, j>(in.inv_def_grad);
221 16396 : const RankFourTensor d_invdf_d_dL = d_invdf_d_F * d_F_d_dL;
222 16396 : const RankFourTensor d_dO_d_invdf = dR.times<i, k, j, l>(I);
223 :
224 : const RankFourTensor d_dO_d_dL =
225 16396 : T.times<m, j, i, m, k, l>(d_R_d_F * d_F_d_dL) + d_dO_d_invdf * d_invdf_d_dL;
226 16396 : const RankFourTensor U = stressAdvectionDerivative(S) * d_dO_d_dL;
227 16396 : out.cauchy_jacobian = cauchyJacobian(Jinv, in.small_jacobian, U);
228 : return out;
229 : }
230 :
231 : // ============================================================================
232 : // Rashid
233 : // ============================================================================
234 : Outputs
235 3333960 : rashid(const Inputs & in, bool need_jacobian)
236 : {
237 : usingTensorIndices(i_, j_, k_, l_, m_, n_);
238 :
239 : // r_hat = exp(Deltaw). Skip its R4 derivative when we don't need the Jacobian.
240 3333960 : RankFourTensor d_rhat_d_dW;
241 6588035 : const RankTwoTensor rhat = rotationFromVorticity(in.dW, need_jacobian ? &d_rhat_d_dW : nullptr);
242 3333960 : const RankTwoTensor rhatT = rhat.transpose();
243 :
244 : // sigma_{n+1} = r_hat (sigma_n + Deltasigma) r_hat^T (eq. 22)
245 3333960 : const RankTwoTensor S_inner = in.cauchy_stress_old + in.dS;
246 3333960 : Outputs out;
247 3333960 : out.cauchy_stress = rhat * S_inner * rhatT;
248 :
249 3333960 : if (!need_jacobian)
250 : return out;
251 :
252 : // J^{-1}_{ijkl} = r_hat_ik r_hat_jl (eq. 24), the rank-4 (r_hat (x) r_hat) sandwich. Used for
253 : // the eigenstrain Jacobian and as the "outer" operator for the constitutive piece. We name the
254 : // dummy "middle pair" k, l here (output indices i, j, k, l) -- when chained via
255 : // RankFourTensor::operator*, the k, l contract against small_jacobian's first pair, as desired.
256 79885 : const RankFourTensor J_inv = rhat.times<i_, k_, j_, l_>(rhat);
257 :
258 : // dsigma/d(eigenstrain) = -J^{-1} : small_jacobian.
259 : // (mechanical_strain = total_strain - eigenstrain, hence the minus.)
260 79885 : out.dcauchy_stress_d_eigenstrain = -(J_inv * in.small_jacobian);
261 :
262 : // Chain rule pieces for dsigma/d(dL):
263 : // d(Deltaw)/d(dL) = d(dW)/dF * inverse(d(dL)/dF)
264 : // d(Deltad)/d(dL) = I^(4) - d(Deltaw)/d(dL) (Deltad + Deltaw = dL by construction).
265 : // d(r_hat)/d(dL) = d(r_hat)/d(Deltaw) * d(Deltaw)/d(dL).
266 79885 : const RankFourTensor d_F_d_dL = in.d_dL_d_F.inverse();
267 79885 : const RankFourTensor d_dW_d_dL = in.d_dW_d_F * d_F_d_dL;
268 79885 : const RankFourTensor d_dD_d_dL = RankFourTensor::IdentityFour() - d_dW_d_dL;
269 79885 : const RankFourTensor d_rhat_d_dL = d_rhat_d_dW * d_dW_d_dL;
270 :
271 : // Rotation pieces of dsigma/d(dL):
272 : // T1_{ijkl} = (d_rhat_d_dL)_{imkl} * (S r_hat^T)_{mj}
273 : // T2_{ijkl} = (r_hat S)_{im} * (d_rhat_d_dL)_{jmkl}
274 79885 : const RankTwoTensor SR = S_inner * rhatT; // S r_hat^T (shape (m, j))
275 79885 : const RankTwoTensor RS = rhat * S_inner; // r_hat S (shape (i, m))
276 79885 : const RankFourTensor T1 = SR.times<m_, j_, i_, m_, k_, l_>(d_rhat_d_dL);
277 79885 : const RankFourTensor T2 = RS.times<i_, m_, j_, m_, k_, l_>(d_rhat_d_dL);
278 :
279 : // Constitutive piece: J^{-1} : (small_jacobian * d(Deltad)/d(dL)).
280 79885 : const RankFourTensor T3 = J_inv * (in.small_jacobian * d_dD_d_dL);
281 :
282 79885 : out.cauchy_jacobian = T1 + T2 + T3;
283 : return out;
284 : }
285 :
286 : // ============================================================================
287 : // Dispatch
288 : // ============================================================================
289 : Outputs
290 33582604 : compute(const MooseEnum & rate, const Inputs & in, bool need_jacobian)
291 : {
292 : const std::string s = rate;
293 33582604 : if (s == "truesdell")
294 23928348 : return truesdell(in, need_jacobian);
295 9654256 : if (s == "jaumann")
296 3160968 : return jaumann(in, need_jacobian);
297 6493288 : if (s == "green_naghdi")
298 3159348 : return greenNaghdi(in, need_jacobian);
299 3333940 : if (s == "rashid")
300 3333940 : return rashid(in, need_jacobian);
301 0 : mooseError("Unknown objective_rate value: ", s);
302 : }
303 : }
|
