Line data    Source code

       1             : /****************************************************************/
       2             : /*               DO NOT MODIFY THIS HEADER                      */
       3             : /*                       BlackBear                              */
       4             : /*                                                              */
       5             : /*           (c) 2017 Battelle Energy Alliance, LLC             */
       6             : /*                   ALL RIGHTS RESERVED                        */
       7             : /*                                                              */
       8             : /*          Prepared by Battelle Energy Alliance, LLC           */
       9             : /*            Under Contract No. DE-AC07-05ID14517              */
      10             : /*            With the U. S. Department of Energy               */
      11             : /*                                                              */
      12             : /*            See COPYRIGHT for full restrictions               */
      13             : /****************************************************************/
      14             : 
      15             : #ifdef NEML_ENABLED
      16             : 
      17             : #include "NEMLStressBase.h"
      18             : #include "Conversion.h"
      19             : 
      20             : #include <limits>
      21             : #include <type_traits>
      22             : 
      23             : InputParameters
      24        1686 : NEMLStressBase::validParams()
      25             : {
      26        1686 :   InputParameters params = ComputeStressBase::validParams();
      27        3372 :   params.addCoupledVar("temperature", 0.0, "Coupled temperature");
      28        3372 :   params.setDocUnit("temperature", "temperature");
      29        3372 :   params.addParam<Real>("target_increment",
      30             :                         "L2 norm of the inelastic strain increment to target by adjusting the "
      31             :                         "timestep");
      32        3372 :   params.setDocUnit("target_increment", "unitless");
      33        3372 :   params.addParam<bool>("debug",
      34        3372 :                         false,
      35             :                         "Print history and strain state at the current quadrature point when a "
      36             :                         "NEML stress update fails.");
      37        1686 :   return params;
      38           0 : }
      39             : 
      40        1290 : NEMLStressBase::NEMLStressBase(const InputParameters & parameters)
      41             :   : ComputeStressBase(parameters),
      42        1290 :     _hist(declareProperty<std::vector<Real>>(_base_name + "hist")),
      43        2580 :     _hist_old(getMaterialPropertyOld<std::vector<Real>>(_base_name + "hist")),
      44        1290 :     _mechanical_strain_old(
      45        1290 :         getMaterialPropertyOldByName<RankTwoTensor>(_base_name + "mechanical_strain")),
      46        2580 :     _stress_old(getMaterialPropertyOld<RankTwoTensor>(_base_name + "stress")),
      47        1290 :     _energy(declareProperty<Real>(_base_name + "energy")),
      48        2580 :     _energy_old(getMaterialPropertyOld<Real>(_base_name + "energy")),
      49        1290 :     _dissipation(declareProperty<Real>(_base_name + "dissipation")),
      50        2580 :     _dissipation_old(getMaterialPropertyOld<Real>(_base_name + "dissipation")),
      51        1290 :     _temperature(coupledValue("temperature")),
      52        1290 :     _temperature_old(coupledValueOld("temperature")),
      53        1290 :     _inelastic_strain(declareProperty<RankTwoTensor>(_base_name + "inelastic_strain")),
      54        2580 :     _compute_dt(isParamValid("target_increment")),
      55        1332 :     _target_increment(_compute_dt ? getParam<Real>("target_increment") : 0.0),
      56        1290 :     _inelastic_strain_old(
      57        1332 :         _compute_dt ? &getMaterialPropertyOld<RankTwoTensor>(_base_name + "inelastic_strain")
      58             :                     : nullptr),
      59        1290 :     _material_dt(_compute_dt ? &declareProperty<Real>("material_timestep_limit") : nullptr),
      60        1290 :     _damage_index(nullptr),
      61        3870 :     _debug(getParam<bool>("debug"))
      62             : {
      63             :   // We're letting NEML write to raw pointers. Best make sure the stored types are
      64             :   // the same on both ends.
      65             :   static_assert(std::is_same<Real, double>::value,
      66             :                 "MOOSE/libMesh must be compiled with double precision Real types");
      67        1290 : }
      68             : 
      69             : void
      70      209924 : NEMLStressBase::computeQpStress()
      71             : {
      72             :   // We must update:
      73             :   // 1) _stress
      74             :   // 2) _Jacobian_mult
      75             :   // 3) _elastic_strain
      76             :   // 4) _history
      77             : 
      78             :   // First do some declaration and translation
      79             :   double s_np1[6];
      80             :   double s_n[6];
      81      209924 :   rankTwoTensorToNeml(_stress_old[_qp], s_n);
      82             : 
      83             :   double e_np1[6];
      84      209924 :   rankTwoTensorToNeml(_mechanical_strain[_qp], e_np1);
      85             :   double e_n[6];
      86      209924 :   rankTwoTensorToNeml(_mechanical_strain_old[_qp], e_n);
      87             : 
      88      209924 :   double t_np1 = _t;
      89      209924 :   double t_n = _t - _dt;
      90             : 
      91      209924 :   double T_np1 = _temperature[_qp];
      92      209924 :   double T_n = _temperature_old[_qp];
      93             : 
      94             :   mooseAssert(_model->nstore() == _hist[_qp].size(), "History data storage size mismatch");
      95      209924 :   double * const h_np1 = (_model->nstore() > 0 ? _hist[_qp].data() : nullptr);
      96             :   mooseAssert(_model->nstore() == _hist_old[_qp].size(), "History data storage size mismatch");
      97      209924 :   const double * const h_n = (_model->nstore() > 0 ? _hist_old[_qp].data() : nullptr);
      98             : 
      99             :   double A_np1[36];
     100             : 
     101             :   double u_np1;
     102      209924 :   double u_n = _energy_old[_qp];
     103             : 
     104             :   double p_np1;
     105      209924 :   double p_n = _dissipation_old[_qp];
     106             : 
     107             :   double estrain[6];
     108             : 
     109             :   // Actually call the update
     110             :   try
     111             :   {
     112      209924 :     _model->update_sd(
     113             :         e_np1, e_n, T_np1, T_n, t_np1, t_n, s_np1, s_n, h_np1, h_n, A_np1, u_np1, u_n, p_np1, p_n);
     114             :   }
     115           0 :   catch (const neml::NEMLError & e)
     116             :   {
     117           0 :     if (_debug)
     118           0 :       mooseException("NEML stress update failed!\n",
     119             :                      "NEML message: ",
     120             :                      e.message(),
     121             :                      "\n",
     122             :                      "_qp=",
     123             :                      _qp,
     124             :                      " _q_point=",
     125             :                      _q_point[_qp],
     126             :                      " _current_elem->id()=",
     127             :                      _current_elem->id(),
     128             :                      '\n',
     129             : 
     130             :                      "Time increment: ",
     131             :                      _dt,
     132             :                      '\n',
     133             : 
     134             :                      "Old temperature: ",
     135             :                      _temperature_old[_qp],
     136             :                      '\n',
     137             :                      "New temperature: ",
     138             :                      _temperature[_qp],
     139             :                      '\n',
     140             : 
     141             :                      "Old echanical strain: ",
     142             :                      _mechanical_strain_old[_qp],
     143             :                      '\n',
     144             :                      "New mechanical strain: ",
     145             :                      _mechanical_strain[_qp],
     146             :                      '\n',
     147             : 
     148             :                      "Old stress: ",
     149             :                      _stress_old[_qp],
     150             :                      '\n',
     151             : 
     152             :                      "Old history: ",
     153             :                      Moose::stringify(_hist_old[_qp]),
     154             : 
     155             :                      "New history: ",
     156             :                      Moose::stringify(_hist[_qp]));
     157             :     else
     158           0 :       mooseException("NEML stress updated failed: ", e.message());
     159           0 :   }
     160             : 
     161             :   // Do more translation, now back to tensors
     162      209924 :   nemlToRankTwoTensor(s_np1, _stress[_qp]);
     163      209924 :   nemlToRankFourTensor(A_np1, _Jacobian_mult[_qp]);
     164             : 
     165             :   // Get the elastic strain
     166             :   try
     167             :   {
     168      209924 :     _model->elastic_strains(s_np1, T_np1, h_np1, estrain);
     169             :   }
     170           0 :   catch (const neml::NEMLError & e)
     171             :   {
     172           0 :     mooseError("Error in NEML call for elastic strain: ", e.message());
     173           0 :   }
     174             : 
     175             :   // Translate
     176      209924 :   nemlToRankTwoTensor(estrain, _elastic_strain[_qp]);
     177             : 
     178             :   // For EPP purposes calculate the inelastic strain
     179             :   double pstrain[6];
     180     1469468 :   for (unsigned int i = 0; i < 6; ++i)
     181     1259544 :     pstrain[i] = e_np1[i] - estrain[i];
     182             : 
     183      209924 :   nemlToRankTwoTensor(pstrain, _inelastic_strain[_qp]);
     184             : 
     185             :   // compute material timestep
     186      209924 :   if (_compute_dt)
     187             :   {
     188        5600 :     const auto increment = (_inelastic_strain[_qp] - (*_inelastic_strain_old)[_qp]).L2norm();
     189        5600 :     (*_material_dt)[_qp] =
     190        5600 :         increment > 0 ? _dt * _target_increment / increment : std::numeric_limits<Real>::max();
     191             :   }
     192             : 
     193             :   // Store dissipation
     194      209924 :   _energy[_qp] = u_np1;
     195      209924 :   _dissipation[_qp] = p_np1;
     196             :   // get damage index
     197      209924 :   if (_damage_index != nullptr)
     198       34094 :     (*_damage_index)[_qp] = _model->get_damage(h_np1);
     199      209924 : }
     200             : 
     201             : void
     202        2504 : NEMLStressBase::initQpStatefulProperties()
     203             : {
     204        2504 :   ComputeStressBase::initQpStatefulProperties();
     205             : 
     206             :   // Figure out initial history
     207        2504 :   _hist[_qp].resize(_model->nstore());
     208             : 
     209        2504 :   if (_model->nstore() > 0)
     210             :   {
     211             :     try
     212             :     {
     213        2504 :       _model->init_store(_hist[_qp].data());
     214             :     }
     215           0 :     catch (const neml::NEMLError & e)
     216             :     {
     217           0 :       mooseError("Error initializing NEML history: ", e.message());
     218           0 :     }
     219             :   }
     220        2504 :   _energy[_qp] = 0.0;
     221        2504 :   _dissipation[_qp] = 0.0;
     222             : 
     223        2504 :   if (_damage_index != nullptr)
     224         924 :     (*_damage_index)[_qp] = 0.0;
     225        2504 : }
     226             : 
     227             : void
     228    10426968 : NEMLStressBase::rankTwoTensorToNeml(const RankTwoTensor & in, double * const out)
     229             : {
     230    10426968 :   double inds[6][2] = {{0, 0}, {1, 1}, {2, 2}, {1, 2}, {0, 2}, {0, 1}};
     231    10426968 :   double mults[6] = {1.0, 1.0, 1.0, sqrt(2.0), sqrt(2.0), sqrt(2.0)};
     232             : 
     233    72988776 :   for (unsigned int i = 0; i < 6; ++i)
     234    62561808 :     out[i] = in(inds[i][0], inds[i][1]) * mults[i];
     235    10426968 : }
     236             : 
     237             : void
     238     7161236 : NEMLStressBase::nemlToRankTwoTensor(const double * const in, RankTwoTensor & out)
     239             : {
     240             :   static const unsigned int inds[6][2] = {{0, 0}, {1, 1}, {2, 2}, {1, 2}, {0, 2}, {0, 1}};
     241             :   static const double mults[6] = {
     242             :       1.0, 1.0, 1.0, 1.0 / std::sqrt(2.0), 1.0 / std::sqrt(2.0), 1.0 / std::sqrt(2.0)};
     243             : 
     244    50128652 :   for (unsigned int i = 0; i < 6; ++i)
     245             :   {
     246    42967416 :     out(inds[i][0], inds[i][1]) = in[i] * mults[i];
     247    42967416 :     out(inds[i][1], inds[i][0]) = in[i] * mults[i];
     248             :   }
     249     7161236 : }
     250             : 
     251             : void
     252      209924 : NEMLStressBase::nemlToRankFourTensor(const double * const in, RankFourTensor & out)
     253             : {
     254             :   static const unsigned int inds[6][2] = {{0, 0}, {1, 1}, {2, 2}, {1, 2}, {0, 2}, {0, 1}};
     255             :   static const double mults[6] = {1.0, 1.0, 1.0, 1.0 / sqrt(2.0), 1.0 / sqrt(2.0), 1.0 / sqrt(2.0)};
     256             : 
     257     1469468 :   for (unsigned int i = 0; i < 6; ++i)
     258     8816808 :     for (unsigned int j = 0; j < 6; ++j)
     259             :     {
     260     7557264 :       out(inds[i][0], inds[i][1], inds[j][0], inds[j][1]) = in[i * 6 + j] * (mults[i] * mults[j]);
     261     7557264 :       out(inds[i][1], inds[i][0], inds[j][0], inds[j][1]) = in[i * 6 + j] * (mults[i] * mults[j]);
     262     7557264 :       out(inds[i][0], inds[i][1], inds[j][1], inds[j][0]) = in[i * 6 + j] * (mults[i] * mults[j]);
     263     7557264 :       out(inds[i][1], inds[i][0], inds[j][1], inds[j][0]) = in[i * 6 + j] * (mults[i] * mults[j]);
     264             :     }
     265      209924 : }
     266             : 
     267             : #endif // NEML_ENABLED