Line data    Source code

       1              : // Copyright 2024, UChicago Argonne, LLC
       2              : // All Rights Reserved
       3              : // Software Name: NEML2 -- the New Engineering material Model Library, version 2
       4              : // By: Argonne National Laboratory
       5              : // OPEN SOURCE LICENSE (MIT)
       6              : //
       7              : // Permission is hereby granted, free of charge, to any person obtaining a copy
       8              : // of this software and associated documentation files (the "Software"), to deal
       9              : // in the Software without restriction, including without limitation the rights
      10              : // to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
      11              : // copies of the Software, and to permit persons to whom the Software is
      12              : // furnished to do so, subject to the following conditions:
      13              : //
      14              : // The above copyright notice and this permission notice shall be included in
      15              : // all copies or substantial portions of the Software.
      16              : //
      17              : // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
      18              : // IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
      19              : // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
      20              : // AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
      21              : // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
      22              : // OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
      23              : // THE SOFTWARE.
      24              : 
      25              : #include "neml2/models/solid_mechanics/crystal_plasticity/OrientationRate.h"
      26              : 
      27              : #include "neml2/tensors/WR2.h"
      28              : #include "neml2/tensors/SR2.h"
      29              : #include "neml2/tensors/R4.h"
      30              : #include "neml2/tensors/WWR4.h"
      31              : #include "neml2/tensors/WSR4.h"
      32              : 
      33              : namespace neml2
      34              : {
      35              : register_NEML2_object(OrientationRate);
      36              : 
      37              : OptionSet
      38            2 : OrientationRate::expected_options()
      39              : {
      40            2 :   OptionSet options = Model::expected_options();
      41              : 
      42            2 :   options.doc() =
      43              :       "Defines the rate of the crystal orientations as a spin given by \\f$ \\Omega^e = "
      44              :       "w - w^p - \\varepsilon d^p + d^p \\varepsilon \\f$ where \\f$ \\Omega^e = \\dot{Q} Q^T "
      45              :       "\\f$, \\f$ Q \\f$ is the orientation, \\f$ w \\f$ is the vorticity, \\f$ w^p \\f$ is the "
      46              :       "plastic vorticity, \\f$ d^p \\f$ is the plastic deformation rate, and \\f$ \\varepsilon "
      47            2 :       "\\f$ is the elastic stretch.";
      48              : 
      49            6 :   options.set_output("orientation_rate") = VariableName(STATE, "orientation_rate");
      50            2 :   options.set("orientation_rate").doc() = "The name of the orientation rate (spin)";
      51              : 
      52            6 :   options.set_input("elastic_strain") = VariableName(STATE, "elastic_strain");
      53            2 :   options.set("elastic_strain").doc() = "The name of the elastic strain tensor";
      54              : 
      55            6 :   options.set_input("vorticity") = VariableName(FORCES, "vorticity");
      56            2 :   options.set("vorticity").doc() = "The name of the voriticty tensor";
      57              : 
      58            4 :   options.set_input("plastic_deformation_rate") =
      59            6 :       VariableName(STATE, "internal", "plastic_deformation_rate");
      60            2 :   options.set("plastic_deformation_rate").doc() = "The name of the plastic deformation rate";
      61              : 
      62            6 :   options.set_input("plastic_vorticity") = VariableName(STATE, "internal", "plastic_vorticity");
      63            2 :   options.set("plastic_vorticity").doc() = "The name of the plastic vorticity";
      64            2 :   return options;
      65            0 : }
      66              : 
      67            2 : OrientationRate::OrientationRate(const OptionSet & options)
      68              :   : Model(options),
      69            2 :     _R_dot(declare_output_variable<WR2>("orientation_rate")),
      70            2 :     _e(declare_input_variable<SR2>("elastic_strain")),
      71            2 :     _w(declare_input_variable<WR2>("vorticity")),
      72            2 :     _dp(declare_input_variable<SR2>("plastic_deformation_rate")),
      73            4 :     _wp(declare_input_variable<WR2>("plastic_vorticity"))
      74              : {
      75            2 : }
      76              : 
      77              : WR2
      78            3 : multiply_and_make_skew(const SR2 & a, const SR2 & b)
      79              : {
      80            3 :   auto A = R2(a);
      81            3 :   auto B = R2(b);
      82              : 
      83            6 :   return WR2(A * B - B * A);
      84            3 : }
      85              : 
      86              : WSR4
      87            2 : d_multiply_and_make_skew_d_first(const SR2 & b)
      88              : {
      89            2 :   auto I = R2::identity(b.options());
      90            2 :   auto B = R2(b);
      91              :   return WSR4(
      92           16 :       R4(at::einsum("...ia,...bj->...ijab", {I, B}) - at::einsum("...ia,...jb->...ijab", {B, I})));
      93            6 : }
      94              : 
      95              : WSR4
      96            2 : d_multiply_and_make_skew_d_second(const SR2 & a)
      97              : {
      98            2 :   auto I = R2::identity(a.options());
      99            2 :   auto A = R2(a);
     100              :   return WSR4(
     101           16 :       R4(at::einsum("...ia,...jb->...ijab", {A, I}) - at::einsum("...ia,...bj->...ijab", {I, A})));
     102            6 : }
     103              : 
     104              : void
     105            4 : OrientationRate::set_value(bool out, bool dout_din, bool /*d2out_din2*/)
     106              : {
     107            4 :   if (out)
     108            3 :     _R_dot = _w - _wp + multiply_and_make_skew(SR2(_dp), SR2(_e));
     109              : 
     110            4 :   if (dout_din)
     111              :   {
     112            2 :     const auto I = WWR4::identity(_w.options());
     113              : 
     114            2 :     if (_e.is_dependent())
     115            2 :       _R_dot.d(_e) = d_multiply_and_make_skew_d_second(SR2(_dp));
     116              : 
     117            2 :     if (_w.is_dependent())
     118            1 :       _R_dot.d(_w) = I;
     119              : 
     120            2 :     if (_dp.is_dependent())
     121            2 :       _R_dot.d(_dp) = d_multiply_and_make_skew_d_first(SR2(_e));
     122              : 
     123            2 :     if (_wp.is_dependent())
     124            2 :       _R_dot.d(_wp) = -I;
     125            2 :   }
     126            4 : }
     127              : } // namespace neml2