From 99c5476e0bee1fd8057e1875d1bb92433eab96f1 Mon Sep 17 00:00:00 2001 From: ralberd Date: Mon, 27 Nov 2023 15:15:42 -0700 Subject: [PATCH 1/9] refactoring to remove the extra mtk_log_sqrt call in energy computation - is not matching previous implementation --- optimism/material/HyperViscoelastic.py | 100 ++++++---------------- optimism/material/test/test_HyperVisco.py | 81 ++++++++++++++++++ 2 files changed, 106 insertions(+), 75 deletions(-) create mode 100644 optimism/material/test/test_HyperVisco.py diff --git a/optimism/material/HyperViscoelastic.py b/optimism/material/HyperViscoelastic.py index caca5a61..5ac28be0 100644 --- a/optimism/material/HyperViscoelastic.py +++ b/optimism/material/HyperViscoelastic.py @@ -109,18 +109,21 @@ def _eq_strain_energy(dispGrad, props): return Wdev + Wvol def _neq_strain_energy(dispGrad, stateOld, dt, props): - I = np.identity(3) - F = dispGrad + I - state_new = _compute_state_new(dispGrad, stateOld, dt, props) - # Fvs = state_new.reshape((NUM_PRONY_TERMS, 3, 3)) - Fv_new = state_new.reshape((3, 3)) - G_neq = props[PROPS_G_neq] tau = props[PROPS_TAU] eta = G_neq * tau - Fe = F @ np.linalg.inv(Fv_new) - Ee = 0.5 * TensorMath.mtk_log_sqrt(Fe.T @ Fe) + Ee_trial = _compute_elastic_logarithmic_strain(dispGrad, stateOld) + state_inc = _compute_state_increment(Ee_trial, dt, props) + + Ee = Ee_trial - state_inc # doesn't work + # # The below works - matches old implementation + # F = dispGrad + np.identity(3) + # Fv_old = stateOld.reshape((3, 3)) + # Fv_new = linalg.expm(state_inc)@Fv_old + # Fe = F @ np.linalg.inv(Fv_new) + # Ee = 0.5 * TensorMath.mtk_log_sqrt(Fe.T @ Fe) + Me = 2. * G_neq * Ee M_bar = TensorMath.norm_of_deviator_squared(Me) gamma_dot = M_bar / eta @@ -128,80 +131,23 @@ def _neq_strain_energy(dispGrad, stateOld, dt, props): visco_energy = 0.5 * dt * eta * gamma_dot**2 W_neq = G_neq * TensorMath.norm_of_deviator_squared(Ee) + visco_energy - # def vmap_body(n, Fv): - # G_neq = props[PROPS_G_neq + 2 * n] - # # tau = props[PROPS_TAU + 2 * n] - - # Fe = F @ np.linalg.inv(Fv) - # Ee = 0.5 * TensorMath.mtk_log_sqrt(Fe.T @ Fe) - - # # viscous shearing - # # Me = 2. * G_neq * Ee - # # M_bar = TensorMath.norm_of_deviator_squared(Me) - # # visco_energy = (dt / (G_neq * tau)) * M_bar**2 - - # # still need another term I think - # W_neq = G_neq * TensorMath.norm_of_deviator_squared(Ee) #+ visco_energy - # return W_neq - - # W_neq = np.sum(vmap(vmap_body, in_axes=(0, 0))(np.arange(NUM_PRONY_TERMS), Fvs)) return W_neq -# state update def _compute_state_new(dispGrad, stateOld, dt, props): - state_inc = _compute_state_increment(dispGrad, stateOld, dt, props) - # Fv_olds = stateOld.reshape((NUM_PRONY_TERMS, 3, 3)) - # Fv_incs = state_inc.reshape((NUM_PRONY_TERMS, 3, 3)) - - # def vmap_body(n, Fv_old, Fv_inc): - # Fv_new = Fv_inc @ Fv_old - # return Fv_new.ravel() - - # state_new = np.hstack(vmap(vmap_body, in_axes=(0, 0, 0))(np.arange(NUM_PRONY_TERMS), Fv_olds, Fv_incs)) + Ee_trial = _compute_elastic_logarithmic_strain(dispGrad, stateOld) + state_inc = _compute_state_increment(Ee_trial, dt, props) Fv_old = stateOld.reshape((3, 3)) - Fv_inc = state_inc.reshape((3, 3)) - state_new = (Fv_inc @ Fv_old).ravel() - return state_new - -def _compute_state_increment(dispGrad, stateOld, dt, props): - I = np.identity(3) - F = dispGrad + I - # Fv_olds = stateOld.reshape((NUM_PRONY_TERMS, 3, 3)) - - # def vmap_body(n, Fv_old): - # # TODO add shift factor - # G_neq = props[PROPS_G_neq + 2 * n] - # tau = props[PROPS_TAU + 2 * n] - - # # kinematics - # Fe_trial = F @ np.linalg.inv(Fv_old) - # Ee_trial = 0.5 * TensorMath.mtk_log_sqrt(Fe_trial.T @ Fe_trial) - # Ee_dev = Ee_trial - (1. / 3.) * np.trace(Ee_trial) * I + Fv_new = linalg.expm(state_inc)@Fv_old + return Fv_new.ravel() - # # updates - # integration_factor = 1. / (1. + dt / tau) - - # Me = 2.0 * G_neq * Ee_dev - # Me = integration_factor * Me - - # Dv = (1. / (2. * G_neq * tau)) * Me - # A = dt * Dv - - # Fv_inc = linalg.expm(A) - - # return Fv_inc.ravel() - - # state_inc = np.hstack(vmap(vmap_body, in_axes=(0, 0))(np.arange(NUM_PRONY_TERMS), Fv_olds)) - - Fv_old = stateOld.reshape((3, 3)) +def _compute_state_increment(elasticStrain, dt, props): G_neq = props[PROPS_G_neq] tau = props[PROPS_TAU] - Fe_trial = F @ np.linalg.inv(Fv_old) - Ee_trial = 0.5 * TensorMath.mtk_log_sqrt(Fe_trial.T @ Fe_trial) - Ee_dev = Ee_trial - (1. / 3.) * np.trace(Ee_trial) * I + # Ee_dev = elasticStrain - (1. / 3.) * np.trace(elasticStrain) * np.identity(3) + Ee_dev = TensorMath.compute_deviatoric_tensor(elasticStrain) integration_factor = 1. / (1. + dt / tau) @@ -209,8 +155,12 @@ def _compute_state_increment(dispGrad, stateOld, dt, props): Me = integration_factor * Me Dv = (1. / (2. * G_neq * tau)) * Me - A = dt * Dv + return dt * Dv - Fv_inc = linalg.expm(A) +def _compute_elastic_logarithmic_strain(dispGrad, stateOld): + F = dispGrad + np.identity(3) + Fv_old = stateOld.reshape((3, 3)) + + Fe_trial = F @ np.linalg.inv(Fv_old) - return Fv_inc.ravel() \ No newline at end of file + return 0.5 * TensorMath.mtk_log_sqrt(Fe_trial.T @ Fe_trial) \ No newline at end of file diff --git a/optimism/material/test/test_HyperVisco.py b/optimism/material/test/test_HyperVisco.py new file mode 100644 index 00000000..0ec47b70 --- /dev/null +++ b/optimism/material/test/test_HyperVisco.py @@ -0,0 +1,81 @@ +import unittest + +import jax +import jax.numpy as np +from jax.scipy import linalg + +from optimism.material import HyperViscoelastic as HyperVisco +from optimism.test.TestFixture import TestFixture + +def make_disp_grad_from_strain(strain): + return linalg.expm(strain) - np.identity(3) + + +class GradOfPlasticityModelFixture(TestFixture): + def setUp(self): + + G_eq = 0.855 # MPa + K_eq = 1000*G_eq # MPa + G_neq_1 = 5.0 + tau_1 = 0.1 + self.props = { + 'equilibrium bulk modulus' : K_eq, + 'equilibrium shear modulus' : G_eq, + 'non equilibrium shear modulus': G_neq_1, + 'relaxation time' : tau_1, + } + + materialModel = HyperVisco.create_material_model_functions(self.props) + + self.energy_density = jax.jit(materialModel.compute_energy_density) + self.compute_state_new = materialModel.compute_state_new + self.compute_initial_state = materialModel.compute_initial_state + + def test_zero_point(self): + dispGrad = np.zeros((3,3)) + initialState = self.compute_initial_state() + dt = 1.0 + + energy = self.energy_density(dispGrad, initialState, dt) + self.assertNear(energy, 0.0, 12) + + state = self.compute_state_new(dispGrad, initialState, dt) + self.assertArrayNear(state, np.eye(3).ravel(), 12) + + def test_regression_nonzero_point(self): + key = jax.random.PRNGKey(1) + dispGrad = jax.random.uniform(key, (3, 3)) + initialState = self.compute_initial_state() + dt = 1.0 + + energy = self.energy_density(dispGrad, initialState, dt) + self.assertNear(energy, 193.6029283822052, 12) + + state = self.compute_state_new(dispGrad, initialState, dt) + stateGold = np.array([0.946976197737, 0.206631668313, 0.220846029827, + 0.206631668313, 1.155826112068, 0.015472488047, + 0.220846029827, 0.015472488047, 1.003179626352]) + self.assertArrayNear(state, stateGold, 12) + + + # def test_elastic_energy(self): + # strainBelowYield = 0.5*self.props['yield strength']/self.props['elastic modulus'] + + # strain = strainBelowYield*np.diag(np.array([1.0, -self.props['poisson ratio'], -self.props['poisson ratio']])) + # dispGrad = make_disp_grad_from_strain(strain) + # dt = 1.0 + + # state = self.compute_initial_state() + + # energy = self.energy_density(dispGrad, state, dt) + # WExact = 0.5*self.props['elastic modulus']*strainBelowYield**2 + # self.assertNear(energy, WExact, 12) + + # F = dispGrad + np.identity(3) + # kirchhoffStress = self.stress_func(dispGrad, state, dt) @ F.T + # kirchhoffstressExact = np.zeros((3,3)).at[0,0].set(self.props['elastic modulus']*strainBelowYield) + # self.assertArrayNear(kirchhoffStress, kirchhoffstressExact, 12) + + +if __name__ == '__main__': + unittest.main() From 175132a24dc7740dd5eb635b4bee84007183a562 Mon Sep 17 00:00:00 2001 From: ralberd Date: Tue, 28 Nov 2023 11:49:30 -0700 Subject: [PATCH 2/9] fixing bug in computation of elastic log strain, cleaning up some other comments and stuff --- optimism/material/HyperViscoelastic.py | 41 ++++---------------------- 1 file changed, 5 insertions(+), 36 deletions(-) diff --git a/optimism/material/HyperViscoelastic.py b/optimism/material/HyperViscoelastic.py index 5ac28be0..6a2489ff 100644 --- a/optimism/material/HyperViscoelastic.py +++ b/optimism/material/HyperViscoelastic.py @@ -78,18 +78,6 @@ def _make_properties(properties): properties['relaxation time'] ]) - # props = np.hstack((props, properties['number of prony terms'])) - - # for n in range(1, properties['number of prony terms'] + 1): - # print('Prony branch %s properties' % n) - # print(' Shear modulus = %s' % properties['non equilibrium shear modulus %s' % n]) - # print(' Relaxation time = %s' % properties['relaxation time %s' % n]) - # props = np.hstack( - # (props, np.array([properties['non equilibrium shear modulus %s' % n], - # properties['relaxation time %s' % n]]))) - - - return props def _energy_density(dispGrad, state, dt, props): @@ -97,7 +85,6 @@ def _energy_density(dispGrad, state, dt, props): W_neq = _neq_strain_energy(dispGrad, state, dt, props) return W_eq + W_neq -# TODO generalize to arbitrary strain energy density def _eq_strain_energy(dispGrad, props): K, G = props[PROPS_K_eq], props[PROPS_G_eq] F = dispGrad + np.eye(3) @@ -115,14 +102,7 @@ def _neq_strain_energy(dispGrad, stateOld, dt, props): Ee_trial = _compute_elastic_logarithmic_strain(dispGrad, stateOld) state_inc = _compute_state_increment(Ee_trial, dt, props) - - Ee = Ee_trial - state_inc # doesn't work - # # The below works - matches old implementation - # F = dispGrad + np.identity(3) - # Fv_old = stateOld.reshape((3, 3)) - # Fv_new = linalg.expm(state_inc)@Fv_old - # Fe = F @ np.linalg.inv(Fv_new) - # Ee = 0.5 * TensorMath.mtk_log_sqrt(Fe.T @ Fe) + Ee = Ee_trial - state_inc Me = 2. * G_neq * Ee M_bar = TensorMath.norm_of_deviator_squared(Me) @@ -130,9 +110,7 @@ def _neq_strain_energy(dispGrad, stateOld, dt, props): # visco_energy = (dt / (G_neq * tau)) * M_bar**2 visco_energy = 0.5 * dt * eta * gamma_dot**2 - W_neq = G_neq * TensorMath.norm_of_deviator_squared(Ee) + visco_energy - - return W_neq + return G_neq * TensorMath.norm_of_deviator_squared(Ee) + visco_energy def _compute_state_new(dispGrad, stateOld, dt, props): Ee_trial = _compute_elastic_logarithmic_strain(dispGrad, stateOld) @@ -143,24 +121,15 @@ def _compute_state_new(dispGrad, stateOld, dt, props): return Fv_new.ravel() def _compute_state_increment(elasticStrain, dt, props): - G_neq = props[PROPS_G_neq] tau = props[PROPS_TAU] - - # Ee_dev = elasticStrain - (1. / 3.) * np.trace(elasticStrain) * np.identity(3) - Ee_dev = TensorMath.compute_deviatoric_tensor(elasticStrain) - integration_factor = 1. / (1. + dt / tau) - Me = 2.0 * G_neq * Ee_dev - Me = integration_factor * Me - - Dv = (1. / (2. * G_neq * tau)) * Me - return dt * Dv + Ee_dev = TensorMath.compute_deviatoric_tensor(elasticStrain) + return dt * integration_factor * Ee_dev / tau def _compute_elastic_logarithmic_strain(dispGrad, stateOld): F = dispGrad + np.identity(3) Fv_old = stateOld.reshape((3, 3)) Fe_trial = F @ np.linalg.inv(Fv_old) - - return 0.5 * TensorMath.mtk_log_sqrt(Fe_trial.T @ Fe_trial) \ No newline at end of file + return TensorMath.mtk_log_sqrt(Fe_trial.T @ Fe_trial) \ No newline at end of file From 9fbd8d140a8bb4a7736d7733cb84cd634872a81a Mon Sep 17 00:00:00 2001 From: ralberd Date: Tue, 28 Nov 2023 11:53:42 -0700 Subject: [PATCH 3/9] removing comments and some other unused stuff --- optimism/material/HyperViscoelastic.py | 30 -------------------------- 1 file changed, 30 deletions(-) diff --git a/optimism/material/HyperViscoelastic.py b/optimism/material/HyperViscoelastic.py index 6a2489ff..7a03eee2 100644 --- a/optimism/material/HyperViscoelastic.py +++ b/optimism/material/HyperViscoelastic.py @@ -1,11 +1,8 @@ import jax.numpy as np -from jax import vmap from jax.scipy import linalg from optimism import TensorMath from optimism.material.MaterialModel import MaterialModel - -# props PROPS_K_eq = 0 PROPS_G_eq = 1 PROPS_G_neq = 2 @@ -13,31 +10,20 @@ NUM_PRONY_TERMS = -1 -# isvs VISCOUS_DISTORTION = slice(0, 9) def create_material_model_functions(properties): - # prop processing density = properties.get('density') props = _make_properties(properties) - # energy function wrapper def energy_density(dispGrad, state, dt): return _energy_density(dispGrad, state, dt, props) - # wrapper for state var ics def compute_initial_state(shape=(1,)): - # num_prony_terms = properties['number of prony terms'] - - # def vmap_body(_): - # return np.identity(3).ravel() - - # state = np.hstack(vmap(vmap_body, in_axes=(0,))(np.arange(num_prony_terms))) state = np.identity(3).ravel() return state - # update state vars wrapper def compute_state_new(dispGrad, state, dt): state = _compute_state_new(dispGrad, state, dt, props) return state @@ -49,28 +35,12 @@ def compute_state_new(dispGrad, state, dt): density ) -# implementation def _make_properties(properties): - # assert properties['number of prony terms'] > 0, 'Need at least 1 prony term' assert 'equilibrium bulk modulus' in properties.keys() assert 'equilibrium shear modulus' in properties.keys() - # for n in range(1, properties['number of prony terms'] + 1): - # assert 'non equilibrium shear modulus %s' % n in properties.keys() - # assert 'relaxation time %s' % n in properties.keys() assert 'non equilibrium shear modulus' in properties.keys() assert 'relaxation time' in properties.keys() - print('Equilibrium properties') - print(' Bulk modulus = %s' % properties['equilibrium bulk modulus']) - print(' Shear modulus = %s' % properties['equilibrium shear modulus']) - print('Prony branch properties') - print(' Shear modulus = %s' % properties['non equilibrium shear modulus']) - print(' Relaxation time = %s' % properties['relaxation time']) - # this is dirty, fuck jax (can't use an int from a jax numpy array or else jit tries to trace that) - # global NUM_PRONY_TERMS - # NUM_PRONY_TERMS = properties['number of prony terms'] - - # first pack equilibrium properties props = np.array([ properties['equilibrium bulk modulus'], properties['equilibrium shear modulus'], From 7e7b20d42f9209b1d0c61a1e32d2f0c6714c1476 Mon Sep 17 00:00:00 2001 From: ralberd Date: Tue, 28 Nov 2023 11:55:34 -0700 Subject: [PATCH 4/9] adding test to check that update of elastic log strain matches previous way of updating --- optimism/material/test/test_HyperVisco.py | 40 ++++++++++++++--------- 1 file changed, 25 insertions(+), 15 deletions(-) diff --git a/optimism/material/test/test_HyperVisco.py b/optimism/material/test/test_HyperVisco.py index 0ec47b70..13214bb3 100644 --- a/optimism/material/test/test_HyperVisco.py +++ b/optimism/material/test/test_HyperVisco.py @@ -7,6 +7,8 @@ from optimism.material import HyperViscoelastic as HyperVisco from optimism.test.TestFixture import TestFixture +from optimism import TensorMath + def make_disp_grad_from_strain(strain): return linalg.expm(strain) - np.identity(3) @@ -57,24 +59,32 @@ def test_regression_nonzero_point(self): 0.220846029827, 0.015472488047, 1.003179626352]) self.assertArrayNear(state, stateGold, 12) + def test_log_strain_updates(self): + key = jax.random.PRNGKey(1) + dispGrad = jax.random.uniform(key, (3, 3)) + initialState = self.compute_initial_state() + dt = 1.0 - # def test_elastic_energy(self): - # strainBelowYield = 0.5*self.props['yield strength']/self.props['elastic modulus'] - - # strain = strainBelowYield*np.diag(np.array([1.0, -self.props['poisson ratio'], -self.props['poisson ratio']])) - # dispGrad = make_disp_grad_from_strain(strain) - # dt = 1.0 - - # state = self.compute_initial_state() + props = np.array([ + self.props['equilibrium bulk modulus'], + self.props['equilibrium shear modulus'], + self.props['non equilibrium shear modulus'], + self.props['relaxation time'] + ]) + + Ee_trial = HyperVisco._compute_elastic_logarithmic_strain(dispGrad, initialState) + state_inc = HyperVisco._compute_state_increment(Ee_trial, dt, props) + + Ee_new_way = Ee_trial - state_inc # doesn't work - # energy = self.energy_density(dispGrad, state, dt) - # WExact = 0.5*self.props['elastic modulus']*strainBelowYield**2 - # self.assertNear(energy, WExact, 12) + # old implementation + F = dispGrad + np.identity(3) + Fv_old = initialState.reshape((3, 3)) + Fv_new = linalg.expm(state_inc)@Fv_old + Fe = F @ np.linalg.inv(Fv_new) + Ee_old_way = TensorMath.mtk_log_sqrt(Fe.T @ Fe) - # F = dispGrad + np.identity(3) - # kirchhoffStress = self.stress_func(dispGrad, state, dt) @ F.T - # kirchhoffstressExact = np.zeros((3,3)).at[0,0].set(self.props['elastic modulus']*strainBelowYield) - # self.assertArrayNear(kirchhoffStress, kirchhoffstressExact, 12) + self.assertArrayNear(Ee_new_way.ravel(), Ee_old_way.ravel(), 12) if __name__ == '__main__': From ae6a48f33c32c311900e73d846d3c45aaf0b07b5 Mon Sep 17 00:00:00 2001 From: ralberd Date: Tue, 28 Nov 2023 12:02:05 -0700 Subject: [PATCH 5/9] updating gold values now that log strain computation bug (extra multiplication by 0.5) was fixed --- optimism/material/test/test_HyperVisco.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/optimism/material/test/test_HyperVisco.py b/optimism/material/test/test_HyperVisco.py index 13214bb3..f7494a3b 100644 --- a/optimism/material/test/test_HyperVisco.py +++ b/optimism/material/test/test_HyperVisco.py @@ -51,12 +51,12 @@ def test_regression_nonzero_point(self): dt = 1.0 energy = self.energy_density(dispGrad, initialState, dt) - self.assertNear(energy, 193.6029283822052, 12) + self.assertNear(energy, 133.3469451269987, 12) state = self.compute_state_new(dispGrad, initialState, dt) - stateGold = np.array([0.946976197737, 0.206631668313, 0.220846029827, - 0.206631668313, 1.155826112068, 0.015472488047, - 0.220846029827, 0.015472488047, 1.003179626352]) + stateGold = np.array([0.988233534321, 0.437922586964, 0.433881277313, + 0.437922586964, 1.378870045574, 0.079038974065, + 0.433881277313, 0.079038974065, 1.055381729505]) self.assertArrayNear(state, stateGold, 12) def test_log_strain_updates(self): From 44f3fcd40d8a8a19388a775b68c801d298980b4b Mon Sep 17 00:00:00 2001 From: ralberd Date: Tue, 28 Nov 2023 15:40:04 -0700 Subject: [PATCH 6/9] comment to make strain increment return value expression clear --- optimism/material/HyperViscoelastic.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/optimism/material/HyperViscoelastic.py b/optimism/material/HyperViscoelastic.py index 7a03eee2..8e7535df 100644 --- a/optimism/material/HyperViscoelastic.py +++ b/optimism/material/HyperViscoelastic.py @@ -95,7 +95,7 @@ def _compute_state_increment(elasticStrain, dt, props): integration_factor = 1. / (1. + dt / tau) Ee_dev = TensorMath.compute_deviatoric_tensor(elasticStrain) - return dt * integration_factor * Ee_dev / tau + return dt * integration_factor * Ee_dev / tau # dt * D def _compute_elastic_logarithmic_strain(dispGrad, stateOld): F = dispGrad + np.identity(3) From b5b4f3196d253dc28301692e08724d158340e676 Mon Sep 17 00:00:00 2001 From: ralberd Date: Tue, 28 Nov 2023 15:47:01 -0700 Subject: [PATCH 7/9] incorporating Brandon's comments --- optimism/material/HyperViscoelastic.py | 19 +++++++++++-------- 1 file changed, 11 insertions(+), 8 deletions(-) diff --git a/optimism/material/HyperViscoelastic.py b/optimism/material/HyperViscoelastic.py index 8e7535df..38c6eea9 100644 --- a/optimism/material/HyperViscoelastic.py +++ b/optimism/material/HyperViscoelastic.py @@ -36,10 +36,13 @@ def compute_state_new(dispGrad, state, dt): ) def _make_properties(properties): - assert 'equilibrium bulk modulus' in properties.keys() - assert 'equilibrium shear modulus' in properties.keys() - assert 'non equilibrium shear modulus' in properties.keys() - assert 'relaxation time' in properties.keys() + + print('Equilibrium properties') + print(' Bulk modulus = %s' % properties['equilibrium bulk modulus']) + print(' Shear modulus = %s' % properties['equilibrium shear modulus']) + print('Prony branch properties') + print(' Shear modulus = %s' % properties['non equilibrium shear modulus']) + print(' Relaxation time = %s' % properties['relaxation time']) props = np.array([ properties['equilibrium bulk modulus'], @@ -71,8 +74,8 @@ def _neq_strain_energy(dispGrad, stateOld, dt, props): eta = G_neq * tau Ee_trial = _compute_elastic_logarithmic_strain(dispGrad, stateOld) - state_inc = _compute_state_increment(Ee_trial, dt, props) - Ee = Ee_trial - state_inc + delta_Ev = _compute_state_increment(Ee_trial, dt, props) + Ee = Ee_trial - delta_Ev Me = 2. * G_neq * Ee M_bar = TensorMath.norm_of_deviator_squared(Me) @@ -84,10 +87,10 @@ def _neq_strain_energy(dispGrad, stateOld, dt, props): def _compute_state_new(dispGrad, stateOld, dt, props): Ee_trial = _compute_elastic_logarithmic_strain(dispGrad, stateOld) - state_inc = _compute_state_increment(Ee_trial, dt, props) + delta_Ev = _compute_state_increment(Ee_trial, dt, props) Fv_old = stateOld.reshape((3, 3)) - Fv_new = linalg.expm(state_inc)@Fv_old + Fv_new = linalg.expm(delta_Ev)@Fv_old return Fv_new.ravel() def _compute_state_increment(elasticStrain, dt, props): From d6b3a154f19a58e4502aa6c5ab9b7ad608e9bc2e Mon Sep 17 00:00:00 2001 From: ralberd Date: Tue, 28 Nov 2023 16:08:20 -0700 Subject: [PATCH 8/9] removing test that was used for debugging --- optimism/material/test/test_HyperVisco.py | 29 ----------------------- 1 file changed, 29 deletions(-) diff --git a/optimism/material/test/test_HyperVisco.py b/optimism/material/test/test_HyperVisco.py index f7494a3b..f3d23646 100644 --- a/optimism/material/test/test_HyperVisco.py +++ b/optimism/material/test/test_HyperVisco.py @@ -7,8 +7,6 @@ from optimism.material import HyperViscoelastic as HyperVisco from optimism.test.TestFixture import TestFixture -from optimism import TensorMath - def make_disp_grad_from_strain(strain): return linalg.expm(strain) - np.identity(3) @@ -59,33 +57,6 @@ def test_regression_nonzero_point(self): 0.433881277313, 0.079038974065, 1.055381729505]) self.assertArrayNear(state, stateGold, 12) - def test_log_strain_updates(self): - key = jax.random.PRNGKey(1) - dispGrad = jax.random.uniform(key, (3, 3)) - initialState = self.compute_initial_state() - dt = 1.0 - - props = np.array([ - self.props['equilibrium bulk modulus'], - self.props['equilibrium shear modulus'], - self.props['non equilibrium shear modulus'], - self.props['relaxation time'] - ]) - - Ee_trial = HyperVisco._compute_elastic_logarithmic_strain(dispGrad, initialState) - state_inc = HyperVisco._compute_state_increment(Ee_trial, dt, props) - - Ee_new_way = Ee_trial - state_inc # doesn't work - - # old implementation - F = dispGrad + np.identity(3) - Fv_old = initialState.reshape((3, 3)) - Fv_new = linalg.expm(state_inc)@Fv_old - Fe = F @ np.linalg.inv(Fv_new) - Ee_old_way = TensorMath.mtk_log_sqrt(Fe.T @ Fe) - - self.assertArrayNear(Ee_new_way.ravel(), Ee_old_way.ravel(), 12) - if __name__ == '__main__': unittest.main() From b1c8519eddf088cd39346f2969336ee831e1c600 Mon Sep 17 00:00:00 2001 From: ralberd Date: Wed, 29 Nov 2023 17:38:08 -0700 Subject: [PATCH 9/9] fixing test harness name --- optimism/material/test/test_HyperVisco.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/optimism/material/test/test_HyperVisco.py b/optimism/material/test/test_HyperVisco.py index f3d23646..8f800a3e 100644 --- a/optimism/material/test/test_HyperVisco.py +++ b/optimism/material/test/test_HyperVisco.py @@ -11,7 +11,7 @@ def make_disp_grad_from_strain(strain): return linalg.expm(strain) - np.identity(3) -class GradOfPlasticityModelFixture(TestFixture): +class HyperViscoModelFixture(TestFixture): def setUp(self): G_eq = 0.855 # MPa