GeneralizedRadialReturnStressUpdate.h

Go to the documentation of this file.

//* This file is part of the MOOSE framework
//* https://mooseframework.inl.gov
//*
//* All rights reserved, see COPYRIGHT for full restrictions
//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
//*
//* Licensed under LGPL 2.1, please see LICENSE for details
//* https://www.gnu.org/licenses/lgpl-2.1.html

#pragma once

#include "StressUpdateBase.h"
#include "GeneralizedReturnMappingSolution.h"
#include "MooseTypes.h"

#include "libmesh/ignore_warnings.h"
#include "Eigen/Dense"
#include "Eigen/Eigenvalues"
#include "libmesh/restore_warnings.h"

#include <limits>

namespace Eigen
{
namespace internal
{
template <>
struct cast_impl<ADReal, int>
{
  static inline int run(const ADReal & x) { return static_cast<int>(MetaPhysicL::raw_value(x)); }
};
} // namespace internal
} // namespace Eigen

// Instantiation of ADReal will be enabled shortly, once
// https://github.com/roystgnr/MetaPhysicL/issues/70 makes it to MOOSE
// typedef Eigen::Matrix<ADReal, 6, 6> AnisotropyMatrix;

typedef Eigen::Matrix<Real, 6, 6> AnisotropyMatrixReal;
typedef Eigen::Matrix<Real, 3, 3> AnisotropyMatrixRealBlock;

template <bool is_ad>
class GeneralizedRadialReturnStressUpdateTempl : public StressUpdateBaseTempl<is_ad>,
                                                 public GeneralizedReturnMappingSolutionTempl<is_ad>
{
public:
  static InputParameters validParams();

  GeneralizedRadialReturnStressUpdateTempl(const InputParameters & parameters);

  virtual void
  updateState(GenericRankTwoTensor<is_ad> & strain_increment,
              GenericRankTwoTensor<is_ad> & inelastic_strain_increment,
              const GenericRankTwoTensor<is_ad> & rotation_increment,
              GenericRankTwoTensor<is_ad> & stress_new,
              const RankTwoTensor & stress_old,
              const GenericRankFourTensor<is_ad> & elasticity_tensor,
              const RankTwoTensor & elastic_strain_old,
              bool /*compute_full_tangent_operator = false*/,
              RankFourTensor & /*tangent_operator = StressUpdateBaseTempl<is_ad>::_identityTensor*/)
      override;

  virtual Real
  computeReferenceResidual(const GenericDenseVector<is_ad> & effective_trial_stress,
                           const GenericDenseVector<is_ad> & stress_new,
                           const GenericReal<is_ad> & residual,
                           const GenericReal<is_ad> & scalar_effective_inelastic_strain) override;

  virtual GenericReal<is_ad> minimumPermissibleValue(
      const GenericDenseVector<is_ad> & /*effective_trial_stress*/) const override
  {
    return 0.0;
  }

  virtual GenericReal<is_ad>
  maximumPermissibleValue(const GenericDenseVector<is_ad> & effective_trial_stress) const override;

  virtual Real computeTimeStepLimit() override;

  virtual Real computeIntegrationErrorTimeStep() { return std::numeric_limits<Real>::max(); }

  bool requiresIsotropicTensor() override { return true; }

  bool isBlockDiagonal(const AnisotropyMatrixReal & A);

protected:
  using Material::_current_elem;
  using Material::_dt;
  using Material::_q_point;
  using Material::_qp;

  virtual void initQpStatefulProperties() override;

  void propagateQpStatefulPropertiesRadialReturn();

  virtual void
  computeStressInitialize(const GenericDenseVector<is_ad> & /*stress_dev*/,
                          const GenericDenseVector<is_ad> & /*stress*/,
                          const GenericRankFourTensor<is_ad> & /*elasticity_tensor*/) = 0;

  virtual GenericReal<is_ad> computeStressDerivative(const Real /*effective_trial_stress*/,
                                                     const Real /*scalar*/)
  {
    return 0.0;
  }

  virtual void computeStressFinalize(const GenericRankTwoTensor<is_ad> & inelasticStrainIncrement,
                                     const GenericReal<is_ad> & delta_gamma,
                                     GenericRankTwoTensor<is_ad> & stress,
                                     const GenericDenseVector<is_ad> & /*stress_dev*/,
                                     const GenericRankTwoTensor<is_ad> & stress_old,
                                     const GenericRankFourTensor<is_ad> & elasticity_tensor) = 0;
  virtual void computeStrainFinalize(GenericRankTwoTensor<is_ad> & inelasticStrainIncrement,
                                     const GenericRankTwoTensor<is_ad> & stress,
                                     const GenericDenseVector<is_ad> & stress_dev,
                                     const GenericReal<is_ad> & delta_gamma) = 0;

  /*
   * Method that determines the tangent calculation method. For anisotropic creep only models, the
   * tangent calculation method is ELASTIC for now
   */
  virtual TangentCalculationMethod getTangentCalculationMethod() override
  {
    return TangentCalculationMethod::ELASTIC;
  }

  void outputIterationSummary(std::stringstream * iter_output,
                              const unsigned int total_it) override;

  void computeTangentOperator(Real /*effective_trial_stress*/,
                              RankTwoTensor & /*stress_new*/,
                              bool /*compute_full_tangent_operator*/,
                              RankFourTensor & /*tangent_operator*/);

  GenericReal<is_ad> _three_shear_modulus;

  GenericMaterialProperty<Real, is_ad> & _effective_inelastic_strain;
  const MaterialProperty<Real> & _effective_inelastic_strain_old;

  MaterialProperty<Real> & _inelastic_strain_rate;
  const MaterialProperty<Real> & _inelastic_strain_rate_old;

  Real _max_inelastic_increment;

  Real _max_integration_error;

  Real _max_integration_error_time_step;

  const bool _use_transformation;
};

typedef GeneralizedRadialReturnStressUpdateTempl<false> GeneralizedRadialReturnStressUpdate;
typedef GeneralizedRadialReturnStressUpdateTempl<true> ADGeneralizedRadialReturnStressUpdate;