Skip to content

Commit

Permalink
fixes bug in the dissipation potential and updates test to better check
Browse files Browse the repository at this point in the history
  • Loading branch information
@ralberd
ralberd committed Sep 11, 2024
1 parent 78dee3a commit fb0b3e8
Show file tree
Hide file tree
Showing 2 changed files with 25 additions and 37 deletions.
21 changes: 10 additions & 11 deletions optimism/material/HyperViscoelastic.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,7 @@ def compute_state_new(dispGrad, state, dt):
return state

def compute_material_qoi(dispGrad, state, dt):
return _compute_dissipation(dispGrad, state, dt, props)
return _compute_dissipated_energy(dispGrad, state, dt, props)

return MaterialModel(compute_energy_density = energy_density,
compute_initial_state = compute_initial_state,
Expand Down Expand Up @@ -63,9 +63,11 @@ def _energy_density(dispGrad, state, dt, props):
Ee = Ee_trial - delta_Ev

W_neq = _neq_strain_energy(Ee, props)
Psi = _incremental_dissipated_energy(Ee, dt, props)

Dv = delta_Ev / dt
Psi = _dissipation_potential(Dv, props)

return W_eq + W_neq + Psi
return W_eq + W_neq + dt * Psi

def _eq_strain_energy(dispGrad, props):
K, G = props[PROPS_K_eq], props[PROPS_G_eq]
Expand All @@ -81,22 +83,19 @@ def _neq_strain_energy(elasticStrain, props):
G_neq = props[PROPS_G_neq]
return G_neq * TensorMath.norm_of_deviator_squared(elasticStrain)

def _incremental_dissipated_energy(elasticStrain, dt, props):
def _dissipation_potential(Dv, props):
G_neq = props[PROPS_G_neq]
tau = props[PROPS_TAU]
eta = G_neq * tau

Me = 2. * G_neq * elasticStrain
M_bar = TensorMath.norm_of_deviator_squared(Me)

return 0.5 * dt * M_bar / eta
return eta * TensorMath.norm_of_deviator_squared(Dv)

def _compute_dissipation(dispGrad, state, dt, props):
def _compute_dissipated_energy(dispGrad, state, dt, props):
Ee_trial = _compute_elastic_logarithmic_strain(dispGrad, state)
delta_Ev = _compute_state_increment(Ee_trial, dt, props)
Ee = Ee_trial - delta_Ev

return _incremental_dissipated_energy(Ee, dt, props)
Dv = delta_Ev / dt
return dt * _dissipation_potential(Dv, props)

def _compute_state_new(dispGrad, stateOld, dt, props):
Ee_trial = _compute_elastic_logarithmic_strain(dispGrad, stateOld)
Expand Down
41 changes: 15 additions & 26 deletions test/material/test_HyperVisco.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,8 @@ def setUp(self):
G_neq_1 = 5.0
tau_1 = 25.0
self.G_neq_1 = G_neq_1
self.tau_1 = tau_1 # needed for analytic calculation below
self.tau_1 = tau_1
self.eta = G_neq_1 * tau_1
self.props = {
'equilibrium bulk modulus' : K_eq,
'equilibrium shear modulus' : G_eq,
Expand Down Expand Up @@ -52,26 +53,23 @@ def test_loading_only(self):
state_old = self.compute_initial_state()
energies = np.zeros(n_steps)
states = np.zeros((n_steps, state_old.shape[0]))
dissipated_energies = np.zeros(n_steps)

# numerical solution
for n, F in enumerate(Fs):
dispGrad = F - np.eye(3)
energies = energies.at[n].set(self.energy_density(dispGrad, state_old, dt))
state_new = self.compute_state_new(dispGrad, state_old, dt)
states = states.at[n, :].set(state_new)
dissipated_energies = dissipated_energies.at[n].set(self.compute_material_qoi(dispGrad, state_old, dt))
state_old = state_new


Fvs = jax.vmap(lambda Fv: Fv.reshape((3, 3)))(states)
Fes = jax.vmap(lambda F, Fv: F @ np.linalg.inv(Fv), in_axes=(0, 0))(Fs, Fvs)

Evs = jax.vmap(lambda Fv: log_symm(Fv))(Fvs)
Ees = jax.vmap(lambda Fe: log_symm(Fe))(Fes)

Dvs = jax.vmap(lambda Ee: (1. / self.tau_1) * deviator(Ee))(Ees)
Mes = jax.vmap(lambda Ee: 2. * self.G_neq_1 * deviator(Ee))(Ees)
dissipations = jax.vmap(lambda D, M: np.tensordot(D, M), in_axes=(0, 0))(Dvs, Mes)

# analytic solution
e_v_11 = (2. / 3.) * strain_rate * times - \
(2. / 3.) * strain_rate * self.tau_1 * (1. - np.exp(-times / self.tau_1))
Expand All @@ -86,38 +84,27 @@ def test_loading_only(self):
[0., 0., e_22]]
), in_axes=(0, 0)
)(e_e_11, e_e_22)

Me_analytic = jax.vmap(lambda Ee: 2. * self.G_neq_1 * deviator(Ee))(Ee_analytic)
Dv_analytic = jax.vmap(lambda Me: 1. / (2. * self.G_neq_1 * self.tau_1) * deviator(Me))(Me_analytic)
dissipations_analytic = jax.vmap(lambda Me, Dv: np.tensordot(Me, Dv), in_axes=(0, 0))(Me_analytic, Dv_analytic)
Dv_analytic = jax.vmap(lambda Me: 1. / (2. * self.eta) * deviator(Me))(Me_analytic)
dissipated_energies_analytic = jax.vmap(lambda Dv: dt * self.eta * np.tensordot(deviator(Dv), deviator(Dv)) )(Dv_analytic)

# test
self.assertArrayNear(Evs[:, 0, 0], e_v_11, 3)
self.assertArrayNear(Ees[:, 0, 0], e_e_11, 3)
self.assertArrayNear(Ees[:, 1, 1], e_e_22, 3)
self.assertArrayNear(dissipations, dissipations_analytic, 3)
self.assertArrayNear(dissipated_energies, dissipated_energies_analytic, 3)

if plotting:
plt.figure(1)
plt.plot(times, np.log(Fs[:, 0, 0]))
plt.xlabel('Time (s)')
plt.ylabel('F 11')
plt.savefig('times_Fs.png')

plt.figure(2)
plt.plot(times, energies)
plt.xlabel('Time (s)')
plt.ylabel('Algorithmic energy')
plt.savefig('times_energies.png')

plt.figure(3)
plt.plot(times, Evs[:, 0, 0], marker='o', linestyle='None', markevery=10)
plt.plot(times, e_v_11, linestyle='--')
plt.xlabel('Time (s)')
plt.ylabel('Viscous Logarithmic Strain 11 component')
plt.legend(["Numerical", "Analytic"], loc=0, frameon=True)
plt.savefig('times_viscous_stretch.png')

plt.figure(4)
plt.figure(2)
plt.plot(times, Ees[:, 0, 0], marker='o', linestyle='None', markevery=10, label='11', color='blue')
plt.plot(times, Ees[:, 1, 1], marker='o', linestyle='None', markevery=10, label='22', color='red')
plt.plot(times, e_e_11, linestyle='--', color='blue')
Expand All @@ -127,11 +114,13 @@ def test_loading_only(self):
plt.legend(["11 Numerical", "22 Numerical", "11 Analytic", "22 Analytic"], loc=0, frameon=True)
plt.savefig('times_elastic_stretch.png')

plt.figure(5)
plt.plot(times, dissipations, marker='o', linestyle='None', markevery=10)
plt.plot(times, dissipations_analytic, linestyle='--')
plt.figure(3)
plt.plot(times, np.cumsum(dissipated_energies), marker='o', linestyle='None', markevery=10)
plt.plot(times, np.cumsum(dissipated_energies_analytic), linestyle='--')
plt.xlabel('Time (s)')
plt.ylabel('Dissipated Energy Density (MPa)')
plt.legend(["Numerical", "Analytic"], loc=0, frameon=True)
plt.savefig('times_dissipation.png')
plt.savefig('times_dissipated_energy.png')



Expand Down

0 comments on commit fb0b3e8

Please sign in to comment.