Merge branch 'main' into ralberd/add_snap_example

Author: ralberd

optimism/material/HyperViscoelastic.py

@@ -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,
@@ -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]
@@ -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)

test/material/test_HyperVisco.py

@@ -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,
@@ -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))
@@ -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')
@@ -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')