MultiParameterPlasticityStressUpdate.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 <array>

class MultiParameterPlasticityStressUpdate : public StressUpdateBase
{
public:
  static InputParameters validParams();

  MultiParameterPlasticityStressUpdate(const InputParameters & parameters,
                                       unsigned num_sp,
                                       unsigned num_yf,
                                       unsigned num_intnl);

protected:
  virtual void initQpStatefulProperties() override;
  using StressUpdateBase::updateState;
  virtual void updateState(RankTwoTensor & strain_increment,
                           RankTwoTensor & inelastic_strain_increment,
                           const RankTwoTensor & rotation_increment,
                           RankTwoTensor & stress_new,
                           const RankTwoTensor & stress_old,
                           const RankFourTensor & elasticity_tensor,
                           const RankTwoTensor & elastic_strain_old,
                           bool compute_full_tangent_operator,
                           RankFourTensor & tangent_operator) override;

  virtual void propagateQpStatefulProperties() override;

  virtual TangentCalculationMethod getTangentCalculationMethod() override
  {
    return TangentCalculationMethod::FULL;
  }

  constexpr static unsigned _tensor_dimensionality = 3;

  const unsigned _num_sp;

  const std::vector<Real> _definitely_ok_sp;

  std::vector<std::vector<Real>> _Eij;

  Real _En;

  std::vector<std::vector<Real>> _Cij;

  const unsigned _num_yf;

  const unsigned _num_intnl;

  const unsigned _max_nr_its;

  const bool _perform_finite_strain_rotations;

  const Real _smoothing_tol;

  const Real _smoothing_tol2;

  const Real _f_tol;

  const Real _f_tol2;

  const Real _min_step_size;

  bool _step_one;

  const bool _warn_about_precision_loss;

  MaterialProperty<RankTwoTensor> & _plastic_strain;

  const MaterialProperty<RankTwoTensor> & _plastic_strain_old;

  MaterialProperty<std::vector<Real>> & _intnl;

  const MaterialProperty<std::vector<Real>> & _intnl_old;

  MaterialProperty<std::vector<Real>> & _yf;

  MaterialProperty<Real> & _iter;

  MaterialProperty<Real> & _max_iter_used;

  const MaterialProperty<Real> & _max_iter_used_old;

  MaterialProperty<Real> & _linesearch_needed;

  struct yieldAndFlow
  {
    Real f;                                // yield function value
    std::vector<Real> df;                  // df/d(stress_param[i])
    std::vector<Real> df_di;               // df/d(intnl_variable[a])
    std::vector<Real> dg;                  // d(flow)/d(stress_param[i])
    std::vector<std::vector<Real>> d2g;    // d^2(flow)/d(sp[i])/d(sp[j])
    std::vector<std::vector<Real>> d2g_di; // d^2(flow)/d(sp[i])/d(intnl[a])

    yieldAndFlow() : yieldAndFlow(0, 0) {}

    yieldAndFlow(unsigned num_var, unsigned num_intnl)
      : f(0.0),
        df(num_var),
        df_di(num_intnl),
        dg(num_var),
        d2g(num_var, std::vector<Real>(num_var, 0.0)),
        d2g_di(num_var, std::vector<Real>(num_intnl, 0.0))
    {
    }

    // Zero everything without resizing/reallocating
    void reset()
    {
      f = 0.0;
      std::fill(df.begin(), df.end(), 0.0);
      std::fill(df_di.begin(), df_di.end(), 0.0);
      std::fill(dg.begin(), dg.end(), 0.0);
      for (auto & row : d2g)
        std::fill(row.begin(), row.end(), 0.0);
      for (auto & row : d2g_di)
        std::fill(row.begin(), row.end(), 0.0);
    }

    // this may be involved in the smoothing of a group of yield functions
    bool operator<(const yieldAndFlow & fd) const { return f < fd.f; }
  };

  Real yieldF(const std::vector<Real> & stress_params, const std::vector<Real> & intnl) const;

  Real yieldF(const std::vector<Real> & yfs) const;

  Real ismoother(Real f_diff) const;

  Real smoother(Real f_diff) const;

  Real dsmoother(Real f_diff) const;

  yieldAndFlow smoothAllQuantities(const std::vector<Real> & stress_params,
                                   const std::vector<Real> & intnl) const;

  int lineSearch(Real & res2,
                 std::vector<Real> & stress_params,
                 Real & gaE,
                 const std::vector<Real> & trial_stress_params,
                 yieldAndFlow & smoothed_q,
                 const std::vector<Real> & intnl_ok,
                 std::vector<Real> & intnl,
                 std::vector<Real> & rhs,
                 Real & linesearch_needed) const;

  int nrStep(const yieldAndFlow & smoothed_q,
             const std::vector<Real> & trial_stress_params,
             const std::vector<Real> & stress_params,
             const std::vector<Real> & intnl,
             Real gaE,
             std::vector<Real> & rhs) const;

  Real calculateRes2(const std::vector<Real> & rhs) const;

  void calculateRHS(const std::vector<Real> & trial_stress_params,
                    const std::vector<Real> & stress_params,
                    Real gaE,
                    const yieldAndFlow & smoothed_q,
                    std::vector<Real> & rhs) const;

  void dnRHSdVar(const yieldAndFlow & smoothed_q,
                 const std::vector<std::vector<Real>> & dintnl,
                 const std::vector<Real> & stress_params,
                 Real gaE,
                 std::vector<double> & jac) const;

  virtual void errorHandler(const std::string & message) const;

  virtual void yieldFunctionValuesV(const std::vector<Real> & stress_params,
                                    const std::vector<Real> & intnl,
                                    std::vector<Real> & yf) const = 0;

  virtual void computeAllQV(const std::vector<Real> & stress_params,
                            const std::vector<Real> & intnl,
                            std::vector<yieldAndFlow> & all_q) const = 0;

  virtual void preReturnMapV(const std::vector<Real> & trial_stress_params,
                             const RankTwoTensor & stress_trial,
                             const std::vector<Real> & intnl_old,
                             const std::vector<Real> & yf,
                             const RankFourTensor & Eijkl);

  virtual void initializeVarsV(const std::vector<Real> & trial_stress_params,
                               const std::vector<Real> & intnl_old,
                               std::vector<Real> & stress_params,
                               Real & gaE,
                               std::vector<Real> & intnl) const;

  virtual void setIntnlValuesV(const std::vector<Real> & trial_stress_params,
                               const std::vector<Real> & current_stress_params,
                               const std::vector<Real> & intnl_old,
                               std::vector<Real> & intnl) const = 0;

  virtual void setIntnlDerivativesV(const std::vector<Real> & trial_stress_params,
                                    const std::vector<Real> & current_stress_params,
                                    const std::vector<Real> & intnl,
                                    std::vector<std::vector<Real>> & dintnl) const = 0;

  virtual void computeStressParams(const RankTwoTensor & stress,
                                   std::vector<Real> & stress_params) const = 0;

  virtual void initializeReturnProcess();

  virtual void finalizeReturnProcess(const RankTwoTensor & rotation_increment);

  virtual void setStressAfterReturnV(const RankTwoTensor & stress_trial,
                                     const std::vector<Real> & stress_params,
                                     Real gaE,
                                     const std::vector<Real> & intnl,
                                     const yieldAndFlow & smoothed_q,
                                     const RankFourTensor & Eijkl,
                                     RankTwoTensor & stress) const = 0;

  virtual void setInelasticStrainIncrementAfterReturn(const RankTwoTensor & stress_trial,
                                                      Real gaE,
                                                      const yieldAndFlow & smoothed_q,
                                                      const RankFourTensor & elasticity_tensor,
                                                      const RankTwoTensor & returned_stress,
                                                      RankTwoTensor & inelastic_strain_increment);

  virtual void consistentTangentOperatorV(const RankTwoTensor & stress_trial,
                                          const std::vector<Real> & trial_stress_params,
                                          const RankTwoTensor & stress,
                                          const std::vector<Real> & stress_params,
                                          Real gaE,
                                          const yieldAndFlow & smoothed_q,
                                          const RankFourTensor & Eijkl,
                                          bool compute_full_tangent_operator,
                                          const std::vector<std::vector<Real>> & dvar_dtrial,
                                          RankFourTensor & cto);

  virtual void dstressparam_dstress(const RankTwoTensor & stress,
                                    std::vector<RankTwoTensor> & dsp) const = 0;

  virtual void d2stressparam_dstress(const RankTwoTensor & stress,
                                     std::vector<RankFourTensor> & d2sp) const = 0;

  virtual void setEffectiveElasticity(const RankFourTensor & Eijkl) = 0;

  void dVardTrial(bool elastic_only,
                  const std::vector<Real> & trial_stress_params,
                  const std::vector<Real> & stress_params,
                  Real gaE,
                  const std::vector<Real> & intnl,
                  const yieldAndFlow & smoothed_q,
                  Real step_size,
                  bool compute_full_tangent_operator,
                  std::vector<std::vector<Real>> & dvar_dtrial) const;

  bool precisionLoss(const std::vector<Real> & solution,
                     const std::vector<Real> & stress_params,
                     Real gaE) const;

private:
  std::vector<Real> _trial_sp;

  RankTwoTensor _stress_trial;

  std::vector<Real> _rhs;

  std::vector<std::vector<Real>> _dvar_dtrial;

  std::vector<Real> _ok_sp;

  std::vector<Real> _ok_intnl;

  std::vector<Real> _del_stress_params;

  std::vector<Real> _current_sp;

  std::vector<Real> _current_intnl;

  yieldAndFlow _smoothed_q;

  mutable std::vector<RankTwoTensor> _dsp_scratch;

  mutable std::vector<RankTwoTensor> _dsp_trial_scratch;

  mutable std::vector<RankFourTensor> _d2sp_scratch;

  mutable std::vector<Real> _yfs_scratch;

  mutable std::vector<Real> _sp_params_old_scratch;

  mutable std::vector<Real> _delta_nr_params_scratch;

  mutable std::vector<std::vector<Real>> _dintnl_scratch;

  mutable std::vector<double> _jac_scratch;

  mutable std::vector<PetscBLASInt> _ipiv_scratch;

  mutable std::vector<yieldAndFlow> _all_q_scratch;

private:
  enum class SmootherFunctionType
  {
    cos,
    poly1,
    poly2,
    poly3
  } _smoother_function_type;
};