From 5cc148da109a95e06957c35ded7b14516896df0d Mon Sep 17 00:00:00 2001 From: Mark Messner Date: Wed, 20 Sep 2023 11:10:09 -0500 Subject: [PATCH] Fixes --- src/damage.cxx | 13 ++++++++++--- test/test_damage.py | 12 ++++++------ 2 files changed, 16 insertions(+), 9 deletions(-) diff --git a/src/damage.cxx b/src/damage.cxx index 3db22fe0..e9586a86 100644 --- 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 index 9cdffb75..ae0b4f60 100644 --- 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)