www.mooseframework.org
RateDepSmearCrackModel.C
Go to the documentation of this file.
1 /****************************************************************/
2 /* MOOSE - Multiphysics Object Oriented Simulation Environment */
3 /* */
4 /* All contents are licensed under LGPL V2.1 */
5 /* See LICENSE for full restrictions */
6 /****************************************************************/
8 
9 template <>
10 InputParameters
12 {
13 
14  InputParameters params = validParams<ConstitutiveModel>();
15 
16  params.addRequiredParam<Real>("ref_damage_rate", "Reference damage rate");
17  params.addRequiredParam<unsigned int>("nstate", "Number of state variables");
18  params.addParam<Real>("exponent", 1.0, "Power law exponent");
19  params.addParam<unsigned int>("maxiter", 20, "Constitutive update iteration");
20  params.addParam<Real>("tol", 1e-5, "Constitutive update tolerance");
21  params.addParam<Real>("zero_tol", 1e-8, "Tolerance for zero");
22  params.addParam<Real>(
23  "intvar_incr_tol", 0.1, "Allowable relative increment size for state variables");
24  params.addParam<bool>("input_random_scaling_var",
25  false,
26  "Flag to specify scaling parameter to generate random stress");
27  params.addParam<Real>("random_scaling_var",
28  1e9,
29  "Scaling value: Too large a value can cause "
30  "non-positive definiteness - use 0.1 of young's "
31  "modulus");
32 
33  return params;
34 }
35 
36 RateDepSmearCrackModel::RateDepSmearCrackModel(const InputParameters & parameters)
37  : ConstitutiveModel(parameters),
38  _ref_damage_rate(getParam<Real>("ref_damage_rate")),
39  _nstate(getParam<unsigned int>("nstate")),
40  _exponent(getParam<Real>("exponent")),
41  _maxiter(getParam<unsigned int>("maxiter")),
42  _tol(getParam<Real>("tol")),
43  _zero_tol(getParam<Real>("zero_tol")),
44  _intvar_incr_tol(getParam<Real>("intvar_incr_tol")),
45  _input_rndm_scale_var(getParam<bool>("input_random_scaling_var")),
46  _rndm_scale_var(getParam<Real>("random_scaling_var")),
47  _intvar(declareProperty<std::vector<Real>>("intvar")),
48  _intvar_old(getMaterialPropertyOld<std::vector<Real>>("intvar")),
49  _stress_undamaged(declareProperty<SymmTensor>("stress_undamaged")),
50  _stress_undamaged_old(getMaterialPropertyOld<SymmTensor>("stress_undamaged"))
51 {
52 
53  _intvar_incr.resize(_nstate, 0.0);
54  _intvar_tmp.resize(_nstate, 0.0);
55  _intvar_old_tmp.resize(_nstate, 0.0);
56  _resid.resize(_nstate, 0.0);
57  _jac.resize(_nstate * _nstate, 0.0);
58  _dvar.resize(_nstate, 0.0);
59 
61 }
62 
63 void
65 {
66  _intvar[_qp].resize(_nstate, 0.0);
67  const_cast<MaterialProperty<std::vector<Real>> &>(_intvar_old)[_qp].resize(_nstate, 0.0);
68 }
69 
70 void
71 RateDepSmearCrackModel::computeStress(const Elem & /*current_elem*/,
72  const SymmElasticityTensor & elasticityTensor,
73  const SymmTensor & stress_old,
74  SymmTensor & strain_increment,
75  SymmTensor & stress_new)
76 {
77  _elasticity = elasticityTensor;
78  _stress_old = stress_old;
79  _dstrain = strain_increment;
80 
81  initVariables();
82  solve();
83 
85  _rndm_scale_var = elasticityTensor.valueAtIndex(0);
86 
87  if (_nconv || _err_tol)
88  {
89  mooseWarning("RateDepSmearCrackModel: Constitutive cutback");
91  }
92  else
93  {
96 
97  stress_new = _stress_new;
99  }
100 }
101 
102 void
104 {
105 
108 
109  for (unsigned int i = 0; i < _nstate; ++i)
110  {
111  _intvar_tmp[i] = _intvar_old[_qp][i];
112  _intvar_old_tmp[i] = _intvar_old[_qp][i];
113  }
114 }
115 
116 void
118 {
119 
120  unsigned int iter = 0;
121  _err_tol = false;
122 
123  calcResidual();
124 
126 
127  while (_nconv && iter < _maxiter && !_err_tol)
128  {
129  calcJacobian();
130  updateVariables();
131  calcResidual();
132 
134 
135  iter++;
136  }
137 }
138 
139 void
141 {
142  int error = matrixInversion(_jac, _nstate);
143  if (error != 0)
144  mooseError("Error in Matrix Inversion in RankFourTensor");
145 
146  for (unsigned int i = 0; i < _nstate; i++)
147  {
148  _dvar[i] = 0.0;
149  for (unsigned int j = 0; j < _nstate; j++)
150  _dvar[i] += _jac[i * _nstate + j] * _resid[j];
151  }
152 
153  for (unsigned int i = 0; i < _nstate; i++)
154  _intvar_tmp[i] -= _dvar[i];
155 }
156 
157 bool
159 {
160  Real vold, r;
161 
162  for (unsigned int i = 0; i < _nstate; i++)
163  {
164  vold = std::abs(_intvar_old_tmp[i]);
165  r = std::abs(_resid[i]);
166 
167  if (vold > _zero_tol)
168  {
169  if (r > _tol * vold)
170  return true;
171  }
172  else
173  {
174  if (r > _zero_tol)
175  return true;
176  }
177  }
178  return false;
179 }
180 
181 void
183 {
184  for (unsigned int i = 0; i < _nstate; ++i)
185  _intvar[_qp][i] = _intvar_tmp[i];
186 }
187 
188 void
190 {
191 }
192 
193 void
195 {
196  calcStateIncr();
197 
198  if (_err_tol)
199  return;
200 
201  for (unsigned int i = 0; i < _nstate; ++i)
203 }
204 
205 void
207 {
208 }
209 
210 void
212 {
213 }
214 
215 int
216 RateDepSmearCrackModel::matrixInversion(std::vector<Real> & A, int n) const
217 {
218  int return_value, buffer_size = n * 64;
219  std::vector<PetscBLASInt> ipiv(n);
220  std::vector<PetscScalar> buffer(buffer_size);
221 
222  // Following does a LU decomposition of "square matrix A"
223  // upon return "A = P*L*U" if return_value == 0
224  // Here i use quotes because A is actually an array of length n^2, not a matrix of size n-by-n
225  LAPACKgetrf_(&n, &n, &A[0], &n, &ipiv[0], &return_value);
226 
227  if (return_value != 0)
228  // couldn't LU decompose because: illegal value in A; or, A singular
229  return return_value;
230 
231 // get the inverse of A
232 #if PETSC_VERSION_LESS_THAN(3, 5, 0)
233  FORTRAN_CALL(dgetri)(&n, &A[0], &n, &ipiv[0], &buffer[0], &buffer_size, &return_value);
234 #else
235  LAPACKgetri_(&n, &A[0], &n, &ipiv[0], &buffer[0], &buffer_size, &return_value);
236 #endif
237 
238  return return_value;
239 }
240 
This class defines a basic set of capabilities any elasticity tensor should have. ...
void FORTRAN_CALL() dgetri(...)
static void initRandom()
Definition: SymmTensor.h:440
virtual void initQpStatefulProperties()
Real _tol
Maximum number of Newton Raphson iteration.
Real _rndm_scale_var
Flag to specify scaling parameter to generate random stress.
virtual void computeStress(const Elem &current_elem, const SymmElasticityTensor &elasticity_tensor, const SymmTensor &stress_old, SymmTensor &strain_increment, SymmTensor &stress_new)
std::vector< Real > _intvar_tmp
std::vector< Real > _intvar_incr
RateDepSmearCrackModel(const InputParameters &parameters)
Real valueAtIndex(int i) const
int matrixInversion(std::vector< Real > &A, int n) const
InputParameters validParams< RateDepSmearCrackModel >()
static SymmTensor genRandomSymmTensor(Real scalefactor)
Definition: SymmTensor.h:447
virtual void solve()
This function solves the state variables.
virtual void updateVariables()
This function updates variables during solve a(i+1) = a(i) + da(i+1)
const MaterialProperty< SymmTensor > & _stress_undamaged_old
bool _err_tol
Convergence flag.
const MaterialProperty< std::vector< Real > > & _intvar_old
MaterialProperty< SymmTensor > & _stress_undamaged
virtual void calcResidual()
This function calculates the residual as r = v - v_old - dv.
virtual void postSolveStress()
This function updates the stress after solve.
unsigned int _nstate
reference damage rate
InputParameters validParams< ConstitutiveModel >()
virtual void calcStateIncr()
This function calculated thw increment of the state variables (dv) used to form the residual...
bool _input_rndm_scale_var
Allowable relative increment size of state variables (dv)
virtual void postSolveVariables()
This function updates the internal variables after solve.
MaterialProperty< std::vector< Real > > & _intvar
Variable value.
virtual void calcJacobian()
This function calculates the Jacobian.
std::vector< Real > _intvar_old_tmp
SymmElasticityTensor _elasticity
Real _zero_tol
Relative tolerance factor for convergence of the Newton Raphson solve.
bool getConvergeVar()
This function returns true if convergence is not achieved.