adding test comparing dissipated energy to difference between interna…

…l and free energies at a material point

optimism/material/HyperViscoelastic.py

@@ -99,9 +99,7 @@ def _compute_dissipation(dispGrad, stateOld, dt, props):
     Ee = Ee_trial - delta_Ev 
 
     Me = 2. * G_neq * Ee
-    Me_dev = TensorMath.compute_deviatoric_tensor(Me)
-
-    return 0.5 * np.tensordot(Me, Me_dev) / eta
+    return 0.5 * np.tensordot(Me, Me) / eta
 
 def _compute_state_new(dispGrad, stateOld, dt, props):
     Ee_trial = _compute_elastic_logarithmic_strain(dispGrad, stateOld)

optimism/material/MaterialUniaxialSimulator.py

@@ -29,7 +29,7 @@
 """
 
 
-def run(materialModel, strain_history, maxTime, steps=10):
+def run(materialModel, strain_history, maxTime, steps=10, tol=1e-3):
     """ Generates the uniaxial response of a given material
     Args
     ----
@@ -57,7 +57,7 @@ def obj_func(freeStrains, p):
         strain = makeStrainTensor_(freeStrains, p)
         return energy_density(strain, p[1], dt)
 
-    uniaxialTolerance=1e-3
+    uniaxialTolerance=tol
     solverSettings = EqSolver.get_settings(tol=uniaxialTolerance)
     internalVariables =  materialModel.compute_initial_state()
     freeStrains = np.zeros(2)

optimism/material/test/test_HyperVisco.py

@@ -3,13 +3,13 @@
 import jax
 import jax.numpy as np
 from jax.scipy import linalg
+from matplotlib import pyplot as plt
 
 from optimism.material import HyperViscoelastic as HyperVisco
 from optimism.test.TestFixture import TestFixture
+from optimism.material import MaterialUniaxialSimulator
 
-def make_disp_grad_from_strain(strain):
-    return linalg.expm(strain) - np.identity(3)
-
+plotting=False
 
 class HyperViscoModelFixture(TestFixture):
     def setUp(self):
@@ -62,7 +62,80 @@ def test_regression_nonzero_point(self):
         self.assertArrayNear(state, stateGold, 12)
 
         dissipatedEnergy = self.compute_material_qoi(dispGrad, initialState, dt)
-        self.assertNear(dissipatedEnergy, 0.0, 12)
+        self.assertNear(dissipatedEnergy, 8.186677722954062, 12)
+
+class HyperViscoUniaxial(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,
+        } 
+        self.mat = HyperVisco.create_material_model_functions(self.props)
+        self.totalTime = 1.0e1*tau_1
+
+    def test_dissipated_energy_cycle(self):
+        # cyclic loading
+        stages = 2
+        maxStrain = 2.0e-2
+        maxTime = self.totalTime
+        steps = 20
+        strainInc = maxStrain/(steps-1)
+        dt=maxTime/(steps-1)
+
+        maxTime *= stages
+        steps *= stages
+
+        t0 = self.totalTime
+        def cyclic_constant_true_strain_rate(t):
+            outval = []
+            for tval in t:
+                if tval <= t0:
+                    outval.append( (tval / dt) * strainInc )
+                else:
+                    outval.append( ( (maxTime - tval) / dt) * strainInc )
+            return np.array(outval)
+
+        uniaxial = MaterialUniaxialSimulator.run(self.mat, cyclic_constant_true_strain_rate, maxTime, steps=steps, tol=1e-8)
+
+        # compute internal energy 
+        intEnergies = []
+        intEnergies.append(0.0)
+        for step in range(1,steps): 
+            F_diff = uniaxial.strainHistory[step] - uniaxial.strainHistory[step-1]
+            P_sum = uniaxial.stressHistory[step] + uniaxial.stressHistory[step-1]
+            intEnergies.append(0.5 * np.tensordot(P_sum,F_diff))
+        internalEnergyHistory = np.cumsum(np.array(intEnergies))
+
+        # compute dissipated energy
+        timePoints = np.linspace(0.0, maxTime, num=steps)
+        dt = timePoints[1] - timePoints[0]
+        dissipatedEnergies = []
+        dissipatedEnergies.append(0.0)
+        compute_dissipation = jax.jit(self.mat.compute_material_qoi)
+        for step in range(1,steps): 
+            dissipatedEnergies.append( dt * compute_dissipation(uniaxial.strainHistory[step], uniaxial.internalVariableHistory[step], dt) )
+        dissipatedEnergyHistory = np.cumsum(np.array(dissipatedEnergies))
+
+        # plot energies
+        if plotting:
+            plt.figure()
+            plt.plot(uniaxial.time, internalEnergyHistory, marker='o', fillstyle='none')
+            plt.plot(uniaxial.time, uniaxial.energyHistory, marker='x', fillstyle='none')
+            plt.plot(uniaxial.time, dissipatedEnergyHistory, marker='v', fillstyle='none')
+            plt.plot(uniaxial.time, uniaxial.energyHistory + np.cumsum(np.array(dissipatedEnergies)), marker='s', fillstyle='none')
+            plt.xlabel('Time')
+            plt.ylabel('Energy')
+            plt.legend(["Internal", "Free", "Dissipated", "Free + Dissipated"], loc=0, frameon=True)
+            plt.savefig('energy_histories.png')
+
+        self.assertArrayNear(dissipatedEnergyHistory, internalEnergyHistory - uniaxial.energyHistory, 5)
 
 
 if __name__ == '__main__':