Merge pull request #77 from sandialabs/ralberd/material_qoi

Adding integrated material qois (dissipation for HyperVisco)

examples/viscoelastic_shear_test/plot_results.py

@@ -0,0 +1,49 @@
+from matplotlib import pyplot as plt
+import numpy as np
+
+def plot_total_work(fileName, linestyle='-'):
+    energies = np.load(fileName)
+    totalWork = energies['externalWork']
+    time = energies['time']
+    plt.plot(time, totalWork, linestyle=linestyle)
+
+def plot_algorithmic_potential(fileName, linestyle='-'):
+    energies = np.load(fileName)
+    strainEnergy = energies['algorithmicPotential']
+    time = energies['time']
+    plt.plot(time, strainEnergy, linestyle=linestyle)
+
+def plot_dissipated_energy(fileName, linestyle='-'):
+    energies = np.load(fileName)
+    dissipatedEnergy = energies['dissipation']
+    time = energies['time']
+    plt.plot(time, dissipatedEnergy, linestyle=linestyle)
+
+def plot_force_disp(fileName, linestyle='-'):
+    histories = np.load(fileName)
+    forces = histories['forces']
+    disps = histories['disps']
+    plt.plot(disps, forces, linestyle=linestyle)
+
+# plot energy histories
+energyPrefix = "energy_histories"
+forceDispPrefix = "force_disp_histories"
+fileType = ".npz"
+
+energyFile = energyPrefix + fileType
+plot_total_work(energyFile, '-')
+plot_algorithmic_potential(energyFile, ':')
+plot_dissipated_energy(energyFile, '--')
+
+plt.xlabel('Time')
+plt.ylabel('Energy (mJ)')
+plt.legend(["External Work", "Algorithmic Potential Energy", "Dissipated Energy"], loc=0, frameon=True)
+plt.savefig('energy_histories.png')
+
+plt.figure()
+forceDispFile = forceDispPrefix + fileType
+plot_force_disp(forceDispFile, '-')
+
+plt.xlabel('Displacement')
+plt.ylabel('Force')
+plt.savefig('force_disp.png')

examples/viscoelastic_shear_test/shear_test.py

@@ -0,0 +1,161 @@
+import jax.numpy as np
+from jax import jit, grad
+
+from optimism import EquationSolver as EqSolver
+from optimism import FunctionSpace
+from optimism import QuadratureRule
+from optimism.test.MeshFixture import MeshFixture
+from optimism import Mechanics
+from optimism import Mesh
+from optimism import Objective
+from optimism.material import HyperViscoelastic
+
+# This example recreates the shear test considered in
+#   Reese and Govindjee (1998). A theory of finite viscoelasticity and numerical aspects. 
+#       https://doi.org/10.1016/S0020-7683(97)00217-5
+# and in variational minimization form in 
+#   Fancello, Ponthot and Stainier (2006). A variational formulation of constitutive models and updates in non-linear finite viscoelasticity.
+#        https://doi.org/10.1002/nme.1525
+
+class ShearTest(MeshFixture):
+
+    def setUp(self):
+        dummyDispGrad = np.eye(2)
+        self.mesh = self.create_mesh_and_disp(3, 3, [0.,1.], [0.,1.],
+                                                lambda x: dummyDispGrad.dot(x))[0]
+        self.mesh = Mesh.create_higher_order_mesh_from_simplex_mesh(self.mesh, order=2, createNodeSetsFromSideSets=True)
+
+        self.quadRule = QuadratureRule.create_quadrature_rule_on_triangle(degree=2)
+        self.fs = FunctionSpace.construct_function_space(self.mesh, self.quadRule)
+
+        G_eq = 30.25 
+        K_eq = 1000*G_eq 
+        G_neq_1 = 77.77 
+        tau_1   = 17.5 
+        props = {
+            'equilibrium bulk modulus'     : K_eq,
+            'equilibrium shear modulus'    : G_eq,
+            'non equilibrium shear modulus': G_neq_1,
+            'relaxation time'              : tau_1,
+        } 
+        self.mat = HyperViscoelastic.create_material_model_functions(props)
+
+        self.EBCs = [FunctionSpace.EssentialBC(nodeSet='top', component=1),
+                FunctionSpace.EssentialBC(nodeSet='top', component=0),
+                FunctionSpace.EssentialBC(nodeSet='bottom', component=0),
+                FunctionSpace.EssentialBC(nodeSet='bottom', component=1)]
+
+        self.freq = 0.3 / 2.0 / np.pi 
+        self.cycles = 1
+        self.steps_per_cycle = 64
+        self.steps = self.cycles*self.steps_per_cycle
+        totalTime = self.cycles / self.freq
+        self.dt = totalTime / self.steps
+        self.maxDisp = 1.0
+
+    def run(self):
+        mechFuncs = Mechanics.create_mechanics_functions(self.fs,
+                                                         "plane strain",
+                                                         self.mat)
+        dofManager = FunctionSpace.DofManager(self.fs, 2, self.EBCs)
+
+        def create_field(Uu, disp):
+            def get_ubcs(disp):
+                V = np.zeros(self.mesh.coords.shape)
+                index = (self.mesh.nodeSets['top'], 0)
+                V = V.at[index].set(disp)
+                return dofManager.get_bc_values(V)
+
+            return dofManager.create_field(Uu, get_ubcs(disp))
+
+        def energy_function_all_dofs(U, p):
+            internalVariables = p.state_data
+            return mechFuncs.compute_strain_energy(U, internalVariables, self.dt)
+
+        def compute_energy(Uu, p):
+            U = create_field(Uu, p.bc_data)
+            return energy_function_all_dofs(U, p)
+
+        nodal_forces = jit(grad(energy_function_all_dofs, argnums=0))
+        integrate_dissipation = jit(mechFuncs.integrated_material_qoi)
+        integrate_free_energy = jit(mechFuncs.compute_strain_energy)
+
+        def write_output(Uu, p, step):
+            from optimism import VTKWriter
+            U = create_field(Uu, p.bc_data)
+            plotName = 'mechanics-'+str(step).zfill(3)
+            writer = VTKWriter.VTKWriter(self.mesh, baseFileName=plotName)
+
+            writer.add_nodal_field(name='displ', nodalData=U, fieldType=VTKWriter.VTKFieldType.VECTORS)
+
+            energyDensities, stresses = mechFuncs.\
+                compute_output_energy_densities_and_stresses(U, p.state_data, self.dt)
+            cellEnergyDensities = FunctionSpace.project_quadrature_field_to_element_field(self.fs, energyDensities)
+            cellStresses = FunctionSpace.project_quadrature_field_to_element_field(self.fs, stresses)
+            dissipationDensities = mechFuncs.compute_output_material_qoi(U, p.state_data, self.dt)
+            cellDissipationDensities = FunctionSpace.project_quadrature_field_to_element_field(self.fs, dissipationDensities)
+            writer.add_cell_field(name='strain_energy_density',
+                                  cellData=cellEnergyDensities,
+                                  fieldType=VTKWriter.VTKFieldType.SCALARS)
+            writer.add_cell_field(name='piola_stress',
+                                  cellData=cellStresses,
+                                  fieldType=VTKWriter.VTKFieldType.TENSORS)
+            writer.add_cell_field(name='dissipation_density',
+                                  cellData=cellDissipationDensities,
+                                  fieldType=VTKWriter.VTKFieldType.SCALARS)
+            writer.write()
+
+        Uu = dofManager.get_unknown_values(np.zeros(self.mesh.coords.shape))
+        ivs = mechFuncs.compute_initial_state()
+        p = Objective.Params(bc_data=0., state_data=ivs)
+        U = create_field(Uu, p.bc_data)
+        self.objective = Objective.Objective(compute_energy, Uu, p)
+
+        index = (self.mesh.nodeSets['top'], 0)
+
+        time = 0.0
+        times = []
+        externalWorkStore = []
+        dissipationStore = []
+        incrementalPotentialStore = []
+        forceHistory = []
+        dispHistory = []
+        for step in range(1, self.steps+1):
+            force_prev = np.array(nodal_forces(U, p).at[index].get())
+            applied_disp_prev = U.at[index].get()
+
+            disp = self.maxDisp * np.sin(2.0 * np.pi * self.freq * time)
+
+            p = Objective.param_index_update(p, 0, disp)
+            Uu, solverSuccess = EqSolver.nonlinear_equation_solve(self.objective, Uu, p, EqSolver.get_settings(), useWarmStart=False)
+            U = create_field(Uu, p.bc_data)
+            ivs = mechFuncs.compute_updated_internal_variables(U, p.state_data, self.dt)
+            p = Objective.param_index_update(p, 1, ivs)
+
+            write_output(Uu, p, step)
+
+            force = np.array(nodal_forces(U, p).at[index].get())
+            applied_disp = U.at[index].get()
+            externalWorkStore.append( 0.5*np.tensordot((force + force_prev),(applied_disp - applied_disp_prev), axes=1) )
+            incrementalPotentialStore.append(integrate_free_energy(U, ivs, self.dt))
+
+            forceHistory.append( np.sum(force) )
+            dispHistory.append(disp)
+
+            dissipationStore.append( integrate_dissipation(U, ivs, self.dt) )
+
+            times.append(time)
+            time += self.dt
+
+        # storing for plots
+        with open("energy_histories.npz",'wb') as f:
+            np.savez(f, externalWork=np.cumsum(np.array(externalWorkStore)), dissipation=np.cumsum(np.array(dissipationStore)), algorithmicPotential=np.array(incrementalPotentialStore), time=np.array(times))
+
+        with open("force_disp_histories.npz",'wb') as f:
+            np.savez(f, forces=np.array(forceHistory), disps=np.array(dispHistory))
+
+
+
+app = ShearTest()
+app.setUp()
+app.run()

optimism/Mechanics.py

@@ -13,7 +13,9 @@
                                  'compute_updated_internal_variables',
                                  'compute_element_stiffnesses',
                                  'compute_output_energy_densities_and_stresses',
-                                 'compute_initial_state'])
+                                 'compute_initial_state',
+                                 'integrated_material_qoi',
+                                 'compute_output_material_qoi'])
 
 
 DynamicsFunctions = namedtuple('DynamicsFunctions',
@@ -240,7 +242,7 @@ def compute_initial_state():
         return _compute_initial_state_multi_block(fs, materialModels)
 
 
-    return MechanicsFunctions(compute_strain_energy, jit(compute_updated_internal_variables), jit(compute_element_stiffnesses), jit(compute_output_energy_densities_and_stresses), compute_initial_state)
+    return MechanicsFunctions(compute_strain_energy, jit(compute_updated_internal_variables), jit(compute_element_stiffnesses), jit(compute_output_energy_densities_and_stresses), compute_initial_state, None, None)
 
 
 ######
@@ -292,7 +294,20 @@ def compute_initial_state():
         shape = Mesh.num_elements(fs.mesh), QuadratureRule.len(fs.quadratureRule), 1
         return np.tile(materialModel.compute_initial_state(), shape)
 
-    return MechanicsFunctions(compute_strain_energy, jit(compute_updated_internal_variables), jit(compute_element_stiffnesses), jit(compute_output_energy_densities_and_stresses), compute_initial_state)
+    def lagrangian_qoi(U, gradU, Q, X, dt):
+        return materialModel.compute_material_qoi(gradU, Q, dt)
+
+    def integrated_material_qoi(U, stateVariables, dt=0.0):
+
+        return FunctionSpace.integrate_over_block(fs, U, stateVariables, dt, 
+                                                  lagrangian_qoi,
+                                                  slice(None),
+                                                  modify_element_gradient=modify_element_gradient)
+
+    def compute_output_material_qoi(U, stateVariables, dt=0.0):
+        return FunctionSpace.evaluate_on_block(fs, U, stateVariables, dt, lagrangian_qoi, slice(None), modify_element_gradient=modify_element_gradient)
+
+    return MechanicsFunctions(compute_strain_energy, jit(compute_updated_internal_variables), jit(compute_element_stiffnesses), jit(compute_output_energy_densities_and_stresses), compute_initial_state, integrated_material_qoi, jit(compute_output_material_qoi))
 
 
 def _compute_kinetic_energy(functionSpace, V, internals, density):

optimism/inverse/AdjointFunctionSpace.py

@@ -6,9 +6,7 @@
 from optimism.FunctionSpace import map_element_shape_grads
 from jax import vmap
 
-def construct_function_space_for_adjoint(coords, mesh, quadratureRule, mode2D='cartesian'):
-
-    shapeOnRef = Interpolants.compute_shapes(mesh.parentElement, quadratureRule.xigauss)
+def construct_function_space_for_adjoint(coords, shapeOnRef, mesh, quadratureRule, mode2D='cartesian'):
 
     shapes = vmap(lambda elConns, elShape: elShape, (0, None))(mesh.conns, shapeOnRef.values)
 

optimism/inverse/test/test_Hyperelastic_gradient_checks.py

@@ -8,6 +8,7 @@
 from optimism import EquationSolver as EqSolver
 from optimism import QuadratureRule
 from optimism import FunctionSpace
+from optimism import Interpolants
 from optimism import Mechanics
 from optimism import Objective
 from optimism import Mesh
@@ -17,6 +18,11 @@
 
 from optimism.inverse import MechanicsInverse
 from optimism.inverse import AdjointFunctionSpace
+from collections import namedtuple
+
+EnergyFunctions = namedtuple('EnergyFunctions',
+                            ['energy_function_coords',
+                             'nodal_forces'])
 
 class NeoHookeanGlobalMeshAdjointSolveFixture(FiniteDifferenceFixture):
     def setUp(self):
@@ -78,17 +84,28 @@ def compute_energy(Uu, p):
 
         return storedState
 
-    def strain_energy_objective(self, storedState, parameters):
-        coords = parameters.reshape(self.mesh.coords.shape)
-        adjoint_func_space = AdjointFunctionSpace.construct_function_space_for_adjoint(coords, self.mesh, self.quadRule)
-        mech_funcs = Mechanics.create_mechanics_functions(adjoint_func_space, mode2D='plane strain', materialModel=self.materialModel)
+    def setup_energy_functions(self):
+        shapeOnRef = Interpolants.compute_shapes(self.mesh.parentElement, self.quadRule.xigauss)
+
+        def energy_function_all_dofs(U, p, coords):
+            adjoint_func_space = AdjointFunctionSpace.construct_function_space_for_adjoint(coords, shapeOnRef, self.mesh, self.quadRule)
+            mech_funcs = Mechanics.create_mechanics_functions(adjoint_func_space, mode2D='plane strain', materialModel=self.materialModel)
+            ivs = p.state_data
+            return mech_funcs.compute_strain_energy(U, ivs)
 
         def energy_function_coords(Uu, p, coords):
             U = self.dofManager.create_field(Uu, p.bc_data)
-            internal_vars = p[1]
-            return mech_funcs.compute_strain_energy(U, internal_vars)
+            return energy_function_all_dofs(U, p, coords)
+
+        nodal_forces = jax.grad(energy_function_all_dofs, argnums=0)
+
+        return EnergyFunctions(energy_function_coords, jax.jit(nodal_forces))
+
+    def strain_energy_objective(self, storedState, parameters):
+        parameters = parameters.reshape(self.mesh.coords.shape)
+        energyFuncs = self.setup_energy_functions()
 
-        return energy_function_coords(storedState[-1][0], storedState[-1][1], parameters) 
+        return energyFuncs.energy_function_coords(storedState[-1][0], storedState[-1][1], parameters) 
 
     def strain_energy_gradient(self, storedState, parameters):
         return jax.grad(self.strain_energy_objective, argnums=1)(storedState, parameters)
@@ -108,45 +125,24 @@ def total_work_increment(self, Uu, Uu_prev, p, p_prev, coords, nodal_forces):
 
     def total_work_objective(self, storedState, parameters):
         parameters = parameters.reshape(self.mesh.coords.shape)
-
-        def energy_function_all_dofs(U, p, coords):
-            adjoint_func_space = AdjointFunctionSpace.construct_function_space_for_adjoint(coords, self.mesh, self.quadRule)
-            mech_funcs = Mechanics.create_mechanics_functions(adjoint_func_space, mode2D='plane strain',materialModel=self.materialModel)
-            internal_vars = p[1]
-            return mech_funcs.compute_strain_energy(U, internal_vars)
-
-        nodal_forces = jax.jit(jax.grad(energy_function_all_dofs, argnums=0))
+        energyFuncs = self.setup_energy_functions()
 
         val = 0.0
         for step in range(1, self.steps+1):
             Uu = storedState[step][0]
             Uu_prev = storedState[step-1][0]
             p = storedState[step][1]
             p_prev = storedState[step-1][1]
-            val += self.total_work_increment(Uu, Uu_prev, p, p_prev, parameters, nodal_forces)
+            val += self.total_work_increment(Uu, Uu_prev, p, p_prev, parameters, energyFuncs.nodal_forces)
 
         return val
 
     def total_work_gradient_just_jax(self, storedState, parameters):
         return jax.grad(self.total_work_objective, argnums=1)(storedState, parameters)
 
     def total_work_gradient_with_adjoint(self, storedState, parameters):
-        def energy_function_coords(Uu, p, coords):
-            adjoint_func_space = AdjointFunctionSpace.construct_function_space_for_adjoint(coords, self.mesh, self.quadRule)
-            mech_funcs = Mechanics.create_mechanics_functions(adjoint_func_space, mode2D='plane strain', materialModel=self.materialModel)
-            U = self.dofManager.create_field(Uu, p.bc_data)
-            internal_vars = p[1]
-            return mech_funcs.compute_strain_energy(U, internal_vars)
-
-        def energy_function_all_dofs(U, p, coords):
-            adjoint_func_space = AdjointFunctionSpace.construct_function_space_for_adjoint(coords, self.mesh, self.quadRule)
-            mech_funcs = Mechanics.create_mechanics_functions(adjoint_func_space, mode2D='plane strain',materialModel=self.materialModel)
-            internal_vars = p[1]
-            return mech_funcs.compute_strain_energy(U, internal_vars)
-
-        nodal_forces = jax.jit(jax.grad(energy_function_all_dofs, argnums=0))
-
-        residualInverseFuncs = MechanicsInverse.create_residual_inverse_functions(energy_function_coords)
+        energyFuncs = self.setup_energy_functions()
+        residualInverseFuncs = MechanicsInverse.create_residual_inverse_functions(energyFuncs.energy_function_coords)
 
         compute_df = jax.grad(self.total_work_increment, (0, 1, 4))
 
@@ -161,7 +157,7 @@ def energy_function_all_dofs(U, p, coords):
             Uu_prev = storedState[step-1][0]
             p_prev = storedState[step-1][1]
 
-            df_du, df_dun, df_dx = compute_df(Uu, Uu_prev, p, p_prev, parameters, nodal_forces)
+            df_du, df_dun, df_dx = compute_df(Uu, Uu_prev, p, p_prev, parameters, energyFuncs.nodal_forces)
 
             adjointLoad -= df_du