Interpolator for n for WorkDamage

include/damage.h

@@ -665,10 +665,16 @@ class NEML_EXPORT WorkDamage: public ScalarDamage {
 
   /// Get the derivative of the critical work, accounting for log scaling
   double dWcrit(double Wdot) const;
+
+  /// Get the damage exponent, accounting for log scaling if requested
+  double get_n(double Wdot) const;
+
+  /// Get the derivative of the damage exponent, accounting for log scaling
+  double get_dn(double Wdot) const;
 
  protected:
   std::shared_ptr<Interpolate> Wcrit_;
-  double n_;
+  std::shared_ptr<Interpolate> n_;
   double eps_;
   double work_scale_;
   bool log_;

src/damage.cxx

@@ -1291,7 +1291,7 @@ double StandardScalarDamage::dep(
 WorkDamage::WorkDamage(ParameterSet & params) :
       ScalarDamage(params), 
       Wcrit_(params.get_object_parameter<Interpolate>("Wcrit")),
-      n_(params.get_parameter<double>("n")),
+      n_(params.get_object_parameter<Interpolate>("n")),
       eps_(params.get_parameter<double>("eps")),
       work_scale_(params.get_parameter<double>("work_scale")),
       log_(params.get_parameter<bool>("log"))
@@ -1309,7 +1309,7 @@ ParameterSet WorkDamage::parameters()
 
   pset.add_parameter<NEMLObject>("elastic");
   pset.add_parameter<NEMLObject>("Wcrit");
-  pset.add_parameter<double>("n");
+  pset.add_parameter<NEMLObject>("n");
   pset.add_optional_parameter<double>("eps", 1e-30);
   pset.add_optional_parameter<double>("work_scale", 1.0);
   pset.add_optional_parameter<bool>("log", false);
@@ -1343,10 +1343,11 @@ void WorkDamage::damage(double d_np1, double d_n,
 
   double dt = t_np1 - t_n;
 
-  double val = Wcrit(wrate / work_scale_);
+  double val_Wc = Wcrit(wrate / work_scale_);
+  double val_n = get_n(wrate / work_scale_);
 
-  *dd = d_n + n_ * std::pow(std::fabs(d_np1), (n_-1.0)/n_) *
-        wrate * dt / val / work_scale_;
+  *dd = d_n + val_n * std::pow(std::fabs(d_np1), (val_n-1.0)/val_n) *
+        wrate * dt / val_Wc / work_scale_;
 }
 
 void WorkDamage::ddamage_dd(double d_np1, double d_n,
@@ -1377,13 +1378,21 @@ void WorkDamage::ddamage_dd(double d_np1, double d_n,
   double e[6];
   mat_vec(S, 6, s_np1, 6, e);
   double f = dot_vec(s_np1, e, 6);
-  double x = -wrate / (1.0 - std::fabs(d_np1)) / work_scale_ + 
-      (1.0 - std::fabs(d_np1)) * f / dt;
+  double ddf = -wrate / (1.0 - std::fabs(d_np1)) / work_scale_;
+  double x = ddf + (1.0 - std::fabs(d_np1)) * f / dt;
+
+  double val_n = get_n(wrate / work_scale_);
+  double dval_n = get_dn(wrate / work_scale_);
+  double val_Wc = Wcrit(wrate / work_scale_);
 
-  double other = n_*std::pow(fabs(d_np1), (n_-1.0)/n_) * 
+  double other = val_n*std::pow(fabs(d_np1), (val_n-1.0)/val_n) * 
       dt / val * (1.0 - wrate / val * deriv / work_scale_) * x;
-  *dd = (n_-1.0)*std::pow(std::fabs(d_np1), -1.0/n_) *
-      wrate * dt / val / work_scale_ + other;
+
+  double n_partial = (dt * wrate * std::pow(std::fabs(d_np1), (val_n - 1.0) / val_n) * (val_n + std::log(fabs(d_np1)))) 
+      / (val_n * work_scale_ * val_Wc) * dval_n * ddf;
+
+  *dd = (val_n-1.0)*std::pow(std::fabs(d_np1), -1.0/val_n) *
+      wrate * dt / val / work_scale_ + other + n_partial;
 }
 
 void WorkDamage::ddamage_de(double d_np1, double d_n,
@@ -1393,6 +1402,8 @@ void WorkDamage::ddamage_de(double d_np1, double d_n,
                    double t_np1, double t_n,
                    double * const dd) const
 {
+
+  // Get work rate
   double wrate = workrate(e_np1, e_n, s_np1, s_n, T_np1, T_n, t_np1, t_n,
                           std::fabs(d_np1), d_n);
 
@@ -1403,14 +1414,26 @@ void WorkDamage::ddamage_de(double d_np1, double d_n,
     return;
   }
 
-  double val = Wcrit(wrate / work_scale_);
-  double dval = dWcrit(wrate / work_scale_);
-
-  double fact = n_ * std::pow(std::fabs(d_np1), (n_-1.0)/n_) / val * 
-      (1.0 - wrate / val * dval / work_scale_) * (1.0 - std::fabs(d_np1));
-
+  // Initialise interpolation variables
+  double val_Wc  = Wcrit(wrate / work_scale_);
+  double val_dWc = dWcrit(wrate / work_scale_);
+  double val_n   = get_n(wrate / work_scale_);
+  double val_dn  = get_dn(wrate / work_scale_);
+
+  // Calculate partial derivative
+  //   dR/de = (dwrate/ds) * ( (dR/dwrate) + (dR/dWc)*(dWc/dwrate)
+  //   ... (dR/dn)*(dn/dwrate) )
+  double d_np1_pow = std::pow(fabs(d_np1), (val_n-1.0)/val_n);
+  double dR_dwrate = val_n * d_np1_pow / val_Wc;
+  double dR_dWcrit = -val_n * d_np1_pow * (wrate / work_scale_) /
+    std::pow(val_Wc, 2);
+  double dR_dn     = (wrate / work_scale_) / val_Wc * d_np1_pow *
+    (val_n + std::log(fabs(d_np1))) / val_n;
+  double factor    = dR_dwrate + dR_dWcrit*val_dWc + dR_dn*val_dn;
+
+  // Multiply everything by dwrate/ds
   for (size_t i = 0; i < 6; i++) {
-    dd[i] = fact * s_np1[i];
+    dd[i] = factor * s_np1[i] * (1.0 - std::fabs(d_np1));
   }
 }
 
@@ -1421,38 +1444,45 @@ void WorkDamage::ddamage_ds(double d_np1, double d_n,
                    double t_np1, double t_n,
                    double * const dd) const
 {
+  // Get work rate and check
   double wrate = workrate(e_np1, e_n, s_np1, s_n, T_np1, T_n, t_np1, t_n,
-                          std::fabs(d_np1), d_n);
-
+    std::fabs(d_np1), d_n);
   if ((d_np1 <= 0.0) || (wrate == 0.0)) {
     std::fill(dd, dd+6, 0.0);
     return;
   }
 
-  double val = Wcrit(wrate / work_scale_);
-  double dval = dWcrit(wrate / work_scale_);
-
+  // Initialise interpolation variables
+  double val_Wc  = Wcrit(wrate / work_scale_);
+  double val_dWc = dWcrit(wrate / work_scale_);
+  double val_n   = get_n(wrate / work_scale_);
+  double val_dn  = get_dn(wrate / work_scale_);
+
+  // Calculate partial derivative
+  //   dR/de = (dwrate/ds) * ( (dR/dwrate) + (dR/dWc)*(dWc/dwrate) +
+  //   ... (dR/dn)*(dn/dwrate) )
+  double d_np1_pow = std::pow(fabs(d_np1), (val_n-1.0)/val_n);
+  double dR_dwrate = val_n * d_np1_pow / val_Wc;
+  double dR_dWcrit = -val_n * d_np1_pow * (wrate / work_scale_) /
+    std::pow(val_Wc, 2);
+  double dR_dn     = (wrate / work_scale_) / val_Wc * d_np1_pow *
+    (val_n + std::log(fabs(d_np1))) / val_n;
+  double factor    = dR_dwrate + dR_dWcrit*val_dWc + dR_dn*val_dn;
+
+  // Get inverse elastic tensor
   double S[36];
   elastic_->S(T_np1, S);
-
-  double ds[6];
-  double de[6];
-  for (int i=0; i<6; i++) {
-    ds[i] = s_np1[i] * (1.0-d_np1) - s_n[i] * (1.0-d_n);
-    de[i] = e_np1[i] - e_n[i];
-  }
 
-  double dee[6];
-  mat_vec(S, 6, ds, 6, dee);
-
-  double e[6];
-  mat_vec(S, 6, s_np1, 6, e);
-
-  double fact = n_ * std::pow(fabs(d_np1), (n_-1.0)/n_) / val * 
-      (1.0 - wrate / val * dval / work_scale_)*(1.0-std::fabs(d_np1));
+  // Calculate S:s_n and S:s_np1
+  double S_s_n[6];
+  mat_vec(S, 6, s_n, 6, S_s_n);
+  double S_s_np1[6];
+  mat_vec(S, 6, s_np1, 6, S_s_np1);
 
+  // Multiply everything by dwrate/de
   for (size_t i = 0; i < 6; i++) {
-    dd[i] = fact * (de[i] - dee[i] - (1.0-std::fabs(d_np1))*e[i]);
+    dd[i] = factor * (e_np1[i] - e_n[i] + (S_s_n[i] - 2*S_s_np1[i]) *
+      (1.0 - std::fabs(d_np1))) * (1.0 - std::fabs(d_np1));
   }
 }
 
@@ -1492,7 +1522,6 @@ double WorkDamage::Wcrit(double Wdot) const
 {
   if (log_)
     return std::pow(10.0, Wcrit_->value(std::log10(Wdot)));
-
   return Wcrit_->value(Wdot);
 }
 
@@ -1501,10 +1530,23 @@ 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);
 }
 
+double WorkDamage::get_n(double Wdot) const
+{
+  if (log_)
+    return n_->value(std::log10(Wdot));
+  return n_->value(Wdot);
+}
+
+double WorkDamage::get_dn(double Wdot) const
+{
+  if (log_)
+    return n_->derivative(std::log10(Wdot)) / std::log(10) / Wdot;
+  return n_->derivative(Wdot);
+}
+
 PowerLawDamage::PowerLawDamage(ParameterSet & params) :
       StandardScalarDamage(params), 
       A_(params.get_object_parameter<Interpolate>("A")),
@@ -1575,8 +1617,7 @@ double PowerLawDamage::se(const double * const s) const
                                               pow(s[5], 2.0))) / 2.0);
 }
 
-ExponentialWorkDamage::ExponentialWorkDamage(ParameterSet
-                                                                       & params) :
+ExponentialWorkDamage::ExponentialWorkDamage(ParameterSet & params) :
       StandardScalarDamage(params), 
       W0_(params.get_object_parameter<Interpolate>("W0")),
       k0_(params.get_object_parameter<Interpolate>("k0")),

test/test_damage.py

@@ -88,7 +88,7 @@ def test_ddamage_ddamage(self):
         self.s_np1, self.s_n, self.T_np1, self.T_n, self.t_np1, self.t_n)
     dd_calcd = differentiate(dfn, self.d_np1)
 
-    self.assertTrue(np.isclose(dd_model, dd_calcd))
+    self.assertTrue(np.isclose(dd_model, dd_calcd, rtol = 1.0e-4))
 
   def test_ddamage_dstrain(self):
     dd_model = self.dmodel.ddamage_de(self.d_np1, self.d_n, self.e_np1, self.e_n,
@@ -97,7 +97,7 @@ def test_ddamage_dstrain(self):
         self.s_np1, self.s_n, self.T_np1, self.T_n, self.t_np1, self.t_n)
     dd_calcd = differentiate(dfn, self.e_np1)[0]
 
-    self.assertTrue(np.allclose(dd_model, dd_calcd, rtol = 1.0e-3))
+    self.assertTrue(np.allclose(dd_model, dd_calcd, rtol = 1.0e-4))
 
   def test_ddamage_dstress(self):
     dd_model = self.dmodel.ddamage_ds(self.d_np1, self.d_n, self.e_np1, self.e_n,
@@ -106,7 +106,7 @@ def test_ddamage_dstress(self):
         s, self.s_n, self.T_np1, self.T_n, self.t_np1, self.t_n)
     dd_calcd = differentiate(dfn, self.s_np1)[0]
 
-    self.assertTrue(np.allclose(dd_model, dd_calcd, rtol = 1.0e-3))
+    self.assertTrue(np.allclose(dd_model, dd_calcd, rtol = 1.0e-4))
 
   def test_nparams(self):
     self.assertEqual(self.model.nparams, 7)
@@ -155,8 +155,8 @@ def test_jacobian(self):
     R, J = self.model.RJ(self.x_trial, trial_state)
     Jnum = differentiate(lambda x: self.model.RJ(x, trial_state)[0], 
         self.x_trial)
-    
-    self.assertTrue(np.allclose(J, Jnum, rtol = 1.0e-3))
+
+    self.assertTrue(np.allclose(J, Jnum, rtol = 1.0e-2))
 
 class CommonDamagedModel(object):
   def test_nstore(self):
@@ -186,8 +186,6 @@ def test_tangent_proportional_strain(self):
       A_num = differentiate(lambda e: self.model.update_sd(e, e_n,
         self.T, self.T, t_np1, t_n, s_n, hist_n, u_n, p_n)[0], e_np1)
 
-      print(A_num)
-      print(A_np1)
       self.assertTrue(np.allclose(A_num, A_np1, rtol = 5.0e-2, atol = 1.0e-1))
 
       e_n = np.copy(e_np1)
@@ -220,8 +218,8 @@ def setUp(self):
     self.bmodel = models.SmallStrainRateIndependentPlasticity(self.elastic, 
         flow)
 
-    self.fn = interpolate.PolynomialInterpolate([0.1,5.0, 1e-8])
-    self.n = 2.1
+    self.fn = interpolate.PolynomialInterpolate([0.1, 5.0, 1e-8])
+    self.n = interpolate.PolynomialInterpolate([0.5, 1, 2.1])
 
     self.dmodel = damage.WorkDamage(self.elastic, self.fn, self.n)
 
@@ -263,12 +261,11 @@ def setUp(self):
     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 / self.fn.value(self.Wdot) 
-
+    should = self.d_n + self.n.value(self.Wdot) * self.d_np1**((self.n.value(self.Wdot)-1)/self.n.value(self.Wdot)) * self.Wdot * self.dt / self.fn.value(self.Wdot) 
     self.assertAlmostEqual(damage, should)
 
 class TestWorkDamageLog(unittest.TestCase, CommonScalarDamageModel,
@@ -295,7 +292,7 @@ def setUp(self):
         flow)
 
     self.fn = interpolate.PolynomialInterpolate([0.1])
-    self.n = 2.1
+    self.n = interpolate.PolynomialInterpolate([0.1, 2.1])
 
     self.dmodel = damage.WorkDamage(self.elastic, self.fn, self.n, log = True)
 
@@ -341,13 +338,14 @@ def setUp(self):
   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 / 10.0**(self.fn.value(np.log10(self.Wdot)))
+    should = self.d_n + self.n.value(np.log10(self.Wdot)) * self.d_np1**((self.n.value(np.log10(self.Wdot))-1)/self.n.value(np.log10(self.Wdot))) * self.Wdot * self.dt / 10.0**(self.fn.value(np.log10(self.Wdot)))
 
     self.assertAlmostEqual(damage, should)
 
 class TestClassicalDamage(unittest.TestCase, CommonScalarDamageModel,
     CommonDamagedModel, CommonScalarDamageRate):
   def setUp(self):
+
     self.E = 92000.0
     self.nu = 0.3