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 : #pragma once
11 :
12 : #include "MultiPlasticityRawComponentAssembler.h"
13 :
14 : /**
15 : * MultiPlasticityLinearSystem computes the linear system
16 : * and handles linear-dependence removal
17 : * for use in FiniteStrainMultiPlasticity
18 : *
19 : * Note that if run in debug mode you might have to use
20 : * the --no-trap-fpe flag because PETSc-LAPACK-BLAS
21 : * explicitly compute 0/0 and 1/0, and this causes
22 : * Libmesh to trap the floating-point exceptions
23 : *
24 : * These routines are quite complicated, so here is an
25 : * extended explanation.
26 : *
27 : *
28 : * SURFACES AND MODELS
29 : *
30 : * Each plasticity model can have multiple surfaces
31 : * (eg, Mohr-Coulomb has 6 surfaces), and one internal
32 : * parameter. This is also described in
33 : * MultiPlasticityRawComponentAssembler. The
34 : *
35 : *
36 : * VARIABLE NAMES
37 : *
38 : * _num_surfaces = total number of surfaces
39 : * _num_models = total number of plasticity models
40 : * pm = plasticity multiplier. pm.size() = _num_surfaces
41 : * intnl = internal variable. intnl.size() = _num_models
42 : *
43 : *
44 : * DEGREES OF FREEDOM
45 : *
46 : * The degrees of freedom are:
47 : * the 6 components of stress (it is assumed to be symmetric)
48 : * the plasticity multipliers, pm.
49 : * the internal parameters, intnl.
50 : *
51 : * Note that in any single Newton-Raphson (NR) iteration, the number
52 : * of pm and intnl may be different from any other NR iteration.
53 : * This is because of deactivating surfaces because they
54 : * are deemed unimportant by the calling program (eg, their
55 : * yield function is < 0), or because their flow direction is
56 : * linearly dependent on other surfaces. Therefore:
57 : *
58 : * - there are "not active" surfaces, and these *never enter*
59 : * any right-hand-side (RHS), or jacobian, residual, etc calculations.
60 : * It is as if they never existed.
61 : *
62 : * - there are "active" surfaces, and *all these* enter RHS,
63 : * jacobian, etc, calculations. This set of active surfaces
64 : * may further decompose into "linearly independent" and "linearly
65 : * dependent" surfaces. The latter are surfaces whose flow direction
66 : * is linearly dependent on the former. Only the dofs from
67 : * "linearly independent" constraints are changed during the NR step.
68 : *
69 : * Hence, more exactly, the degrees of freedom, whose changes will be
70 : * provided in the NR step are:
71 : * the 6 components of stress (it is assumed to be symmetric)
72 : * the plasticity multipliers, pm, belonging to linearly independent surfaces
73 : * the internal parameters, intnl, belonging to models with at least one linearly-independent
74 : * surface
75 : *
76 : *
77 : * THE CONSTRAINTS AND RHS
78 : *
79 : * The variables calculated by calculateConstraints and
80 : * calculateRHS are:
81 : * epp = pm*r - E_inv*(trial_stress - stress) = pm*r - delta_dp
82 : * f = yield function [all the active constraints, including the deactivated_due_to_ld. The
83 : * latter ones will not be put into the linear system]
84 : * ic = intnl - intnl_old + pm*h [only for models that contain active surfaces]
85 : *
86 : * Here pm*r = sum_{active_alpha} pm[alpha]*r[alpha]. Note that this contains all the "active"
87 : * surfaces,
88 : * even the ones that have been deactivated_due_to_ld. in calculateConstraints, etc, r
89 : * is a
90 : * std::vector containing only all the active flow directions (including
91 : * deactivated_due_to_ld, but not the "not active").
92 : * f = all the "active" surfaces, even the ones that have been deactivated_due_to_ld. However, the
93 : * latter are not put into the RHS
94 : * pm*h = sum_{active_alpha} pm[alpha]*h[alpha]. Note that this only contains the "active"
95 : * hardening potentials, even the ones that have been deactivated_due_to_ld.
96 : * In calculateConstraints, calculateRHS and calculateJacobian,
97 : * h is a std::vector containing only these "active" ones.
98 : * Hence, the sum_{active_alpha} contains deactivated_due_to_ld contributions.
99 : * HOWEVER, if all the surfaces belonging to a model are either
100 : * "not active" or deactivated_due_to_ld, then this ic is not included in the RHS
101 : *
102 : * The RHS is
103 : * rhs = -(epp(0,0), epp(1,0), epp(1,1), epp(2,0), epp(2,1), epp(2,2), f[0], f[1], ...,
104 : * f[_num_active_f], ic[0], ic[1], ..., ic[num_active_ic])
105 : * Notice the appearance of only the i>=j "epp" components.
106 : *
107 : *
108 : * THE JACOBIAN
109 : *
110 : * This is d(-rhs)/d(dof). Remember that the dofs are dependent
111 : * on what is deactivated_due_to_ld, as specified above.
112 : * In matrix form, the Jacobian is:
113 : * ( depp_dstress depp_dpm depp_dintnl )
114 : * ( df_dstress 0 df_dintnl )
115 : * ( dic_dstress dic_dpm dic_dintnl )
116 : * For the "epp" terms, only the i>=j components are kept in the RHS, so only these terms are kept
117 : * here too
118 : */
119 0 : class MultiPlasticityLinearSystem : public MultiPlasticityRawComponentAssembler
120 : {
121 : public:
122 : static InputParameters validParams();
123 :
124 : MultiPlasticityLinearSystem(const MooseObject * moose_object);
125 :
126 : protected:
127 : /// Tolerance on the minimum ratio of singular values before flow-directions are deemed linearly dependent
128 : Real _svd_tol;
129 :
130 : /// Minimum value of the _f_tol parameters for the Yield Function User Objects
131 : Real _min_f_tol;
132 :
133 : /**
134 : * The constraints. These are set to zero (or <=0 in the case of the yield functions)
135 : * by the Newton-Raphson process, except in the case of linear-dependence which complicates
136 : * things.
137 : * @param stress The stress
138 : * @param intnl_old old values of the internal parameters
139 : * @param intnl internal parameters
140 : * @param pm Current value(s) of the plasticity multiplier(s) (consistency parameters)
141 : * @param delta_dp Change in plastic strain incurred so far during the return
142 : * @param[out] f Active yield function(s)
143 : * @param[out] r Active flow directions
144 : * @param[out] epp Plastic-strain increment constraint
145 : * @param[out] ic Active internal-parameter constraint
146 : * @param active The active constraints.
147 : */
148 : virtual void calculateConstraints(const RankTwoTensor & stress,
149 : const std::vector<Real> & intnl_old,
150 : const std::vector<Real> & intnl,
151 : const std::vector<Real> & pm,
152 : const RankTwoTensor & delta_dp,
153 : std::vector<Real> & f,
154 : std::vector<RankTwoTensor> & r,
155 : RankTwoTensor & epp,
156 : std::vector<Real> & ic,
157 : const std::vector<bool> & active);
158 :
159 : /**
160 : * Calculate the RHS which is
161 : * rhs = -(epp(0,0), epp(1,0), epp(1,1), epp(2,0), epp(2,1), epp(2,2), f[0], f[1], ..., f[num_f],
162 : * ic[0], ic[1], ..., ic[num_ic])
163 : *
164 : * Note that the 'epp' components only contain the upper diagonal. These contain flow directions
165 : * and plasticity-multipliers for all active surfaces, even the deactivated_due_to_ld surfaces.
166 : * Note that the 'f' components only contain the active and not deactivated_due_to_ld surfaces
167 : * Note that the 'ic' components only contain the internal constraints for models which contain
168 : * active and not deactivated_due_to_ld surfaces. They contain hardening-potentials and
169 : * plasticity-multipliers for the active surfaces, even the deactivated_due_to_ld surfaces
170 : *
171 : * @param stress The stress
172 : * @param intnl_old old values of the internal parameters
173 : * @param intnl internal parameters
174 : * @param pm Current value(s) of the plasticity multiplier(s) (consistency parameters)
175 : * @param delta_dp Change in plastic strain incurred so far during the return
176 : * @param[out] rhs the rhs
177 : * @param active The active constraints.
178 : * @param eliminate_ld Check for linear dependence of constraints and put the results into
179 : * deactivated_due_to_ld. Usually this should be true, but for certain debug operations it should
180 : * be false
181 : * @param[out] deactivated_due_to_ld constraints deactivated due to linear-dependence of flow
182 : * directions
183 : */
184 : virtual void calculateRHS(const RankTwoTensor & stress,
185 : const std::vector<Real> & intnl_old,
186 : const std::vector<Real> & intnl,
187 : const std::vector<Real> & pm,
188 : const RankTwoTensor & delta_dp,
189 : std::vector<Real> & rhs,
190 : const std::vector<bool> & active,
191 : bool eliminate_ld,
192 : std::vector<bool> & deactivated_due_to_ld);
193 :
194 : /**
195 : * d(rhs)/d(dof)
196 : */
197 : virtual void calculateJacobian(const RankTwoTensor & stress,
198 : const std::vector<Real> & intnl,
199 : const std::vector<Real> & pm,
200 : const RankFourTensor & E_inv,
201 : const std::vector<bool> & active,
202 : const std::vector<bool> & deactivated_due_to_ld,
203 : std::vector<std::vector<Real>> & jac);
204 :
205 : /**
206 : * Performs one Newton-Raphson step. The purpose here is to find the
207 : * changes, dstress, dpm and dintnl according to the Newton-Raphson procedure
208 : * @param stress Current value of stress
209 : * @param intnl_old The internal variables at the previous "time" step
210 : * @param intnl Current value of the internal variables
211 : * @param pm Current value of the plasticity multipliers (consistency parameters)
212 : * @param E_inv inverse of the elasticity tensor
213 : * @param delta_dp Current value of the plastic-strain increment (ie plastic_strain -
214 : * plastic_strain_old)
215 : * @param[out] dstress The change in stress for a full Newton step
216 : * @param[out] dpm The change in all plasticity multipliers for a full Newton step
217 : * @param[out] dintnl The change in all internal variables for a full Newton step
218 : * @param active The active constraints
219 : * @param[out] deactivated_due_to_ld The constraints deactivated due to linear-dependence of the
220 : * flow directions
221 : */
222 : virtual void nrStep(const RankTwoTensor & stress,
223 : const std::vector<Real> & intnl_old,
224 : const std::vector<Real> & intnl,
225 : const std::vector<Real> & pm,
226 : const RankFourTensor & E_inv,
227 : const RankTwoTensor & delta_dp,
228 : RankTwoTensor & dstress,
229 : std::vector<Real> & dpm,
230 : std::vector<Real> & dintnl,
231 : const std::vector<bool> & active,
232 : std::vector<bool> & deactivated_due_to_ld);
233 :
234 : private:
235 : /**
236 : * Performs a singular-value decomposition of r and returns the singular values
237 : *
238 : * Example: If r has size 5 then the singular values of the following matrix are returned:
239 : * ( r[0](0,0) r[0](0,1) r[0](0,2) r[0](1,1) r[0](1,2) r[0](2,2) )
240 : * ( r[1](0,0) r[1](0,1) r[1](0,2) r[1](1,1) r[1](1,2) r[1](2,2) )
241 : * a = ( r[2](0,0) r[2](0,1) r[2](0,2) r[2](1,1) r[2](1,2) r[2](2,2) )
242 : * ( r[3](0,0) r[3](0,1) r[3](0,2) r[3](1,1) r[3](1,2) r[3](2,2) )
243 : * ( r[4](0,0) r[4](0,1) r[4](0,2) r[4](1,1) r[4](1,2) r[4](2,2) )
244 : *
245 : * @param r The flow directions
246 : * @param[out] s The singular values
247 : * @return The return value from the PETSc LAPACK gesvd reoutine
248 : */
249 : virtual int singularValuesOfR(const std::vector<RankTwoTensor> & r, std::vector<Real> & s);
250 :
251 : /**
252 : * Performs a number of singular-value decompositions
253 : * to check for linear-dependence of the active directions "r"
254 : * If linear dependence is found, then deactivated_due_to_ld will contain 'true' entries where
255 : * surfaces need to be deactivated_due_to_ld
256 : * @param stress the current stress
257 : * @param intnl the current values of internal parameters
258 : * @param f Active yield function values
259 : * @param r the flow directions that for those yield functions that are active upon entry to this
260 : * function
261 : * @param active true if active
262 : * @param[out] deactivated_due_to_ld Yield functions deactivated due to linearly-dependent flow
263 : * directions
264 : */
265 : virtual void eliminateLinearDependence(const RankTwoTensor & stress,
266 : const std::vector<Real> & intnl,
267 : const std::vector<Real> & f,
268 : const std::vector<RankTwoTensor> & r,
269 : const std::vector<bool> & active,
270 : std::vector<bool> & deactivated_due_to_ld);
271 : };
|
