Fixes

Argonne-National-Laboratory · Sep 20, 2023 · 5cc148d · 5cc148d
1 parent 80c78c0
commit 5cc148d
Show file tree

Hide file tree

Showing 2 changed files with 16 additions and 9 deletions.
diff --git a/src/damage.cxx b/src/damage.cxx
@@ -1378,14 +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 = val_n*std::pow(fabs(d_np1), (val_n-1.0)/val_n) * 
       dt / val * (1.0 - wrate / val * deriv / work_scale_) * x;
+
+  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;
+      wrate * dt / val / work_scale_ + other + n_partial;
 }
 
 void WorkDamage::ddamage_de(double d_np1, double d_n,

diff --git a/test/test_damage.py b/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, rtol = 1.0e-3))
+    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,7 +155,7 @@ 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-2))
 
 class CommonDamagedModel(object):
@@ -292,7 +292,7 @@ def setUp(self):
         flow)
 
     self.fn = interpolate.PolynomialInterpolate([0.1])
-    self.n = interpolate.PolynomialInterpolate([0.5, 1, 2.1])
+    self.n = interpolate.PolynomialInterpolate([0.1, 2.1])
 
     self.dmodel = damage.WorkDamage(self.elastic, self.fn, self.n, log = True)
 
@@ -338,7 +338,7 @@ 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.value(self.Wdot) * self.d_np1**((self.n.value(self.Wdot)-1)/self.n.value(self.Wdot)) * 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)