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