Added an option to log transform the critical work relation

doc/sphinx/damage/work.rst

@@ -22,6 +22,11 @@ In principle, the critical work might also depend on temperatures.
 However, at least for one material (Alloy 617) the temperature dependence
 is relatively negligible.
 
+If requested the model will first take the log of the work rate
+before passing it to the :math:`W_{crit}` function and uses :math:`10^f` of the
+returned value.  This means the user is providing the function on a log-log
+scale.
+
 Parameters
 ----------
 
@@ -32,7 +37,8 @@ Parameters
    ``elastic``, :cpp:class:`neml::LinearElasticModel`, Elasticity model, No
    ``Wcrit``, :cpp:class:`neml::Interpolate`, Critical work, No
    ``n``, :code:`double`, Damage exponent, No
-   ``eps``, :code:`double`, Numerical offset``, ``1.0e-30``
+   ``eps``, :code:`double`, Numerical offset, ``1.0e-30``
+   ``log``, :code:`bool`, Log transform the work to failure relation, ``false``
 
 Class description
 -----------------

include/damage.h

@@ -653,11 +653,18 @@ class NEML_EXPORT WorkDamage: public ScalarDamage {
       const double * const stress_np1, const double * const stress_n,
       double T_np1, double T_n, double t_np1, double t_n,
       double d_np1, double d_n) const;
+
+  /// Get the critical work, accounting for log scaling if requested
+  double Wcrit(double Wdot) const;
+
+  /// Get the derivative of the critical work, accounting for log scaling
+  double dWcrit(double Wdot) const;
 
  protected:
   std::shared_ptr<Interpolate> Wcrit_;
   double n_;
   double eps_;
+  bool log_;
 };
 
 static Register<WorkDamage> regWorkDamage;

src/damage.cxx

@@ -1274,7 +1274,8 @@ WorkDamage::WorkDamage(ParameterSet & params) :
       ScalarDamage(params), 
       Wcrit_(params.get_object_parameter<Interpolate>("Wcrit")),
       n_(params.get_parameter<double>("n")),
-      eps_(params.get_parameter<double>("eps"))
+      eps_(params.get_parameter<double>("eps")),
+      log_(params.get_parameter<bool>("log"))
 {
 }
 
@@ -1291,6 +1292,7 @@ ParameterSet WorkDamage::parameters()
   pset.add_parameter<NEMLObject>("Wcrit");
   pset.add_parameter<double>("n");
   pset.add_optional_parameter<double>("eps", 1e-30);
+  pset.add_optional_parameter<bool>("log", false);
 
   return pset;
 }
@@ -1321,7 +1323,7 @@ void WorkDamage::damage(double d_np1, double d_n,
 
   double dt = t_np1 - t_n;
 
-  double val = Wcrit_->value(wrate);
+  double val = Wcrit(wrate);
 
   *dd = d_n + n_ * std::pow(std::fabs(d_np1), (n_-1.0)/n_) * 
       wrate * dt / val;
@@ -1346,8 +1348,8 @@ void WorkDamage::ddamage_dd(double d_np1, double d_n,
     return;
   }
 
-  double val = Wcrit_->value(wrate);
-  double deriv = Wcrit_->derivative(wrate);
+  double val = Wcrit(wrate);
+  double deriv = dWcrit(wrate);
   double dt = t_np1 - t_n;
 
   double S[36];
@@ -1381,8 +1383,8 @@ void WorkDamage::ddamage_de(double d_np1, double d_n,
     return;
   }
 
-  double val = Wcrit_->value(wrate);
-  double dval = Wcrit_->derivative(wrate);
+  double val = Wcrit(wrate);
+  double dval = dWcrit(wrate);
 
   double fact = n_ * std::pow(std::fabs(d_np1), (n_-1.0)/n_) / val * 
       (1.0 - wrate / val * dval) * (1.0 - std::fabs(d_np1));
@@ -1407,8 +1409,8 @@ void WorkDamage::ddamage_ds(double d_np1, double d_n,
     return;
   }
 
-  double val = Wcrit_->value(wrate);
-  double dval = Wcrit_->derivative(wrate);
+  double val = Wcrit(wrate);
+  double dval = dWcrit(wrate);
 
   double S[36];
   elastic_->S(T_np1, S);
@@ -1466,6 +1468,23 @@ double WorkDamage::workrate(
   return fabs(dot_vec(stress_np1, dp, 6) / dt * (1.0 - d_np1));
 }
 
+double WorkDamage::Wcrit(double Wdot) const
+{
+  if (log_)
+    return std::pow(10.0, Wcrit_->value(std::log10(Wdot)));
+
+  return Wcrit_->value(Wdot);
+}
+
+double WorkDamage::dWcrit(double Wdot) const
+{
+  if (log_)
+    return std::pow(10.0, Wcrit_->value(std::log10(Wdot))) * 
+        Wcrit_->derivative(std::log10(Wdot)) / Wdot;
+
+  return Wcrit_->derivative(Wdot);
+}
+
 PowerLawDamage::PowerLawDamage(ParameterSet & params) :
       StandardScalarDamage(params), 
       A_(params.get_object_parameter<Interpolate>("A")),

test/test_damage.py

@@ -275,6 +275,80 @@ def test_definition(self):
 
     self.assertAlmostEqual(damage, should)
 
+class TestWorkDamageLog(unittest.TestCase, CommonScalarDamageModel,
+    CommonDamagedModel):
+  def setUp(self):
+    self.E = 92000.0
+    self.nu = 0.3
+
+    self.s0 = 180.0
+    self.Kp = 1000.0
+    self.H = 1000.0
+
+    self.elastic = elasticity.IsotropicLinearElasticModel(self.E, "youngs",
+        self.nu, "poissons")
+
+    surface = surfaces.IsoKinJ2()
+    iso = hardening.LinearIsotropicHardeningRule(self.s0, self.Kp)
+    kin = hardening.LinearKinematicHardeningRule(self.H)
+    hrule = hardening.CombinedHardeningRule(iso, kin)
+
+    flow = ri_flow.RateIndependentAssociativeFlow(surface, hrule)
+
+    self.bmodel = models.SmallStrainRateIndependentPlasticity(self.elastic, 
+        flow)
+
+    self.fn = interpolate.PolynomialInterpolate([0.1])
+    self.n = 2.1
+
+    self.dmodel = damage.WorkDamage(self.elastic, self.fn, self.n, log = True)
+
+    self.model = damage.NEMLScalarDamagedModel_sd(self.elastic, self.bmodel, 
+            self.dmodel)
+
+    self.stress = np.array([100,-50.0,300.0,-99,50.0,125.0]) * 0.75
+    self.T = 100.0
+
+    self.s_np1 = self.stress
+    self.s_n = np.array([-25,150,250,-25,-100,25]) * 0.0
+
+    self.d_np1 = 0.2
+    self.d_n = 0.1
+
+    self.e_np1 = np.array([0.1,-0.01,0.15,-0.05,-0.1,0.15]) * 0.75
+    self.e_n = np.array([-0.05,0.025,-0.1,0.2,0.11,0.13]) * 0.0
+
+    self.T_np1 = self.T
+    self.T_n = 90.0
+
+    self.t_np1 = 5000.0
+    self.t_n = 0.9
+
+    self.u_n = 0.0
+    self.p_n = 0.0
+
+    # This is a rather boring baseline history state to probe, but I can't
+    # think of a better way to get a "generic" history from a generic model
+    self.hist_n = np.array([self.d_n] + list(self.bmodel.init_store()))
+    self.x_trial = np.array([50,-25,150,-150,190,100.0] + [0.41])
+
+    self.nsteps = 10
+    self.etarget = np.array([0.1,-0.025,0.02,0.015,-0.02,-0.05])
+    self.ttarget = 2.0
+
+    self.ee = np.dot(self.elastic.S(self.T_np1), self.s_np1*(1.0-self.d_np1) - self.s_n*(1.0-self.d_n))
+    self.de = self.e_np1 - self.e_n
+    self.dp = self.de - self.ee
+    self.dt = self.t_np1 - self.t_n
+    self.Wdot = np.dot(self.s_np1*(1.0-self.d_np1), self.dp) / self.dt
+
+  def test_definition(self):
+    damage = self.dmodel.damage(self.d_np1, self.d_n, self.e_np1, self.e_n,
+        self.s_np1, self.s_n, self.T_np1, self.T_n, self.t_np1, self.t_n)
+    should = self.d_n + self.n * self.d_np1**((self.n-1)/self.n) * self.Wdot * self.dt / np.exp(self.fn.value(np.log(self.Wdot)))
+
+    self.assertAlmostEqual(damage, should)
+
 class TestClassicalDamage(unittest.TestCase, CommonScalarDamageModel,
     CommonDamagedModel, CommonScalarDamageRate):
   def setUp(self):