LCOV - code coverage report
Current view: top level - src/Model/Solid/Affine - J2.H (source / functions) Coverage Total Hit
Test: coverage_merged.info Lines: 81.2 % 64 52
Test Date: 2025-02-27 04:17:48 Functions: 100.0 % 9 9

            Line data    Source code
       1              : //
       2              : // This models an isotropic elastic-perfectly-plastic, non-time-dependent solid model.
       3              : // 
       4              : // The energy and derivatives are:
       5              : //
       6              : // .. math::
       7              : //    :nowrap:
       8              : //
       9              : //    \begin{gather}
      10              : //    W = \frac{1}{2}(\varepsilon - \varepsilon_p):\mathbb{C}(\varepsilon-\varepsilon_p) \\
      11              : //    DW = \mathbb{C}(\varepsilon-\varepsilon_p) \\
      12              : //    DDW = \mathbb{C}
      13              : //    \end{gather}
      14              : //
      15              : // where :math:`\mathbb{C}` is an isotropic :code:`Set::Matrix4` and :math:`\varepsilon_p` is 
      16              : // is stored in the :code:`F0` eigenstrain.
      17              : //
      18              : // The plastic strain is evolved according to the following:
      19              : //
      20              : // #. Calculate the deviatoric stress :math:`\sigma_v=\sigma - \frac{1}{3}tr(\sigma)\mathbf{I}`
      21              : // #. Calculate :math:`J_2=\sqrt{\frac{3}{2}\sigma_v:\sigma_v}`
      22              : // #. If :math:`J_2<\sigma_0` then quit - no plasticity occurs
      23              : // #. Calculate :math:`\Delta\sigma = (1-\frac{\sigma_0}{J_2})`, which projects the stress 
      24              : //    back on the yield surface.
      25              : // #. Convert to change in plastic strain, :math:`\Delta\varepsilon=\mathbb{C}^{-1}\Delta\sigma`
      26              : // #. Update plastic strain: :math:`\varepsilon_p += \Delta\varepsilon`
      27              : //
      28              : // Notes:
      29              : //
      30              : // * This does not implement any kind of hardening model. Rate hardening, isotropic hardening,
      31              : //   and kinematic hardening have yet to be implemneted.
      32              : //
      33              : #ifndef MODEL_SOLID_AFFINE_J2_H_
      34              : #define MODEL_SOLID_AFFINE_J2_H_
      35              : 
      36              : #include "AMReX.H"
      37              : #include "IO/ParmParse.H"
      38              : #include "Model/Solid/Affine/Isotropic.H"
      39              : #include "Set/Set.H"
      40              : #include <cmath>
      41              : 
      42              : namespace Model
      43              : {
      44              : namespace Solid
      45              : {
      46              : namespace Affine
      47              : {
      48              : class J2 : public Isotropic
      49              : {
      50              : public:
      51              : 
      52       244445 :     J2() {};
      53              : 
      54       178966 :     void Advance(Set::Scalar dt, Set::Matrix eps, Set::Matrix /*sig*/, Set::Scalar /*time*/) override
      55              :     {
      56       178966 :         Set::Matrix sigma = DW(eps);
      57       178966 :         Set::Matrix sig_dev = sigma - sigma.trace() * Set::Matrix::Identity() / 3.0;
      58       178966 :         if ((hardening<0.0) || (ratecoeff<0.0)) // rate independent
      59              :         {
      60       178966 :             Set::Scalar J2 = sqrt(1.5 * (sig_dev.transpose() * sig_dev).trace());
      61       178966 :             if (J2 < sigma0) return; // No plasticity happens.
      62        13552 :             Set::Matrix dsigp = (1. - sigma0/J2) * sig_dev;
      63        27104 :             F0 += ddw.Inverse() * dsigp; 
      64        13552 :         }
      65              :         else // rate-dependent
      66              :         {
      67            0 :             Set::Matrix eps_dev = eps - eps.trace() * Set::Matrix::Identity() / (double)AMREX_SPACEDIM;
      68            0 :             Set::Scalar m_mu = ddw.Mu();
      69            0 :             Set::Matrix zeta_trial = sig_dev + 2.0*m_mu*eps_dev;
      70            0 :             Set::Scalar f_trial = zeta_trial.norm() - sqrt(1.0 - 1.0/(double)AMREX_SPACEDIM)*(sigma0 + hardening*alpha);
      71            0 :             if (f_trial <= 0.0) return;
      72            0 :             Set::Matrix n_new = zeta_trial / zeta_trial.norm();
      73            0 :             Set::Scalar dGamma = f_trial / (2.0*m_mu*(1 + hardening/(3.0*m_mu)));
      74            0 :             dGamma *= ratecoeff * dt;
      75            0 :             alpha = alpha + sqrt(1.0 - 1.0/(double)AMREX_SPACEDIM)*dGamma;
      76            0 :             F0 += dGamma * n_new;
      77              :         }
      78              :     }
      79              : 
      80              :     // Parameters
      81              :     Set::Scalar sigma0 = 1.0;
      82              :     Set::Scalar hardening = NAN;
      83              :     Set::Scalar ratecoeff = NAN;
      84              : 
      85              :     // Internal variables to track state
      86              :     Set::Scalar alpha = 0.0;
      87              : 
      88          706 :     static J2 Zero()
      89              :     {
      90          706 :         J2 ret;
      91          706 :         ret.ddw = Set::Matrix4<AMREX_SPACEDIM,Set::Sym::Isotropic>::Zero();
      92          706 :         ret.F0 = Set::Matrix::Zero();
      93          706 :         ret.sigma0 = 0.0;
      94          706 :         ret.hardening = 0.0;
      95          706 :         ret.ratecoeff = 0.0;
      96          706 :         return ret;
      97            0 :     }
      98              : 
      99              :     // See also inputs to :ref:`Model::Solid::Affine::Isotropic`
     100            3 :     static void Parse(J2 & value, IO::ParmParse & pp)
     101              :     {
     102            3 :         Isotropic::Parse(value,pp);
     103              :         // J2 Yield criterion
     104           18 :         pp_query_default("sigma0",value.sigma0,1.0);
     105              :         // Hardening coefficient (negative value disables rate hardening)
     106           18 :         pp_query_default("hardening",value.hardening,-1.0);
     107              :         // Rate coefficient (negative value disables rate hardening)
     108           15 :         pp_query_default("ratecoeff",value.ratecoeff,-1.0);
     109            3 :     }
     110              : 
     111              :     #define OP_CLASS J2
     112              :     #define OP_VARS  X(ddw) X(F0) X(sigma0) X(hardening) X(ratecoeff)
     113              :     #include "Model/Solid/InClassOperators.H"
     114              : };
     115              : #include "Model/Solid/ExtClassOperators.H"
     116              : }
     117              : }
     118              : }
     119              : 
     120              : 
     121              : template<>
     122              : AMREX_ATTRIBUTE_WEAK
     123           48 : int Set::Field<Model::Solid::Affine::J2>::NComp() const 
     124              : {
     125           48 :     return AMREX_SPACEDIM*AMREX_SPACEDIM;
     126              : }
     127              : 
     128              : template<>
     129              : AMREX_ATTRIBUTE_WEAK
     130           30 : std::string Set::Field<Model::Solid::Affine::J2>::Name(int i) const 
     131              : {
     132              : #if AMREX_SPACEDIM==2
     133           12 :     if (i==0) return name + ".epsp_xx";
     134            9 :     if (i==1) return name + ".epsp_xy";
     135            6 :     if (i==2) return name + ".epsp_yx";
     136            3 :     if (i==3) return name + ".epsp_yy";
     137              : #endif
     138              : #if AMREX_SPACEDIM==3
     139           18 :     if (i==0) return name + ".epsp_xx";
     140           16 :     if (i==1) return name + ".epsp_xy";
     141           14 :     if (i==2) return name + ".epsp_xz";
     142           12 :     if (i==3) return name + ".epsp_yx";
     143           10 :     if (i==4) return name + ".epsp_yy";
     144            8 :     if (i==5) return name + ".epsp_yz";
     145            6 :     if (i==6) return name + ".epsp_zx";
     146            4 :     if (i==7) return name + ".epsp_zy";
     147            2 :     if (i==8) return name + ".epsp_zz";
     148              : #endif
     149            0 :     return name;
     150              : }
     151              : 
     152              : template<>
     153              : AMREX_ATTRIBUTE_WEAK
     154            8 : void Set::Field<Model::Solid::Affine::J2>::Copy(int a_lev, amrex::MultiFab &a_dst, int a_dstcomp, int a_nghost) const
     155              : {
     156           31 :     for (amrex::MFIter mfi(a_dst, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
     157              :     {
     158           23 :         const amrex::Box& bx = mfi.growntilebox(amrex::IntVect(a_nghost));
     159           23 :         if (bx.ok())
     160              :         {
     161           23 :             amrex::Array4<const Model::Solid::Affine::J2> const & src = ((*this)[a_lev])->array(mfi);
     162           23 :             amrex::Array4<Set::Scalar> const & dst = a_dst.array(mfi);
     163          125 :             for (int n = 0; n < AMREX_SPACEDIM*AMREX_SPACEDIM; n++)
     164              :             {
     165          102 :                 amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
     166        27068 :                     dst(i,j,k,a_dstcomp + n) = src(i,j,k).F0(n/AMREX_SPACEDIM,n%AMREX_SPACEDIM);
     167        13534 :                 });
     168              :             }
     169              :         }
     170            8 :     }    
     171            8 : }
     172              : 
     173              : 
     174              : #endif
        

Generated by: LCOV version 2.0-1