Skip to content

Commit

Permalink
Assert history material property sizing (idaholab#193)
Browse files Browse the repository at this point in the history
  • Loading branch information
dschwen committed Mar 31, 2021
1 parent eb87310 commit 113f4d3
Show file tree
Hide file tree
Showing 2 changed files with 33 additions and 30 deletions.
31 changes: 17 additions & 14 deletions include/materials/NEMLStressBase.h
Original file line number Diff line number Diff line change
Expand Up @@ -35,30 +35,33 @@ class NEMLStressBase : public ComputeStressBase
protected:
/// NEML model
std::unique_ptr<neml::NEMLModel> _model;
/// History variables used by NEML model
///@{

///@{ History variables used by NEML model
MaterialProperty<std::vector<Real>> & _hist;
const MaterialProperty<std::vector<Real>> & _hist_old;
///@}

/// Old mechanical strain
const MaterialProperty<RankTwoTensor> & _mechanical_strain_old;

/// Old stress
const MaterialProperty<RankTwoTensor> & _stress_old;
/// Strain energy
///@{

///@{ Strain energy
MaterialProperty<Real> & _energy;
const MaterialProperty<Real> & _energy_old;
///@}
/// Dissipation
///@{

///@{ Dissipation
MaterialProperty<Real> & _dissipation;
const MaterialProperty<Real> & _dissipation_old;
///@}
/// Coupled temperature variable (defaults to zero if not specified)
///@{

///@{ Coupled temperature variable (defaults to zero if not specified)
const VariableValue & _temperature;
const VariableValue & _temperature_old;
///@}

/// Inelastic strain tensor
MaterialProperty<RankTwoTensor> & _inelastic_strain;

Expand All @@ -72,20 +75,20 @@ class NEMLStressBase : public ComputeStressBase
* format.
* @param in RankTwoTensor to be translated
* @param out NEML vector output
**/
void RankTwoTensorToNeml(const RankTwoTensor & in, double * const out);
*/
void rankTwoTensorToNeml(const RankTwoTensor & in, double * const out);

/**
* Translates a NEML tensor stored in vector format to a RankTwoTensor.
* @param in NEML vector to be translated
* @param out RankTwoTensor output
**/
void NemlToRankTwoTensor(const double * const in, RankTwoTensor & out);
*/
void nemlToRankTwoTensor(const double * const in, RankTwoTensor & out);

/**
* Translates a NEML elasticity tensor to a RankFourTensor.
* @param in NEML elasticity tensor to be translated
* @param out RankFourTensor output
**/
void NemlToRankFourTensor(const double * const in, RankFourTensor & out);
*/
void nemlToRankFourTensor(const double * const in, RankFourTensor & out);
};
32 changes: 16 additions & 16 deletions src/materials/NEMLStressBase.C
Original file line number Diff line number Diff line change
Expand Up @@ -69,21 +69,23 @@ NEMLStressBase::computeQpStress()
// First do some declaration and translation
double s_np1[6];
double s_n[6];
RankTwoTensorToNeml(_stress_old[_qp], s_n);
rankTwoTensorToNeml(_stress_old[_qp], s_n);

double e_np1[6];
RankTwoTensorToNeml(_mechanical_strain[_qp], e_np1);
rankTwoTensorToNeml(_mechanical_strain[_qp], e_np1);
double e_n[6];
RankTwoTensorToNeml(_mechanical_strain_old[_qp], e_n);
rankTwoTensorToNeml(_mechanical_strain_old[_qp], e_n);

double t_np1 = _t;
double t_n = _t - _dt;

double T_np1 = _temperature[_qp];
double T_n = _temperature_old[_qp];

double * const h_np1 = (_model->nhist() > 0 ? &(_hist[_qp][0]) : nullptr);
const double * const h_n = (_model->nhist() > 0 ? &(_hist_old[_qp][0]) : nullptr);
mooseAssert(_model->nhist() == _hist[_qp].size(), "History data storage size mismatch");
double * const h_np1 = (_model->nhist() > 0 ? _hist[_qp].data() : nullptr);
mooseAssert(_model->nhist() == _hist_old[_qp].size(), "History data storage size mismatch");
const double * const h_n = (_model->nhist() > 0 ? _hist_old[_qp].data() : nullptr);

double A_np1[36];

Expand All @@ -95,18 +97,16 @@ NEMLStressBase::computeQpStress()

double estrain[6];

int ier;

// Actually call the update
ier = _model->update_sd(
auto ier = _model->update_sd(
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);

if (ier != neml::SUCCESS)
mooseException("NEML stress update failed!");

// Do more translation, now back to tensors
NemlToRankTwoTensor(s_np1, _stress[_qp]);
NemlToRankFourTensor(A_np1, _Jacobian_mult[_qp]);
nemlToRankTwoTensor(s_np1, _stress[_qp]);
nemlToRankFourTensor(A_np1, _Jacobian_mult[_qp]);

// Get the elastic strain
ier = _model->elastic_strains(s_np1, T_np1, h_np1, estrain);
Expand All @@ -115,14 +115,14 @@ NEMLStressBase::computeQpStress()
mooseError("Error in NEML call for elastic strain!");

// Translate
NemlToRankTwoTensor(estrain, _elastic_strain[_qp]);
nemlToRankTwoTensor(estrain, _elastic_strain[_qp]);

// For EPP purposes calculate the inelastic strain
double pstrain[6];
for (unsigned int i = 0; i < 6; ++i)
pstrain[i] = e_np1[i] - estrain[i];

NemlToRankTwoTensor(pstrain, _inelastic_strain[_qp]);
nemlToRankTwoTensor(pstrain, _inelastic_strain[_qp]);

// compute material timestep
if (_compute_dt)
Expand All @@ -147,7 +147,7 @@ NEMLStressBase::initQpStatefulProperties()

if (_model->nhist() > 0)
{
int ier = _model->init_hist(&(_hist[_qp][0]));
int ier = _model->init_hist(_hist[_qp].data());
if (ier != neml::SUCCESS)
mooseError("Error initializing NEML history!");
}
Expand All @@ -157,7 +157,7 @@ NEMLStressBase::initQpStatefulProperties()
}

void
NEMLStressBase::RankTwoTensorToNeml(const RankTwoTensor & in, double * const out)
NEMLStressBase::rankTwoTensorToNeml(const RankTwoTensor & in, double * const out)
{
double inds[6][2] = {{0, 0}, {1, 1}, {2, 2}, {1, 2}, {0, 2}, {0, 1}};
double mults[6] = {1.0, 1.0, 1.0, sqrt(2.0), sqrt(2.0), sqrt(2.0)};
Expand All @@ -167,7 +167,7 @@ NEMLStressBase::RankTwoTensorToNeml(const RankTwoTensor & in, double * const out
}

void
NEMLStressBase::NemlToRankTwoTensor(const double * const in, RankTwoTensor & out)
NEMLStressBase::nemlToRankTwoTensor(const double * const in, RankTwoTensor & out)
{
static const unsigned int inds[6][2] = {{0, 0}, {1, 1}, {2, 2}, {1, 2}, {0, 2}, {0, 1}};
static const double mults[6] = {
Expand All @@ -181,7 +181,7 @@ NEMLStressBase::NemlToRankTwoTensor(const double * const in, RankTwoTensor & out
}

void
NEMLStressBase::NemlToRankFourTensor(const double * const in, RankFourTensor & out)
NEMLStressBase::nemlToRankFourTensor(const double * const in, RankFourTensor & out)
{
static const unsigned int inds[6][2] = {{0, 0}, {1, 1}, {2, 2}, {1, 2}, {0, 2}, {0, 1}};
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)};
Expand Down

0 comments on commit 113f4d3

Please sign in to comment.