diff --git a/examples/arch_bc/ArchBc.py b/examples/arch_bc/ArchBc.py index 895a48d8..77bae136 100644 --- a/examples/arch_bc/ArchBc.py +++ b/examples/arch_bc/ArchBc.py @@ -144,7 +144,7 @@ def run(self): p = Objective.Params(disp, ivs) precondStrategy = Objective.PrecondStrategy(self.assemble_sparse) - objective = Objective.Objective(self.energy_function, Uu, p, precondStrategy) + objective = Objective.Objective(self.energy_function, Uu, p, 1.0, precondStrategy) self.write_output(Uu, p, step=0) diff --git a/examples/arch_bc/ArchTraction.py b/examples/arch_bc/ArchTraction.py index 87933396..5a0768ea 100644 --- a/examples/arch_bc/ArchTraction.py +++ b/examples/arch_bc/ArchTraction.py @@ -1,5 +1,7 @@ import jax import jax.numpy as np +import numpy as onp +import h5py as h5p from optimism import EquationSolver as EqSolver from optimism import EquationSolverSubspace as SolverSubspace @@ -15,7 +17,7 @@ from optimism.Timer import Timer from optimism import VTKWriter from optimism.test.MeshFixture import MeshFixture - +from optimism.ReadExodusMesh import read__mesh useNewton=False if useNewton: @@ -28,32 +30,32 @@ class TractionArch(MeshFixture): def setUp(self): - self.w = 0.035 + self.w = 0.3 self.archRadius = 1.5 self.ballRadius = self.archRadius/5.0 self.initialBallLoc = self.archRadius + self.w + self.ballRadius N = 5 M = 65 - mesh, _ = \ - self.create_arch_mesh_disp_and_edges(N, M, - self.w, self.archRadius, self.w) + #mesh, _ = \ + # self.create_arch_mesh_disp_and_edges(N, M, + # self.w, self.archRadius, self.w) - mesh = Mesh.create_higher_order_mesh_from_simplex_mesh(mesh, order=2, copyNodeSets=False) + mesh = read__mesh('./foreground_mesh_shallow_arch_C.exo') + #mesh = Mesh.create_higher_order_mesh_from_simplex_mesh(mesh, order=2, copyNodeSets=False) nodeSets = Mesh.create_nodesets_from_sidesets(mesh) self.mesh = Mesh.mesh_with_nodesets(mesh, nodeSets) quadRule = QuadratureRule.create_quadrature_rule_on_triangle(degree=2) self.fs = FunctionSpace.construct_function_space(self.mesh, quadRule) - ebcs = [EssentialBC(nodeSet='left', component=0), - EssentialBC(nodeSet='left', component=1), - EssentialBC(nodeSet='right', component=0), - EssentialBC(nodeSet='right', component=1)] + ebcs = [EssentialBC(nodeSet='iside_b0_0_b1_3', component=0), + EssentialBC(nodeSet='iside_b0_0_b1_3', component=1)] + self.dofManager = DofManager(self.fs, dim=self.mesh.coords.shape[1], EssentialBCs=ebcs) self.lineQuadRule = QuadratureRule.create_quadrature_rule_1D(degree=2) - self.pushArea = (np.pi*self.archRadius/M)*self.mesh.sideSets['push'].shape[0] + self.pushArea = 0.0396*self.mesh.sideSets['iside_b0_0_b1_2'].shape[0] kappa = 10.0 nu = 0.3 @@ -73,7 +75,7 @@ def compute_energy_from_bcs(Uu, Ubc, p): strainEnergy = self.bvpFuncs.compute_strain_energy(U, internalVariables) F = p[0] loadPotential = Mechanics.compute_traction_potential_energy( - self.fs, U, self.lineQuadRule, self.mesh.sideSets['push'], + self.fs, U, self.lineQuadRule, self.mesh.sideSets['iside_b0_0_b1_2'], lambda x, n: np.array([0.0, -F/self.pushArea])) return strainEnergy + loadPotential @@ -90,7 +92,7 @@ def energy_function(self, Uu, p): internalVariables = p[1] strainEnergy = self.bvpFuncs.compute_strain_energy(U, internalVariables) F = p[0] - loadPotential = Mechanics.compute_traction_potential_energy(self.fs, U, self.lineQuadRule, self.mesh.sideSets['push'], + loadPotential = Mechanics.compute_traction_potential_energy(self.fs, U, self.lineQuadRule, self.mesh.sideSets['iside_b0_0_b1_2'], lambda x, n: np.array([0.0, -F/self.pushArea])) return strainEnergy + loadPotential @@ -106,7 +108,7 @@ def assemble_sparse(self, Uu, p): def write_output(self, Uu, p, step): U = self.create_field(Uu, p) - plotName = 'arch_traction-'+str(step).zfill(3) + plotName = 'arch_traction_C-'+str(step).zfill(3) writer = VTKWriter.VTKWriter(self.mesh, baseFileName=plotName) writer.add_nodal_field(name='displacement', nodalData=U, fieldType=VTKWriter.VTKFieldType.VECTORS) @@ -136,9 +138,11 @@ def write_output(self, Uu, p, step): force = p[0] force2 = np.sum(reactions[:,1]) print("applied force, reaction", force, force2) - disp = np.max(np.abs(U[self.mesh.nodeSets['push'],1])) + disp = np.max(np.abs(U[self.mesh.nodeSets['iside_b0_0_b1_2'],1])) self.outputForce.append(float(force)) self.outputDisp.append(float(disp)) + self.dispmax = disp + self.forcev = force with open('arch_traction_Fd.npz','wb') as f: np.savez(f, force=np.array(self.outputForce), displacement=np.array(self.outputDisp)) @@ -160,12 +164,15 @@ def run(self): p = Objective.Params(force, ivs) precondStrategy = Objective.PrecondStrategy(self.assemble_sparse) - objective = Objective.Objective(self.energy_function, Uu, p, precondStrategy) + objective = Objective.Objective(self.energy_function, Uu, p, 1.0, precondStrategy) + + self.disp = np.array([]) + self.force = np.array([]) self.write_output(Uu, p, step=0) steps = 40 - maxForce = 0.01 + maxForce = 0.1 for i in range(1, steps): print('--------------------------------------') print('LOAD STEP ', i) @@ -175,8 +182,16 @@ def run(self): Uu, solverSuccess = EqSolver.nonlinear_equation_solve(objective, Uu, p, self.trSettings, solver_algorithm=solver) self.write_output(Uu, p, i) + self.disp = onp.array(np.append(self.disp,self.dispmax)) + self.force = onp.array(np.append(self.force,self.forcev)) + + dispf = onp.append(self.disp,self.force,axis=0) + + T_matrix_file_2 = h5p.File('./FDPlot_Buckling_Arch_BF_C.h5','w') + T_matrix_file_2.create_dataset('FD', data=dispf) + T_matrix_file_2.close() - unload = True + unload = False if unload: for i in range(steps, 2*steps - 1): print('--------------------------------------') diff --git a/examples/arch_bc/ArchTraction_immersed.py b/examples/arch_bc/ArchTraction_immersed.py new file mode 100644 index 00000000..d63b1f9b --- /dev/null +++ b/examples/arch_bc/ArchTraction_immersed.py @@ -0,0 +1,230 @@ +import jax +import jax.numpy as np +import h5py as h5p +import numpy as onp +from scipy import sparse +from optimism import EquationSolver_Immersed_2 as EqSolver +from optimism import EquationSolverSubspace as SolverSubspace +from optimism import FunctionSpace +from optimism.material import Neohookean as MatModel +from optimism import Mechanics +from optimism import Mesh +from optimism.FunctionSpace import EssentialBC +from optimism.FunctionSpace import DofManager +from optimism import Objective +from optimism import SparseMatrixAssembler +from optimism import QuadratureRule +from optimism.Timer import Timer +from optimism import VTKWriter +from optimism.test.MeshFixture import MeshFixture +from optimism.ReadExodusMesh import read__mesh + +useNewton=False + +if useNewton: + solver = EqSolver.newton +else: + solver = EqSolver.trust_region_minimize + + + +class TractionArch(MeshFixture): + + def setUp(self): + T_matrix_file = h5p.File('./Extraction_Buckling_Arch_A.h5','r') + T = np.transpose(np.array(T_matrix_file['T'])) + self.Tm = T + self.Tmn = onp.array(T) + self.Ttn = onp.transpose(onp.array(T)) + self.w = 0.3 + self.archRadius = 1.5 + self.ballRadius = self.archRadius/5.0 + self.initialBallLoc = self.archRadius + self.w + self.ballRadius + N = 5 + M = 65 + + #mesh, _ = \ + # self.create_arch_mesh_disp_and_edges(N, M, + # self.w, self.archRadius, self.w) + + mesh = read__mesh('./foreground_mesh_shallow_arch_A.exo') + #mesh = Mesh.create_higher_order_mesh_from_simplex_mesh(mesh, order=2, copyNodeSets=False) + nodeSets = Mesh.create_nodesets_from_sidesets(mesh) + self.mesh = Mesh.mesh_with_nodesets(mesh, nodeSets) + + quadRule = QuadratureRule.create_quadrature_rule_on_triangle(degree=2) + self.fs = FunctionSpace.construct_function_space(self.mesh, quadRule) + + ebcs = [EssentialBC(nodeSet='iside_b0_0_b1_3', component=0), + EssentialBC(nodeSet='iside_b0_0_b1_3', component=1)] + + self.dofManager = DofManager(self.fs, dim=self.mesh.coords.shape[1], EssentialBCs=[]) + + self.lineQuadRule = QuadratureRule.create_quadrature_rule_1D(degree=2) + self.pushArea = 0.0396*self.mesh.sideSets['iside_b0_0_b1_2'].shape[0] + + kappa = 10.0 + nu = 0.3 + E = 3*kappa*(1 - 2*nu) + props = {'elastic modulus': E, + 'poisson ratio': nu, + 'version': 'coupled'} + materialModel = MatModel.create_material_model_functions(props) + + self.bvpFuncs = Mechanics.create_mechanics_functions(self.fs, + mode2D="plane strain", + materialModel=materialModel) + + def compute_energy_from_bcs(Uu, Ubc, p): + U = self.dofManager.create_field(Uu, Ubc) + internalVariables = p[1] + strainEnergy = self.bvpFuncs.compute_strain_energy(U, internalVariables) + F = p[0] + loadPotential = Mechanics.compute_traction_potential_energy( + self.fs, U, self.lineQuadRule, self.mesh.sideSets['iside_b0_0_b1_2'], + lambda x, n: np.array([0.0, -F/self.pushArea])) + return strainEnergy + loadPotential + + #self.compute_bc_reactions = jax.jit(jax.grad(compute_energy_from_bcs, 1)) + + self.trSettings = EqSolver.get_settings(max_trust_iters=400, t1=0.4, t2=1.5, eta1=1e-8, eta2=0.2, eta3=0.8, over_iters=100) + + self.outputForce = [0.0] + self.outputDisp = [0.0] + + + def energy_function(self, Uu, p): + U = self.create_field(Uu, p) + internalVariables = p[1] + strainEnergy = self.bvpFuncs.compute_strain_energy(U, internalVariables) + F = p[0] + loadPotential = Mechanics.compute_traction_potential_energy(self.fs, U, self.lineQuadRule, self.mesh.sideSets['iside_b0_0_b1_2'], + lambda x, n: np.array([0.0, -F/self.pushArea])) + displacementPotential1 = Mechanics.compute_displacement_penalty(self.fs, U, self.lineQuadRule, self.mesh.sideSets['iside_b0_0_b1_3'], + lambda x, n: np.array([0.0, 0.0]),0,100) + displacementPotential2 = Mechanics.compute_displacement_penalty(self.fs, U, self.lineQuadRule, self.mesh.sideSets['iside_b0_0_b1_3'], + lambda x, n: np.array([0.0, 0.0]),1,100) + + return strainEnergy + loadPotential + displacementPotential1 + displacementPotential2 + + + def assemble_sparse(self, Uu, p): + U = self.create_field(Uu, p) + internalVariables = p[1] + elementStiffnesses = self.bvpFuncs.compute_element_stiffnesses(U, internalVariables) + d = SparseMatrixAssembler.assemble_sparse_stiffness_matrix(elementStiffnesses, + self.fs.mesh.conns, + self.dofManager).todense() # Background + + prec = sparse.csc_matrix(np.dot(np.dot(np.transpose(self.Tmn),d),self.Tmn)) + + + return prec + + + def write_output(self, Uu, p, step): + U = self.create_field(Uu, p) + + plotName = 'arch_traction_immersed_A-'+str(step).zfill(3) + writer = VTKWriter.VTKWriter(self.mesh, baseFileName=plotName) + + writer.add_nodal_field(name='displacement', nodalData=U, fieldType=VTKWriter.VTKFieldType.VECTORS) + + bcs = np.array(self.dofManager.isBc, dtype=int) + writer.add_nodal_field(name='bcs', nodalData=bcs, fieldType=VTKWriter.VTKFieldType.VECTORS, dataType=VTKWriter.VTKDataType.INT) + + Ubc = self.get_ubcs(p) + internalVariables = p[1] + #rxnBc = self.compute_bc_reactions(Uu, Ubc, p) + #reactions = np.zeros(U.shape).at[self.dofManager.isBc].set(rxnBc) + #writer.add_nodal_field(name='reactions', nodalData=reactions, fieldType=VTKWriter.VTKFieldType.VECTORS) + + energyDensities, stresses = self.bvpFuncs.\ + compute_output_energy_densities_and_stresses(U, internalVariables) + cellEnergyDensities = FunctionSpace.project_quadrature_field_to_element_field(self.fs, energyDensities) + cellStresses = FunctionSpace.project_quadrature_field_to_element_field(self.fs, stresses) + 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.write() + + force = p[0] + #force2 = np.sum(reactions[:,1]) + print("applied force, reaction", force, 0) + disp = np.max(np.abs(U[self.mesh.nodeSets['iside_b0_0_b1_2'],1])) + self.outputForce.append(float(force)) + self.outputDisp.append(float(disp)) + self.dispmax = disp + self.forcev = force + + with open('arch_traction_Fd.npz','wb') as f: + np.savez(f, force=np.array(self.outputForce), displacement=np.array(self.outputDisp)) + + + def get_ubcs(self, p): + V = np.zeros(self.mesh.coords.shape) + return self.dofManager.get_bc_values(V) + + + def create_field(self, Uu, p): + return self.dofManager.create_field(Uu, self.get_ubcs(p)) + + + def run(self): + Uu = np.zeros(np.shape(self.Tm)[1])#self.dofManager.get_unknown_values(np.zeros(self.mesh.coords.shape)) + force = 0.0 + print('DoFs') + print(np.shape(self.Tm)[0]) + ivs = self.bvpFuncs.compute_initial_state() + p = Objective.Params(force, ivs) + + precondStrategy = Objective.PrecondStrategy(self.assemble_sparse) + objective = Objective.Objective(self.energy_function, np.dot(self.Tm,Uu), p, 1.0, precondStrategy) + + self.write_output(np.dot(self.Tm,Uu), p, step=0) + + steps = 40 + maxForce = 0.1 + self.disp = np.array([]) + self.force = np.array([]) + for i in range(1, steps): + print('--------------------------------------') + print('LOAD STEP ', i) + + force += maxForce/steps + print('Force updated, proceeding to solver') + p = Objective.param_index_update(p, 0, force) + + Uu, solverSuccess, diff = EqSolver.nonlinear_equation_solve(objective, Uu, p, self.Tm, self.trSettings, solver_algorithm=solver) + + self.write_output(np.dot(self.Tm,Uu), p, i) + self.disp = onp.array(np.append(self.disp,self.dispmax)) + self.force = onp.array(np.append(self.force,self.forcev)) + + dispf = onp.append(self.disp,self.force,axis=0) + + T_matrix_file_2 = h5p.File('./FDPlot_Buckling_Arch_A.h5','w') + T_matrix_file_2.create_dataset('FD', data=dispf) + T_matrix_file_2.close() + + unload = False + if unload: + for i in range(steps, 2*steps - 1): + print('--------------------------------------') + print('LOAD STEP ', i) + + force -= maxForce/steps + p = Objective.param_index_update(p, 0, force) + Uu, solverSuccess = EqSolver.nonlinear_equation_solve(objective, Uu, p, self.trSettings, solver_algorithm=solver) + + self.write_output(Uu, p, i) + +app = TractionArch() +app.setUp() +with Timer(name="AppRun"): + app.run() + diff --git a/examples/arch_bc/Extraction_Buckling_Arch.h5 b/examples/arch_bc/Extraction_Buckling_Arch.h5 new file mode 100644 index 00000000..82620084 Binary files /dev/null and b/examples/arch_bc/Extraction_Buckling_Arch.h5 differ diff --git a/examples/arch_bc/Extraction_Buckling_Arch_A.h5 b/examples/arch_bc/Extraction_Buckling_Arch_A.h5 new file mode 100644 index 00000000..451a762b Binary files /dev/null and b/examples/arch_bc/Extraction_Buckling_Arch_A.h5 differ diff --git a/examples/arch_bc/Extraction_Buckling_Arch_B.h5 b/examples/arch_bc/Extraction_Buckling_Arch_B.h5 new file mode 100644 index 00000000..cf9abe25 Binary files /dev/null and b/examples/arch_bc/Extraction_Buckling_Arch_B.h5 differ diff --git a/examples/arch_bc/FDPlot_Buckling_Arch_A.h5 b/examples/arch_bc/FDPlot_Buckling_Arch_A.h5 new file mode 100644 index 00000000..e7bcdef2 Binary files /dev/null and b/examples/arch_bc/FDPlot_Buckling_Arch_A.h5 differ diff --git a/examples/arch_bc/FDPlot_Buckling_Arch_B.h5 b/examples/arch_bc/FDPlot_Buckling_Arch_B.h5 new file mode 100644 index 00000000..afac9cbe Binary files /dev/null and b/examples/arch_bc/FDPlot_Buckling_Arch_B.h5 differ diff --git a/examples/arch_bc/FDPlot_Buckling_Arch_BF_A.h5 b/examples/arch_bc/FDPlot_Buckling_Arch_BF_A.h5 new file mode 100644 index 00000000..2e37cd8a Binary files /dev/null and b/examples/arch_bc/FDPlot_Buckling_Arch_BF_A.h5 differ diff --git a/examples/arch_bc/FDPlot_Buckling_Arch_BF_B.h5 b/examples/arch_bc/FDPlot_Buckling_Arch_BF_B.h5 new file mode 100644 index 00000000..a569dc90 Binary files /dev/null and b/examples/arch_bc/FDPlot_Buckling_Arch_BF_B.h5 differ diff --git a/examples/arch_bc/FDPlot_Buckling_Arch_BF_C.h5 b/examples/arch_bc/FDPlot_Buckling_Arch_BF_C.h5 new file mode 100644 index 00000000..042f68e9 Binary files /dev/null and b/examples/arch_bc/FDPlot_Buckling_Arch_BF_C.h5 differ diff --git a/examples/arch_bc/FDPlot_Buckling_Arch_C.h5 b/examples/arch_bc/FDPlot_Buckling_Arch_C.h5 new file mode 100644 index 00000000..f6a5f14e Binary files /dev/null and b/examples/arch_bc/FDPlot_Buckling_Arch_C.h5 differ diff --git a/examples/arch_bc/foreground_mesh_shallow_arch.exo b/examples/arch_bc/foreground_mesh_shallow_arch.exo new file mode 100644 index 00000000..a795645a Binary files /dev/null and b/examples/arch_bc/foreground_mesh_shallow_arch.exo differ diff --git a/examples/arch_bc/foreground_mesh_shallow_arch_A.exo b/examples/arch_bc/foreground_mesh_shallow_arch_A.exo new file mode 100644 index 00000000..a795645a Binary files /dev/null and b/examples/arch_bc/foreground_mesh_shallow_arch_A.exo differ diff --git a/examples/arch_bc/foreground_mesh_shallow_arch_B.exo b/examples/arch_bc/foreground_mesh_shallow_arch_B.exo new file mode 100644 index 00000000..b09a5c6c Binary files /dev/null and b/examples/arch_bc/foreground_mesh_shallow_arch_B.exo differ diff --git a/examples/arch_bc/foreground_mesh_shallow_arch_C.exo b/examples/arch_bc/foreground_mesh_shallow_arch_C.exo new file mode 100644 index 00000000..2e67a090 Binary files /dev/null and b/examples/arch_bc/foreground_mesh_shallow_arch_C.exo differ diff --git a/examples/arch_bc/foreground_mesh_shallow_arch_D.exo b/examples/arch_bc/foreground_mesh_shallow_arch_D.exo new file mode 100644 index 00000000..ac22f6d8 Binary files /dev/null and b/examples/arch_bc/foreground_mesh_shallow_arch_D.exo differ diff --git a/examples/euler_buckle/EulerBuckleImmersed_Buckling.py b/examples/euler_buckle/EulerBuckleImmersed_Buckling.py new file mode 100644 index 00000000..1b3d8391 --- /dev/null +++ b/examples/euler_buckle/EulerBuckleImmersed_Buckling.py @@ -0,0 +1,243 @@ +import jax +import jax.numpy as np +import h5py as h5p +import numpy as onp +from scipy import sparse + +from optimism import EquationSolver_Immersed_2 as EqSolver +from optimism import FunctionSpace +from optimism.material import Neohookean as MatModel +from optimism import Mechanics +from optimism import Mesh +from optimism.FunctionSpace import EssentialBC +from optimism.FunctionSpace import DofManager +from optimism import Objective +from optimism import SparseMatrixAssembler +from optimism import QuadratureRule +from optimism.Timer import Timer +from optimism import VTKWriter +from optimism.test.MeshFixture import MeshFixture +from optimism.ReadExodusMesh import read__mesh + + +useNewton=False + +if useNewton: + solver = EqSolver.newton +else: + solver = EqSolver.trust_region_minimize + +class EulerBeam(MeshFixture): + + def setUp(self): + # No mesh generation. We are reading the mesh from MORIS, in addition to the extraction operator. + T_matrix_file = h5p.File('./Extraction_Buckling_Optimism_C.h5','r') + T = np.transpose(np.array(T_matrix_file['T'])) + print(np.shape(T)) + self.w = 0.25 + #self.Tm = np.identity(715) + self.Tm = T + self.Tmn = onp.array(self.Tm) + #h = 1.5 + #N = 5 + #M = 45 + #mesh, _ = self.create_mesh_and_disp(N, M, [-self.w/2,self.w/2], [0., h], lambda x: None) # Mesh set-up + #mesh = Mesh.create_higher_order_mesh_from_simplex_mesh(mesh, order=2, copyNodeSets=False) # hp-refinement + #nodeSets = Mesh.create_nodesets_from_sidesets(mesh) # BC Patches? + #self.mesh = Mesh.mesh_with_nodesets(mesh, nodeSets) # Integrate BC Patches ? + mesh = read__mesh('./foreground_mesh_optimism_buckle_C.exo') + nodeSets = Mesh.create_nodesets_from_sidesets(mesh) + self.mesh = Mesh.mesh_with_nodesets(mesh, nodeSets) + + + # No change here. + quadRule = QuadratureRule.create_quadrature_rule_on_triangle(degree=2) # Quadrature + self.fs = FunctionSpace.construct_function_space(self.mesh, quadRule) + + ebcs = [EssentialBC(nodeSet='iside_b0_0_b1_1', component=0), + EssentialBC(nodeSet='iside_b0_0_b1_2', component=0), + EssentialBC(nodeSet='iside_b0_0_b1_2', component=1)] + + # Replace EBCs with empty array, because we are enforcing dirichlet BCs via Penalty. + self.dofManager = DofManager(self.fs, dim=self.mesh.coords.shape[1], EssentialBCs=[]) + + # Quad + self.lineQuadRule = QuadratureRule.create_quadrature_rule_1D(degree=2) # quad + + kappa = 10.0 + nu = 0.3 + E = 3*kappa*(1 - 2*nu) + props = {'elastic modulus': E, + 'poisson ratio': nu, + 'version': 'coupled'} + materialModel = MatModel.create_material_model_functions(props) + + self.bvpFuncs = Mechanics.create_mechanics_functions(self.fs, + mode2D="plane strain", + materialModel=materialModel) # Potential stuff again? + + #self.compute_bc_reactions = jax.jit(jax.grad(self.energy_from_bcs, 1)) # ??? gradient of potential + + self.trSettings = EqSolver.get_settings(max_trust_iters=100, t1=0.4, t2=1.5, eta1=1e-5, eta2=0.2, eta3=0.8, over_iters=500) # Solver + + self.outputForce = [0.0] + self.outputDisp = [0.0] + + # Hit displacement with extraction operator in the main run function, so that the mapped displacement is only passed to all these. + + def energy_function_from_full_field(self, U, p): + # Add penalty here. Replace with correct sideset name. + + internalVariables = p[1] + strainEnergy = self.bvpFuncs.compute_strain_energy(U, internalVariables) + F = p[0] + + loadPotential = Mechanics.compute_traction_potential_energy(self.fs, U, self.lineQuadRule, self.mesh.sideSets['iside_b0_0_b1_1'], + lambda x, n: np.array([0.0, -F/self.w])) + displacementPotential1 = Mechanics.compute_displacement_penalty(self.fs, U, self.lineQuadRule, self.mesh.sideSets['iside_b0_0_b1_2'], + lambda x, n: np.array([0.0, 0.0]),0,100) + displacementPotential2 = Mechanics.compute_displacement_penalty(self.fs, U, self.lineQuadRule, self.mesh.sideSets['iside_b0_0_b1_2'], + lambda x, n: np.array([0.0, 0.0]),1,100) + displacementPotential3 = Mechanics.compute_displacement_penalty(self.fs, U, self.lineQuadRule, self.mesh.sideSets['iside_b0_0_b1_1'], + lambda x, n: np.array([0.0, 0.0]),0,100) + + + + + return strainEnergy + loadPotential + displacementPotential1 + displacementPotential2 + displacementPotential3 # potential to minimize... + + # This one, not sure, but not required in the main forward solve. + + def energy_from_bcs(self, Uu, Ubc, p): + U = self.dofManager.create_field(Uu, Ubc) + return self.energy_function_from_full_field(U, p) # adding BC contributions + + def energy_function(self, Uu, p): + U = self.create_field(Uu, p) + return self.energy_function_from_full_field(U, p) + + def assemble_sparse(self, Uu, p): + # No problem here. Just hit with extraction operator after the global stiffness matrix a.k.a. preconditioner inv computed. + U = self.create_field(Uu, p) + internalVariables = p[1] + elementStiffnesses = self.bvpFuncs.compute_element_stiffnesses(U, internalVariables) + d = SparseMatrixAssembler.assemble_sparse_stiffness_matrix(elementStiffnesses, + self.fs.mesh.conns, + self.dofManager).todense() # Background + + prec = sparse.csc_matrix(np.dot(np.dot(np.transpose(self.Tmn),d),self.Tmn)) + + + return prec + + def write_output(self, Uu, p, step): + U = self.create_field(Uu, p) + plotName = 'euler_column_immersed_2_buckle_optimism_C_clip-'+str(step).zfill(3) + writer = VTKWriter.VTKWriter(self.mesh, baseFileName=plotName) + + writer.add_nodal_field(name='displacement', nodalData=U, fieldType=VTKWriter.VTKFieldType.VECTORS) + + bcs = np.array(self.dofManager.isBc, dtype=int) + writer.add_nodal_field(name='bcs', nodalData=bcs, fieldType=VTKWriter.VTKFieldType.VECTORS, dataType=VTKWriter.VTKDataType.INT) + + Ubc = self.get_ubcs(p) + internalVariables = p[1] + #rxnBc = self.compute_bc_reactions(Uu, Ubc, p) + #reactions = np.zeros(U.shape).at[self.dofManager.isBc].set(rxnBc) + #writer.add_nodal_field(name='reactions', nodalData=reactions, fieldType=VTKWriter.VTKFieldType.VECTORS) + + energyDensities, stresses = self.bvpFuncs.compute_output_energy_densities_and_stresses(U, internalVariables) + cellEnergyDensities = FunctionSpace.project_quadrature_field_to_element_field(self.fs, energyDensities) + cellStresses = FunctionSpace.project_quadrature_field_to_element_field(self.fs, stresses) + 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.write() + + force = p[0] + #force2 = np.sum(reactions[:,1]) + print("applied force, reaction", force) + disp = np.max(np.abs(U[self.mesh.sideSets['iside_b0_0_b1_1'],1])) + print('Max Displacement') + print(disp) + self.outputForce.append(float(force)) + self.outputDisp.append(float(disp)) + + with open('column_Fd_immersed.npz','wb') as f: + np.savez(f, force=np.array(self.outputForce), displacement=np.array(self.outputDisp)) + + + def get_ubcs(self, p): + V = np.zeros(self.mesh.coords.shape) + return self.dofManager.get_bc_values(V) # No issues. + + + def create_field(self, Uu, p): + return self.dofManager.create_field(Uu, self.get_ubcs(p)) + + + + + def run(self): + #T_matrix_file = h5p.File('./Extraction.h5','r') + #T = np.transpose(np.array(T_matrix_file['T'])) + # Now Uu is defined in the background. + #Uu = self.dofManager.get_unknown_values(np.zeros(self.mesh.coords.shape)) + #print(np.shape(self.Tm)[0]) + Uu = np.zeros(np.shape(self.Tm)[1]) + print('DoFs') + print(np.shape(self.Tm)[0]) + force = 0.0 + ivs = self.bvpFuncs.compute_initial_state() + p = Objective.Params(force, ivs) + + precondStrategy = Objective.PrecondStrategy(self.assemble_sparse) + + # Project to foreground + objective = Objective.Objective(self.energy_function, np.dot(self.Tm,Uu), p, self.Tm, precondStrategy) + + #Project to foreground + + self.write_output(np.dot(self.Tm,Uu), p, step=0) + + steps = 5 + maxForce = 0.02 + for i in range(1, steps): + key = jax.random.key(0) + init_guess = Uu+jax.random.uniform(key=key,shape=np.shape(Uu)[0],minval=0,maxval=0.001) + print('--------------------------------------') + print('LOAD STEP ', i) + force += maxForce/steps + p = Objective.param_index_update(p, 0, force) + Uu, solverSuccess, diff = EqSolver.nonlinear_equation_solve(objective, init_guess, p, self.Tm, self.trSettings, solver_algorithm=solver) + # Project to foreground + self.write_output(np.dot(self.Tm,Uu), p, i) + dst = np.linalg.solve(np.dot(np.transpose(self.Tm),(self.Tm)),(np.dot(np.transpose(self.Tm),np.dot(self.Tm,Uu)))) + print("Norm Value") + Tt = np.transpose(self.Tm) + # Test whether solver is in background space or not. + print(np.linalg.norm((dst)-Uu)) + # Check the scaling + print(np.linalg.norm(diff-np.ones(np.shape(self.Tm)[1]))) + + + unload = False + if unload: + for i in range(steps, 2*steps - 1): + print('--------------------------------------') + print('LOAD STEP ', i) + force -= maxForce/steps + p = Objective.param_index_update(p, 0, force) + Uu, solverSuccess = EqSolver.nonlinear_equation_solve(objective, Uu, p, self.Tm, self.trSettings, solver_algorithm=solver) + + self.write_output(Uu, p, i) + +app = EulerBeam() +app.setUp() +with Timer(name="AppRun"): + app.run() + diff --git a/examples/euler_buckle/EulerBuckle.py b/examples/euler_buckle/EulerBuckle_BodyFitted.py similarity index 84% rename from examples/euler_buckle/EulerBuckle.py rename to examples/euler_buckle/EulerBuckle_BodyFitted.py index dc230457..bcb48786 100644 --- a/examples/euler_buckle/EulerBuckle.py +++ b/examples/euler_buckle/EulerBuckle_BodyFitted.py @@ -1,6 +1,6 @@ import jax import jax.numpy as np - +import numpy as onp from optimism import EquationSolver as EqSolver from optimism import FunctionSpace from optimism.material import Neohookean as MatModel @@ -14,6 +14,7 @@ from optimism.Timer import Timer from optimism import VTKWriter from optimism.test.MeshFixture import MeshFixture +from optimism.ReadExodusMesh import read__mesh useNewton=False @@ -25,21 +26,21 @@ class EulerBeam(MeshFixture): def setUp(self): - self.w = 0.1 - h = 1.5 - N = 5 - M = 45 - mesh, _ = self.create_mesh_and_disp(N, M, [-self.w/2,self.w/2], [0., h], lambda x: None) - mesh = Mesh.create_higher_order_mesh_from_simplex_mesh(mesh, order=2, copyNodeSets=False) + self.w = 0.25 + #h = 20 + #N = 5 + #M = 80 + #mesh, _ = self.create_mesh_and_disp(N, M, [-self.w/2,self.w/2], [0., h], lambda x: None) + mesh = read__mesh('./foreground_mesh_optimism_buckle_D.exo')#Mesh.create_higher_order_mesh_from_simplex_mesh(mesh, order=2, copyNodeSets=False) nodeSets = Mesh.create_nodesets_from_sidesets(mesh) self.mesh = Mesh.mesh_with_nodesets(mesh, nodeSets) - + quadRule = QuadratureRule.create_quadrature_rule_on_triangle(degree=2) self.fs = FunctionSpace.construct_function_space(self.mesh, quadRule) - ebcs = [EssentialBC(nodeSet='top', component=0), - EssentialBC(nodeSet='bottom', component=0), - EssentialBC(nodeSet='bottom', component=1)] + ebcs = [EssentialBC(nodeSet='iside_b0_0_b1_1', component=0), + EssentialBC(nodeSet='iside_b0_0_b1_2', component=0), + EssentialBC(nodeSet='iside_b0_0_b1_2', component=1)] self.dofManager = DofManager(self.fs, dim=self.mesh.coords.shape[1], EssentialBCs=ebcs) @@ -68,7 +69,7 @@ def energy_function_from_full_field(self, U, p): internalVariables = p[1] strainEnergy = self.bvpFuncs.compute_strain_energy(U, internalVariables) F = p[0] - loadPotential = Mechanics.compute_traction_potential_energy(self.fs, U, self.lineQuadRule, self.mesh.sideSets['top'], + loadPotential = Mechanics.compute_traction_potential_energy(self.fs, U, self.lineQuadRule, self.mesh.sideSets['iside_b0_0_b1_1'], lambda x, n: np.array([0.0, -F/self.w])) return strainEnergy + loadPotential @@ -90,7 +91,7 @@ def assemble_sparse(self, Uu, p): def write_output(self, Uu, p, step): U = self.create_field(Uu, p) - plotName = 'euler_column-'+str(step).zfill(3) + plotName = 'euler_column_D-'+str(step).zfill(3) writer = VTKWriter.VTKWriter(self.mesh, baseFileName=plotName) writer.add_nodal_field(name='displacement', nodalData=U, fieldType=VTKWriter.VTKFieldType.VECTORS) @@ -119,9 +120,11 @@ def write_output(self, Uu, p, step): force = p[0] force2 = np.sum(reactions[:,1]) print("applied force, reaction", force, force2) - disp = np.max(np.abs(U[self.mesh.nodeSets['top'],1])) + disp = np.max(np.abs(U[self.mesh.nodeSets['iside_b0_0_b1_1'],1])) self.outputForce.append(float(force)) self.outputDisp.append(float(disp)) + print('Max Displacement') + print(disp) with open('column_Fd.npz','wb') as f: np.savez(f, force=np.array(self.outputForce), displacement=np.array(self.outputDisp)) @@ -143,22 +146,24 @@ def run(self): p = Objective.Params(force, ivs) precondStrategy = Objective.PrecondStrategy(self.assemble_sparse) - objective = Objective.Objective(self.energy_function, Uu, p, precondStrategy) + objective = Objective.Objective(self.energy_function, Uu, p, 1.0, precondStrategy) self.write_output(Uu, p, step=0) - steps = 30 - maxForce = 0.05 + steps = 2 + maxForce = 0.0022 for i in range(1, steps): + key = jax.random.key(0) + init_guess = Uu+jax.random.uniform(key=key,shape=np.shape(Uu)[0],minval=0,maxval=0.001) print('--------------------------------------') print('LOAD STEP ', i) force += maxForce/steps p = Objective.param_index_update(p, 0, force) - Uu, solverSuccess = EqSolver.nonlinear_equation_solve(objective, Uu, p, self.trSettings, solver_algorithm=solver) + Uu, solverSuccess = EqSolver.nonlinear_equation_solve(objective, init_guess, p, self.trSettings, solver_algorithm=solver) self.write_output(Uu, p, i) - unload = True + unload = False if unload: for i in range(steps, 2*steps - 1): print('--------------------------------------') diff --git a/examples/euler_buckle/Extraction_Buckling_Optimism_2.h5 b/examples/euler_buckle/Extraction_Buckling_Optimism_2.h5 new file mode 100644 index 00000000..e646711f Binary files /dev/null and b/examples/euler_buckle/Extraction_Buckling_Optimism_2.h5 differ diff --git a/examples/euler_buckle/Extraction_Buckling_Optimism_A.h5 b/examples/euler_buckle/Extraction_Buckling_Optimism_A.h5 new file mode 100644 index 00000000..4a76eeb3 Binary files /dev/null and b/examples/euler_buckle/Extraction_Buckling_Optimism_A.h5 differ diff --git a/examples/euler_buckle/Extraction_Buckling_Optimism_B.h5 b/examples/euler_buckle/Extraction_Buckling_Optimism_B.h5 new file mode 100644 index 00000000..9feceb9d Binary files /dev/null and b/examples/euler_buckle/Extraction_Buckling_Optimism_B.h5 differ diff --git a/examples/euler_buckle/Extraction_Buckling_Optimism_C.h5 b/examples/euler_buckle/Extraction_Buckling_Optimism_C.h5 new file mode 100644 index 00000000..0446b917 Binary files /dev/null and b/examples/euler_buckle/Extraction_Buckling_Optimism_C.h5 differ diff --git a/examples/euler_buckle/foreground_mesh_optimism_buckle_2.exo b/examples/euler_buckle/foreground_mesh_optimism_buckle_2.exo new file mode 100644 index 00000000..9de362c0 Binary files /dev/null and b/examples/euler_buckle/foreground_mesh_optimism_buckle_2.exo differ diff --git a/examples/euler_buckle/foreground_mesh_optimism_buckle_A.exo b/examples/euler_buckle/foreground_mesh_optimism_buckle_A.exo new file mode 100644 index 00000000..d4730c86 Binary files /dev/null and b/examples/euler_buckle/foreground_mesh_optimism_buckle_A.exo differ diff --git a/examples/euler_buckle/foreground_mesh_optimism_buckle_B.exo b/examples/euler_buckle/foreground_mesh_optimism_buckle_B.exo new file mode 100644 index 00000000..1cd2ff3d Binary files /dev/null and b/examples/euler_buckle/foreground_mesh_optimism_buckle_B.exo differ diff --git a/examples/euler_buckle/foreground_mesh_optimism_buckle_C.exo b/examples/euler_buckle/foreground_mesh_optimism_buckle_C.exo new file mode 100644 index 00000000..d069acf7 Binary files /dev/null and b/examples/euler_buckle/foreground_mesh_optimism_buckle_C.exo differ diff --git a/examples/euler_buckle/foreground_mesh_optimism_buckle_D.exo b/examples/euler_buckle/foreground_mesh_optimism_buckle_D.exo new file mode 100644 index 00000000..4de3c356 Binary files /dev/null and b/examples/euler_buckle/foreground_mesh_optimism_buckle_D.exo differ diff --git a/optimism/AlSolver.py b/optimism/AlSolver.py index e573da54..d7e4db77 100644 --- a/optimism/AlSolver.py +++ b/optimism/AlSolver.py @@ -1,5 +1,5 @@ from optimism.JaxConfig import * -import optimism.EquationSolver as EqSolver +import EquationSolver_Immersed_2 as EqSolver from optimism.NewtonSolver import globalized_newton_step from optimism.NewtonSolver import newton_step from optimism import WarmStart diff --git a/optimism/BoundConstrainedSolver.py b/optimism/BoundConstrainedSolver.py index 0216c8c8..66ad02ee 100644 --- a/optimism/BoundConstrainedSolver.py +++ b/optimism/BoundConstrainedSolver.py @@ -3,7 +3,7 @@ #from optimism import ConstrainedObjective from optimism import AlSolver from optimism import WarmStart -from optimism import EquationSolver +from optimism import EquationSolver_Immersed_2 def bound_constrained_solve(boundConstrainedObjective, @@ -13,7 +13,7 @@ def bound_constrained_solve(boundConstrainedObjective, sub_problem_callback=None, useWarmStart=True, updatePrecond=True, - sub_problem_solver=EquationSolver.trust_region_minimize): + sub_problem_solver=EquationSolver_Immersed_2.trust_region_minimize): boundConstrainedObjective.reset_kappa() diff --git a/optimism/ConstrainedObjective.py b/optimism/ConstrainedObjective.py index 0e8cf11b..ea898403 100644 --- a/optimism/ConstrainedObjective.py +++ b/optimism/ConstrainedObjective.py @@ -1,5 +1,5 @@ from optimism.JaxConfig import * -import optimism.EquationSolver as EquationSolver +import EquationSolver_Immersed_2 as EquationSolver_Immersed_2 from optimism.Objective import Objective from optimism.Objective import param_index_update from optimism.SparseCholesky import SparseCholesky diff --git a/optimism/EquationSolver.py b/optimism/EquationSolver.py index 28f8574c..e9af15df 100644 --- a/optimism/EquationSolver.py +++ b/optimism/EquationSolver.py @@ -237,13 +237,14 @@ def print_min_banner(objective, modelObjective, res, modelRes, cgIters, trSize, willAccept = " True" if willAccept else "False" - print('obj=', f"{objective:15.9}", - ', model obj=', f"{modelObjective:15.9}", - ', res=', f"{res:15.9}", ', model res=', f"{modelRes:15.9}", - ', cg iters=', f"{cgIters:3}", - ', tr size=', f"{trSize:12.6}", - ', ', f"{onBoundary:9}", - ', accepted=', willAccept) + with open('./log_regsolve.txt', 'a') as f: + print('obj=', f"{objective:15.9}", + ', model obj=', f"{modelObjective:15.9}", + ', res=', f"{res:15.9}", ', model res=', f"{modelRes:15.9}", + ', cg iters=', f"{cgIters:3}", + ', tr size=', f"{trSize:12.6}", + ', ', f"{onBoundary:9}", + ', accepted=', willAccept) # utility function to save a matrix for later analysis offline @@ -426,7 +427,7 @@ def trust_region_minimize(objective, x, settings, callback=None): if is_converged(objective, y, realObjective, modelObjective, gy, g + Jd, cgIters, trSizeUsed, settings): if callback: callback(y, objective) - return y, True + return d, y, True modelImprove = -modelObjective realImprove = -realObjective @@ -584,7 +585,7 @@ def newton_solve(objective_func, solution, maxSteps=1): def nonlinear_equation_solve(objective, x0, p, settings, solver_algorithm=trust_region_minimize, callback=None, - useWarmStart=True, + useWarmStart=False, updatePrecond=True): xBar0 = objective.scaling * x0 diff --git a/optimism/EquationSolver_Immersed_2.py b/optimism/EquationSolver_Immersed_2.py new file mode 100644 index 00000000..ed4fe618 --- /dev/null +++ b/optimism/EquationSolver_Immersed_2.py @@ -0,0 +1,626 @@ +# Approach 2: Solver operates in BG, only objective computation in FG. + +from optimism.JaxConfig import * +import optimism.LU as LU +from optimism.Timer import Timer +from optimism import WarmStart +from scipy.sparse import save_npz +from scipy.sparse.linalg import LinearOperator, gmres + + +Settings = namedtuple('Settings', ['t1','t2','eta1','eta2','eta3','max_trust_iters','tol','max_cg_iters','max_cumulative_cg_iters','cg_tol','cg_inexact_solve_ratio','tr_size','min_tr_size','check_stability','use_preconditioned_inner_product_for_cg','use_incremental_objective','debug_info','over_iters']) + +#must have cg_tol < tol +# consider changing to 0.5, 1.75, 1e-9, 0.1, 0.75 +def get_settings(t1=0.25, t2=1.75, eta1=1e-10, eta2=0.1, eta3=0.5, + max_trust_iters=100, + tol=1e-8, + max_cg_iters=50, + max_cumulative_cg_iters=1000, + cg_tol=None, + cg_inexact_solve_ratio=1e-5, + tr_size=2.0, + min_tr_size=1e-8, + check_stability=False, + use_preconditioned_inner_product_for_cg=False, + use_incremental_objective=False, + debug_info=True, + over_iters=0): + if cg_tol==None: + cg_tol = 0.2*tol + + return Settings(t1, t2, eta1, eta2, eta3, + max_trust_iters=max_trust_iters, + tol=tol, + max_cg_iters=max_cg_iters, + max_cumulative_cg_iters=max_cumulative_cg_iters, + cg_tol=cg_tol, + cg_inexact_solve_ratio=cg_inexact_solve_ratio, + tr_size=tr_size, + min_tr_size=min_tr_size, + check_stability=check_stability, + use_preconditioned_inner_product_for_cg=use_preconditioned_inner_product_for_cg, + use_incremental_objective=use_incremental_objective, + debug_info=debug_info, + over_iters=over_iters) + + +def settings_with_new_tol(oldSettings, newTol): + newSettings = Settings(t1=oldSettings.t1, + t2=oldSettings.t2, + eta1=oldSettings.eta1, + eta2=oldSettings.eta2, + eta3=oldSettings.eta3, + max_trust_iters=oldSettings.max_trust_iters, + tol=newTol, # here + max_cg_iters=oldSettings.max_cg_iters, + max_cumulative_cg_iters=oldSettings.max_cumulative_cg_iters, + cg_tol=0.2*newTol, # here + cg_inexact_solve_ratio=oldSettings.cg_inexact_solve_ratio, + tr_size=oldSettings.tr_size, + min_tr_size=oldSettings.min_tr_size, + check_stability=oldSettings.check_stability, + use_preconditioned_inner_product_for_cg=oldSettings.use_preconditioned_inner_product_for_cg, + use_incremental_objective=oldSettings.use_incremental_objective, + debug_info=oldSettings.debug_info, + over_iters=oldSettings.over_iters) + + + return newSettings + + +negCurveString='neg curve' +boundaryString='boundary' +interiorString='interior' + + +def is_on_boundary(stepType): + return stepType==boundaryString or stepType==negCurveString + + +def project_to_boundary(z, d, trSize, zz): + # find tau s.t. (z + tau*d)^2 = trSize^2 + dd = np.dot(d,d) + zd = np.dot(z,d) + tau = (np.sqrt( (trSize**2-zz)*dd + zd**2 ) - zd)/dd + # is it ever better to choose the - sqrt() branch? + return z + tau*d + + +def preconditioned_project_to_boundary(z, d, T, trSize, zz, mult_by_approx_hessian): + Tt = np.transpose(T) + # find tau s.t. (z + tau*d)^2 = trSize^2 + Pd = (mult_by_approx_hessian(d)) # in BG + dd = np.dot(d,Pd) + zd = np.dot(z,Pd) + tau = (np.sqrt( (trSize**2-zz)*dd + zd**2 ) - zd)/dd + # is it ever better to choose the - sqrt() branch? + return z + tau*d + + +def project_to_boundary_with_coefs(z, d, trSize, zz, zd, dd): + # find tau s.t. (z + tau*d)^2 = trSize^2 + tau = (np.sqrt( (trSize**2-zz)*dd + zd**2 ) - zd)/dd + return z + tau*d + + +def update_step_length_squared(alpha, zz, zd, dd): + return zz + 2*alpha*zd + alpha*alpha*dd + + +def cg_inner_products_preconditioned(alpha, beta, zd, dd, rPr, z, d): + # recurrence formulas from Gould et al. doi:10.1137/S1052623497322735 + zd = beta * ( zd + alpha*dd ) + dd = rPr + beta*beta*dd + return zd, dd + + +def cg_inner_products_unpreconditioned(alpha, beta, zd, dd, rPr, z, d): + zd = z @ d + dd = d @ d + return zd, dd + + +# returns increment of solution and True if the entire system is solved +def solve_trust_region_equality_constraint(x, g, J, trSize, settings): + # minimuze 0.5*(J*z+g)^2 + + r = J.multiply_by_transpose( J.solve_transpose(J.solve(g)) ) + d = -r + z = 0.*x + zz = 0. + rr = np.dot(r,r) + + cgTolSquared = max(settings.cg_tol**2, 1e-10*rr) + + for i in range(settings.max_cg_iters): + curvature = np.sum( (J.solve(J@d))**2 ) + + if curvature <= 0: + return z, negCurveString, i + + alpha = rr / curvature + + zNp1 = z + alpha*d + zzNp1 = np.dot(zNp1,zNp1) + if zzNp1 > trSize**2: + return project_to_boundary(z, d, trSize, zz), boundaryString, i+1 + + z = zNp1 + zz = zzNp1 + + r += alpha * J.multiply_by_transpose(J.solve_transpose(J.solve(J@d))) + rrNp1 = np.dot(r,r) + + if rrNp1 < cgTolSquared: + return z, interiorString, i+1 + + beta = rrNp1 / rr + rr = rrNp1 + d = -r + beta*d + + return z, interiorString+'_', i+1 + + +def solve_trust_region_minimization(x, r, T, hess_vec_func, precond, trSize, settings): + # minimize r@z + 0.5*z@J@z + z = 0.*x # z is size of BG x + zz = 0. + Tt = np.transpose(T) + + cgInexactRelTol = settings.cg_inexact_solve_ratio + cgTolSquared = max(settings.cg_tol**2, cgInexactRelTol*cgInexactRelTol*r@r) + if np.dot(T,r)@np.dot(T,r) < cgTolSquared: + return z, z, interiorString, 0 + + Pr = (precond(r)) # BG + + d = -Pr # BG + cauchyP = np.array(d) #BG. + + rPr = r@Pr # BG + + zz = 0.0 + zd = 0.0 + if settings.use_preconditioned_inner_product_for_cg: + dd = rPr # BG + cg_inner_products = cg_inner_products_preconditioned + else: + dd = d @ d # BG + cg_inner_products = cg_inner_products_unpreconditioned + + for i in range(settings.max_cg_iters): + curvature = d@( np.dot(Tt,hess_vec_func(np.dot(T,d)) )) # BG. Map hess-vec to BG. + alpha = rPr / curvature # BG + + + + zNp1 = z + alpha*d # BG + zzNp1 = update_step_length_squared(alpha, zz, zd, dd) # All components are in BG. + + if curvature <= 0: + zOut = project_to_boundary_with_coefs(z, d, trSize, + zz, zd, dd) + return zOut, cauchyP, negCurveString, i+1 # ALL in BG. + + if zzNp1 > trSize**2: + zOut = project_to_boundary_with_coefs(z, d, trSize, + zz, zd, dd) # All in BG. + return zOut, cauchyP, boundaryString, i+1 + + z = zNp1 # z + alpha*d # BG. + + r += alpha * np.dot(Tt,hess_vec_func(np.dot(T,d))) # Map hess-vec to BG + Pr = precond(r) # BG + rPrNp1 = r@Pr # BG + + if r@r < cgTolSquared: + return z, cauchyP, interiorString, i+1 + + beta = rPrNp1 / rPr # BG + rPr = rPrNp1 #BG + d = -Pr + beta*d #BG + + zz = zzNp1 + zd, dd = cg_inner_products(alpha, beta, zd, dd, rPr, z, d) # everything in BG. + + return z, cauchyP, interiorString+'_', i+1 + + +# essentially deprecated +def print_banner(objective, modelObjective, cgIters, trSize, onBoundary, willAccept): + willAccept = " True" if willAccept else "False" + + print('obj=', f"{objective:16.11}", + ', model obj=', f"{modelObjective:16.11}", + ', cg iters=', f"{cgIters:3}", + ', tr size=', f"{trSize:12.6}", + ', ', f"{onBoundary:9}", + ', acceped=', willAccept) + + +def print_min_banner(objective, modelObjective, res, modelRes, cgIters, trSize, onBoundary, willAccept, settings): + if settings.debug_info==False: return + + willAccept = " True" if willAccept else "False" + print('obj=', f"{objective:15.9}", + ', model obj=', f"{modelObjective:15.9}", + ', res=', f"{res:15.9}", ', model res=', f"{modelRes:15.9}", + ', cg iters=', f"{cgIters:3}", + ', tr size=', f"{trSize:12.6}", + ', ', f"{onBoundary:9}", + ', accepted=', willAccept) + + +# utility function to save a matrix for later analysis offline +def output_matrix(objective, x): + hess = objective.hessian(x) + hess = SparseCholesky.convert_to_sparse(hess) + print('saving matrix') + save_npz('matrix', hess) + exit() + + +def trust_region_least_squares_solve(objective, x, settings): + trSize = settings.tr_size + + equation = objective.gradient + g = equation(x) + o = 0.5*np.dot(g, g) + + J = LU.LU(objective.hessian(x)) + + for i in range(settings.max_trust_iters): + + dx, stepType, cgIters = solve_trust_region_equality_constraint(x, g, J, trSize, settings) + + modelObjective = 0.5*np.sum( (J@dx + g)**2 ) + + y = x+dx + gy = equation(y) + oy = 0.5*np.dot(gy, gy) + + if 2*oy < settings.tol**2: + print_banner(oy, modelObjective, cgIters, trSize, stepType, True) + if settings.check_stability: + objective.check_stability(y) + return y, True + + modelImprove = o - modelObjective + realImprove = o - oy + + rho = realImprove / modelImprove + + if not (rho >= settings.eta2): + trSize *= settings.t1 + elif rho > settings.eta3 and is_on_boundary(stepType): + trSize *= settings.t2 + + print_banner(oy, modelObjective, cgIters, trSize, stepType, rho >= settings.eta1) + + if trSize < settings.min_tr_size: break + + if rho >= settings.eta1: + x = y + g = gy + o = oy + J.update( objective.hessian(x) ) + + print("Reached the maximum number of trust region iterations or trust region is too small.") + if settings.check_stability: + objective.check_stability(x) + return x, False + + +def is_converged(objective, x, realO, modelO, realRes, modelRes, cgIters, trSize, settings): + gg = realRes@realRes + if gg < settings.tol**2: + modelResNorm = np.linalg.norm(modelRes) + realResNorm = np.sqrt(gg) + print_min_banner(realO, modelO, + realResNorm, + modelResNorm, + cgIters, + trSize, + interiorString, + True, + settings) + if settings.check_stability: + objective.check_stability(x) + + print('') # a bit of output formatting + + return True + return False + + +def dogleg_step(cp, newtonP, T, trSize, mat_mul): + Tt = np.transpose(T) + cc = cp@mat_mul(cp) # Precond is FG, input also FG, so map input to FG and output to BG + nn = newtonP@mat_mul(newtonP) # Precond is FG, input also FG, so map input to FG and output to BG + tt = trSize*trSize + + if cc >= tt: #return cauchy point if it extends outside the tr + #print('cp on boundary') + return cp * np.sqrt(tt/cc) + + if cc > nn: # return cauchy point? seems the preconditioner was not accurate? + print('cp outside newton, preconditioner likely inaccurate') + return cp + + if nn > tt: # on the dogleg (we have nn >= cc, and tt >= cc) + #print('dogleg') + return preconditioned_project_to_boundary(cp, + newtonP-cp, T, + trSize, + cc, + mat_mul) + #print('quasi-newton step') + return newtonP + + +def trust_region_minimize(objective, x, T, settings, callback=None): + # All in background. + trSize = settings.tr_size + triedNewPrecond = False + Tt = np.transpose(T) + gradientAndTanOpt = objective.gradient_and_tangent + gradient = objective.gradient + print(np.linalg.norm(np.dot(T,x))) + + #g,hess_vec_func = gradientAndTanOpt(x) + g = np.dot(Tt,gradient(np.dot(T,x))) # Map to FG for computing gradient then map back to BG + print(np.linalg.norm(g)) + o = objective.value(np.dot(T,x)) # Obj comp on foreground. + gNorm = np.linalg.norm(g) #BG. + print("\nInitial objective, residual = ", o, gNorm) + + # this could potentially return an unstable solution + if is_converged(objective, x, 0.0, 0.0, g, g, 0, trSize, settings): + if callback: callback(x, objective) + return x, True + + cumulativeCgIters=0 + + for i in range(settings.max_trust_iters): + # minimize 0.5*(2*r + J_sd)*d = r + 0.5*dJd + + if settings.use_incremental_objective: + incremental_objective = lambda d: 0.5*((g + objective.gradient(x+d)) @ d) # False. Does not matter + else: + incremental_objective = lambda d: objective.value(np.dot(T,x+d)) - o # Map to FG for obj. + + hess_vec_func = lambda v: objective.hessian_vec(np.dot(T,x), v) # Map to FG since hessian comp there. + mult_by_approx_hessian = objective.multiply_by_approx_hessian if settings.use_preconditioned_inner_product_for_cg else lambda x: x + + gKg = g@np.dot(Tt,hess_vec_func(np.dot(T,g))) # Input conv. in fn definition, hess-vec in FG so map it. + print(np.isnan(gKg)) + if gKg > 0: + alpha = -(g@g) / gKg # In BG + cauchyPoint = alpha * g # In BG + cauchyPointNormSquared = cauchyPoint@(mult_by_approx_hessian(cauchyPoint)) # BG + else: + cauchyPoint = -g * (trSize / np.sqrt((g@(Tt,mult_by_approx_hessian(g))))) # BG + cauchyPointNormSquared = trSize*trSize + print('negative curvature unpreconditioned cauchy point direction found.') + + if cauchyPointNormSquared >= trSize*trSize: + print('unpreconditioned gradient cauchy point outside trust region at dist = ', np.sqrt(cauchyPointNormSquared)) + cauchyPoint *= (trSize / np.sqrt(cauchyPointNormSquared)) + cauchyPointNormSquared = trSize*trSize + qNewtonPoint = cauchyPoint # BG + stepType = boundaryString + cgIters = 1 + else: + qNewtonPointb, _, stepType, cgIters = \ + solve_trust_region_minimization(x, g, T, + hess_vec_func, + objective.apply_precond, + trSize, settings) # All in BG + qNewtonPoint = qNewtonPointb + + + + cumulativeCgIters += cgIters + + trSizeUsed = trSize + happyAboutTrSize=False + while not happyAboutTrSize: + d = dogleg_step(cauchyPoint, qNewtonPoint, T, trSize, mult_by_approx_hessian) # d in BG + + Jd = np.dot(Tt,hess_vec_func(np.dot(T,d))) #BG + dJd = d @ Jd #BG + modelObjective = g @ d + 0.5*dJd + + y = x+d + realObjective = incremental_objective(d) # Fn definition maps it to FG. + gy = np.dot(Tt,gradient(np.dot(T,y))) + + if is_converged(objective, y, realObjective, modelObjective, + gy, g + Jd, cgIters, trSizeUsed, settings): + if callback: callback(y, objective) # Convergence check on foreground. + return d, y, True + + modelImprove = -modelObjective + realImprove = -realObjective + + rho = realImprove / modelImprove + + if modelObjective > 0: + print('Found a positive model objective increase. Debug if you see this.') + rho = realImprove / -modelImprove + #exit(1) + + if not rho >= settings.eta2: # write it this way to handle NaNs + trSize *= settings.t1 + elif rho > settings.eta3 and is_on_boundary(stepType): + trSize *= settings.t2 + + modelRes = g + Jd + modelResNorm = np.linalg.norm(modelRes) + realResNorm = np.linalg.norm(gy) + + willAccept = rho >= settings.eta1 or (rho >= -0 and realResNorm <= gNorm) + + print_min_banner(realObjective, modelObjective, + realResNorm, modelResNorm, + cgIters, trSizeUsed, stepType, + willAccept, + settings) + + if willAccept: + x = y # BG + g = gy # BG + #g,hess_vec_func = gradientAndTanOpt(x) + o = objective.value(np.dot(T,x)) # Map input to FG + gNorm = realResNorm #BG + triedNewPrecond = False + happyAboutTrSize = True + + if callback: callback(x, objective) + else: + # set these for output + # trust region will continue to strink until we find a solution on the boundary + stepType=boundaryString + cgIters = 0 + + if cgIters >= settings.max_cg_iters or cumulativeCgIters >= settings.max_cumulative_cg_iters: + objective.update_precond(np.dot(T,x)) # Preconditioner needs FG input. + cumulativeCgIters=0 + + trSizeUsed = trSize + + if trSize < settings.min_tr_size: + + if not triedNewPrecond: + print("The trust region is too small, updating precond and trying again.") + objective.update_precond(np.dot(T,x)) # Preconditioner needs FG input. + cumulativeCgIters=0 + triedNewPrecond = True + happyAboutTrSize = True + trSize = settings.tr_size + else: + print("The trust region is still too small. Accepting, but be careful.") + if callback: callback(x, objective) + return d, x, False + + print("Reached the maximum number of trust region iterations.") + if settings.check_stability: + objective.check_stability(x) + + if callback: callback(x, objective) + return d, x, False + + +def newton(objective, x, settings, callback=None): + + res_func = objective.gradient + hess_func = objective.hessian + + sz = x.size + precondOp = LinearOperator((sz,sz), + lambda v: objective.apply_precond(v)) + + for i in range(settings.max_trust_iters): + print('Newton step = ', i) + + if callback: callback(x, objective.p) + + sparse=True + if sparse: + def apply_a(v): + return objective.hessian_vec(x, v) + + aOp = LinearOperator((sz,sz), + lambda v: apply_a(v)) + + g = res_func(x) + gNorm = np.linalg.norm(g) + + print('norm = ', gNorm) + + dx, exitcode = gmres(aOp, g, tol=1e-3*gNorm, atol=0, M=precondOp, maxiter=settings.max_cg_iters) + + gTrial = np.linalg.norm( res_func(x-dx) ) + print('trial norm = ', gTrial) + while gTrial > gNorm: + dx *= 0.25 + gTrial = np.linalg.norm( res_func(x-dx) ) + print('trial norm = ', gTrial) + + x -= dx + + else: + g = res_func(x) + H = hess_func(x) + x -= np.linalg.solve(H,g) + + o = objective.value(x) + if is_converged(objective, x, o, o, + g, g, 1, np.inf, settings): + return x, True + + print("Reached the maximum number of newton iterations.") + if settings.check_stability: + objective.check_stability(x) + return x, False + + +# This solve uses a different interface. There is no globalization yet + +def newton_solve(objective_func, solution, maxSteps=1): + grad_func = grad(objective_func) + hess_func = jacfwd(grad_func) + + for step in range(maxSteps): + with Timer(name="Step"): + with Timer(name="Fill Grad"): + gradients = grad_func(solution) + print('step: ', step, 'res = ', np.linalg.norm(gradients)) + + with Timer(name="Fill Hess"): + hess = hess_func(solution) + + nDofs = gradients.size + gradShape = gradients.shape + + gradients = np.reshape(gradients, (nDofs,)) + + hess = np.reshape(hess, (nDofs, nDofs)) + + with Timer(name="LinearSolve"): + solution = solution - np.reshape(np.linalg.solve(hess,gradients), gradShape) + + return solution, True + + +def nonlinear_equation_solve(objective, x0, p, T, settings, + solver_algorithm=trust_region_minimize, + callback=None, + useWarmStart=False, + updatePrecond=True): + Tt = np.transpose(T) + xBar0 = objective.scaling * x0 + + + if useWarmStart: + if updatePrecond: + objective.update_precond(xBar0) + + dxBar = WarmStart.warm_start_increment(objective, + xBar0, p) + xBar0 += dxBar + objective.p = p + else: + objective.p = p + xBar0 = xBar0 + + if updatePrecond: + + objective.update_precond(np.dot(T,xBar0)) # Map to FG for preconditioner. + print('preconditioner updated') + + d, xBar, solverSuccess = solver_algorithm(objective, xBar0, T, settings, callback=callback) # All in background, except the obj part. + + return objective.invScaling *xBar, solverSuccess, np.divide(xBar,objective.invScaling *xBar) + \ No newline at end of file diff --git a/optimism/FunctionSpace.py b/optimism/FunctionSpace.py index 4903de20..fc7e4607 100644 --- a/optimism/FunctionSpace.py +++ b/optimism/FunctionSpace.py @@ -337,15 +337,16 @@ def interpolate_nodal_field_on_edge(functionSpace, U, interpolationPoints, edge) """ edgeShapes = Interpolants.compute_shapes(functionSpace.mesh.parentElement1d, interpolationPoints) edgeU = get_nodal_values_on_edge(functionSpace, U, edge) - return edgeShapes.values.T@edgeU + return edgeShapes.values.T@edgeU, edgeShapes.gradients.T@edgeU def integrate_function_on_edge(functionSpace, func, U, quadRule, edge): - uq = interpolate_nodal_field_on_edge(functionSpace, U, quadRule.xigauss, edge) + uq, uq_eta = interpolate_nodal_field_on_edge(functionSpace, U, quadRule.xigauss, edge) Xq = interpolate_nodal_field_on_edge(functionSpace, functionSpace.mesh.coords, quadRule.xigauss, edge) edgeCoords = Mesh.get_edge_coords(functionSpace.mesh, edge) - _, normal, jac = Mesh.compute_edge_vectors(functionSpace.mesh, edgeCoords) - integrand = jax.vmap(func, (0, 0, None))(uq, Xq, normal) + _, normal, jac, normals = Mesh.compute_edge_vectors(functionSpace.mesh, edgeCoords) + uq_x = np.array([np.multiply(uq_eta[0,:],normals[0]),np.multiply(uq_eta[1,:],normals[1])]) + integrand = jax.vmap(func, (0, 0, None, None))(uq, Xq, normal, uq_x) return np.dot(integrand, jac*quadRule.wgauss) diff --git a/optimism/Mechanics.py b/optimism/Mechanics.py index 263fa1da..e13f8ee9 100644 --- a/optimism/Mechanics.py +++ b/optimism/Mechanics.py @@ -7,6 +7,8 @@ from optimism.TensorMath import tensor_2D_to_3D from optimism import QuadratureRule from optimism import Interpolants +import jax.numpy as np + MechanicsFunctions = namedtuple('MechanicsFunctions', ['compute_strain_energy', @@ -466,7 +468,57 @@ def compute_traction_potential_energy(fs, U, quadRule, edges, load): outward unit normal. time: current time """ - def compute_energy_density(u, X, n): + def compute_energy_density(u, X, n, f): + f = 1 traction = load(X, n) return -np.dot(u, traction) return FunctionSpace.integrate_function_on_edges(fs, compute_energy_density, U, quadRule, edges) + +# Penalty Function +def compute_displacement_penalty(fs, U, quadRule, edge, dispTarget, comp, penalty): + """Compute penalty energy. + + Arguments: + fs: a FunctionSpace object + U: the nodal displacements + quadRule: the 1D quadrature rule to use for the integration + edges: array of edges, each row is an edge. Each edge has two entries, the + element ID, and the permutation of that edge in the triangle (0, 1, + 2). + load: Callable that returns the traction vector. The signature is + load(X, n), where X is coordinates of a material point, and n is the + outward unit normal. + time: current time + """ + def compute_energy_density(u, X, n, f): + f = 1 + traction = dispTarget(X, n) + if comp == 2: + return penalty*((u[1]-(traction)[1])**2 + (u[0]-(traction)[0])**2) + elif comp == 1: + return penalty*(u[1]-(traction)[1])**2 + elif comp == 0: + return penalty*(u[0]-(traction)[0])**2 + return FunctionSpace.integrate_function_on_edges(fs, compute_energy_density, U, quadRule, edge) +# THIS IS NOT COMPLETELY CORRECT AND IS NOT BEING USED ANYWHERE +def compute_nitsche_term(fs, U, quadRule, edge): + """Compute Nitsche. IMPORTANT - ONLY VALID FOR LINEAR ELASTICITY. + + Arguments: + fs: a FunctionSpace object + U: the nodal displacements + quadRule: the 1D quadrature rule to use for the integration + edges: array of edges, each row is an edge. Each edge has two entries, the + element ID, and the permutation of that edge in the triangle (0, 1, + 2). + load: Callable that returns the traction vector. The signature is + load(X, n), where X is coordinates of a material point, and n is the + outward unit normal. + time: current time + """ + def compute_energy_density_nitsche(u, X, n, u_x): + # Compute displacement-normal product. + gradUn = u_x@n + return np.dot(gradUn,u) + + return FunctionSpace.integrate_function_on_edges(fs, compute_energy_density_nitsche, U, quadRule, edge) diff --git a/optimism/Mesh.py b/optimism/Mesh.py index 019d2493..2a2c2747 100644 --- a/optimism/Mesh.py +++ b/optimism/Mesh.py @@ -303,7 +303,7 @@ def num_elements(mesh): def num_nodes(mesh): - return mesh.coords.shape[0] + return np.shape(np.reshape(np.unique(mesh.conns),-1))[0]#mesh.coords.shape[0] def get_edge_node_indices(mesh : Mesh, edge): @@ -347,4 +347,4 @@ def compute_edge_vectors(mesh : Mesh, edgeCoords): tangent = Xv[1] - Xv[0] normal = np.array([tangent[1], -tangent[0]]) jac = np.linalg.norm(tangent) - return tangent/jac, normal/jac, jac \ No newline at end of file + return tangent/jac, normal/jac, jac, normal \ No newline at end of file diff --git a/optimism/Objective.py b/optimism/Objective.py index 8f3d3f8b..0ea13cfa 100644 --- a/optimism/Objective.py +++ b/optimism/Objective.py @@ -77,7 +77,11 @@ def precond_at_attempt(self, attempt): class Objective: - def __init__(self, f, x, p, precondStrategy=None): + def __init__(self, f, x, p, T, precondStrategy=None): + + # Define the transposed extraction operator here - + self.Tt = np.transpose(T) + self.Tm = T self.precond = SparseCholesky() self.precondStrategy = precondStrategy @@ -85,14 +89,15 @@ def __init__(self, f, x, p, precondStrategy=None): self.p = p self.objective=jit(f) - self.grad_x = jit(grad(f,0)) + self.grad_x = jit((grad(f,0))) self.grad_p = jit(grad(f,1)) self.hess_vec = jit(lambda x, p, vx: jvp(lambda z: self.grad_x(z,p), (x,), (vx,))[1]) + self.vec_hess = jit(lambda x, p, vx: - vjp(lambda z: self.grad_x(z,p), x)[1](vx)) + (vjp(lambda z: self.grad_x(z,p), x)[1](vx))) self.jac_xp_vec = jit(lambda x, p, vp0: jvp(lambda q0: self.grad_x(x, param_index_update(p,0,q0)), @@ -116,10 +121,10 @@ def __init__(self, f, x, p, precondStrategy=None): self.vec_jac_xp4 = jit(lambda x, p, vx: vjp(lambda q4: self.grad_x(x, param_index_update(p,4,q4)), p[4])[1](vx)) - + # Doubt here. self.grad_and_tangent = lambda x, p: linearize(lambda z: self.grad_x(z,p), x) - self.hess = jit(jacfwd(self.grad_x, 0)) + self.hess = jit((jacfwd(self.grad_x, 0))) self.scaling = 1.0 self.invScaling = 1.0 @@ -129,19 +134,19 @@ def value(self, x): return self.objective(x, self.p) def gradient(self, x): - return self.grad_x(x, self.p) + return (self.grad_x(x, self.p)) def gradient_p(self, x): return self.grad_p(x, self.p) def hessian_vec(self, x, vx): - return self.hess_vec(x, self.p, vx) + return (self.hess_vec(x, self.p, vx)) def gradient_and_tangent(self, x): - return self.grad_and_tangent(x, self.p) + return (self.grad_and_tangent(x, self.p)) def vec_hessian(self, x, vx): - return self.vec_hess(x, self.p, vx) + return (self.vec_hess(x, self.p, vx)) def hessian(self, x): return self.hess(x, self.p) diff --git a/optimism/ReadExodusMesh.py b/optimism/ReadExodusMesh.py index de6fc27e..16891425 100644 --- a/optimism/ReadExodusMesh.py +++ b/optimism/ReadExodusMesh.py @@ -7,12 +7,36 @@ exodusToNativeTri6NodeOrder = np.array([0, 3, 1, 5, 4, 2]) -def read_exodus_mesh(fileName): +def read__mesh(fileName): with netCDF4.Dataset(fileName) as exData: coords = _read_coordinates(exData) conns, blocks = _read_blocks(exData) nodeSets = _read_node_sets(exData) sideSets = _read_side_sets(exData) + connsn = onp.array(conns) + + + # Rearranging arrays; we are only reading one block so stuff should only come from there. + # Get foreground Node indices first. + FG_Node_Indices = onp.transpose(onp.unique(onp.reshape(connsn,-1))) + + # Now, get the coordinates corresponding to that + coords_FG = onp.array(coords[FG_Node_Indices,:]) + + + + # Now, overwrite conns, because that still uses the old indices + for i in range(onp.shape(connsn)[0]): + for j in range(onp.shape(connsn)[1]): + val = connsn[i,j] + ind = onp.where(FG_Node_Indices == val) + indv = ind[0].item() + connsn[i,j] = indv + + conns = np.array(connsn) + + # Finally, overwrite coords, so that all the other stuff can come from it. + coords = np.array(coords_FG) elementType = _read_element_type(exData).lower() if elementType == "tri3" or elementType == "tri": @@ -22,9 +46,10 @@ def read_exodus_mesh(fileName): basis, basis1d = Interpolants.make_parent_elements(degree = 2) simplexNodesOrdinals = _get_vertex_nodes_from_exodus_tri6_mesh(conns) conns = conns[:, exodusToNativeTri6NodeOrder] - else: - raise - + #else: + # raise Exception('Cannot work with higher than 2nd order elements.') + + return Mesh.Mesh(coords, conns, simplexNodesOrdinals, basis, basis1d, blocks, nodeSets, sideSets) @@ -75,7 +100,7 @@ def _read_blocks(exodusDataset): blockConns = [] blocks = {} firstElemInBlock = 0 - for i in range(nBlocks): + for i in range(1): nodesPerElemInBlock = len(exodusDataset.dimensions['num_nod_per_el' + str(i+1)]) assert nodesPerElemInBlock == nodesPerElem diff --git a/optimism/SparseCholesky.py b/optimism/SparseCholesky.py index e6e43584..09e215b9 100644 --- a/optimism/SparseCholesky.py +++ b/optimism/SparseCholesky.py @@ -12,7 +12,7 @@ class SparseCholesky: def factorize(self): - + print('Assembling preconditioner', 0) stiffnessTryStep = 0 self.A = self.new_stiffness_func(stiffnessTryStep) diff --git a/optimism/WarmStart.py b/optimism/WarmStart.py index 47642bf1..441e8c1e 100644 --- a/optimism/WarmStart.py +++ b/optimism/WarmStart.py @@ -3,7 +3,7 @@ from optimism.JaxConfig import * import optimism.Objective as Objective -def warm_start_increment(objective, x, pNew, index=0): +def warm_start_increment(objective, x, T, pNew, index=0): dp = objective.p[index] - pNew[index] if index==0: diff --git "a/optimism/__init__.py\357\200\272Zone.Identifier" "b/optimism/__init__.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git "a/optimism/contact/Contact.py\357\200\272Zone.Identifier" "b/optimism/contact/Contact.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git "a/optimism/contact/EdgeCpp.py\357\200\272Zone.Identifier" "b/optimism/contact/EdgeCpp.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git "a/optimism/contact/EdgeIntersection.py\357\200\272Zone.Identifier" "b/optimism/contact/EdgeIntersection.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git "a/optimism/contact/Friction.py\357\200\272Zone.Identifier" "b/optimism/contact/Friction.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git "a/optimism/contact/Levelset.py\357\200\272Zone.Identifier" "b/optimism/contact/Levelset.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git "a/optimism/contact/LevelsetConstraint.py\357\200\272Zone.Identifier" "b/optimism/contact/LevelsetConstraint.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git "a/optimism/contact/MortarContact.py\357\200\272Zone.Identifier" "b/optimism/contact/MortarContact.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git "a/optimism/contact/PenaltyContact.py\357\200\272Zone.Identifier" "b/optimism/contact/PenaltyContact.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git "a/optimism/contact/Search.py\357\200\272Zone.Identifier" "b/optimism/contact/Search.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git "a/optimism/contact/__init__.py\357\200\272Zone.Identifier" "b/optimism/contact/__init__.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git "a/optimism/inverse/AdjointFunctionSpace.py\357\200\272Zone.Identifier" "b/optimism/inverse/AdjointFunctionSpace.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git "a/optimism/inverse/DesignOpt.py\357\200\272Zone.Identifier" "b/optimism/inverse/DesignOpt.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git "a/optimism/inverse/MechanicsInverse.py\357\200\272Zone.Identifier" "b/optimism/inverse/MechanicsInverse.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git a/optimism/inverse/NonlinearSolve.py b/optimism/inverse/NonlinearSolve.py index 37b26ccf..65c2ca2b 100644 --- a/optimism/inverse/NonlinearSolve.py +++ b/optimism/inverse/NonlinearSolve.py @@ -1,14 +1,14 @@ from optimism.JaxConfig import * from jax import custom_jvp, custom_vjp -from optimism import EquationSolver +from . import EquationSolver_Immersed_2 from optimism import Objective from optimism import WarmStart @partial(custom_vjp, nondiff_argnums=(0,1)) def nonlinear_solve(mechanicalEnergy, settings, UuGuess, designParams): p = Objective.param_index_update(mechanicalEnergy.p, 2, designParams) - return EquationSolver.nonlinear_equation_solve(mechanicalEnergy, + return EquationSolver_Immersed_2.nonlinear_equation_solve(mechanicalEnergy, UuGuess, p, settings)[0] @@ -23,7 +23,7 @@ def nonlinear_solve_b(mechanicalEnergy, settings, rdata, v): hess_vec_func = lambda w: mechanicalEnergy.hessian_vec(Uu, w) - results,_ = EquationSolver.solve_trust_region_minimization(0.0*Uu, + results,_ = EquationSolver_Immersed_2.solve_trust_region_minimization(0.0*Uu, v, hess_vec_func, mechanicalEnergy.apply_precond, @@ -45,7 +45,7 @@ def nonlinear_solve_with_state(mechanicalEnergy, settings, UuGuess, p): mechanicalEnergy.update_precond(UuGuess) UuGuess += WarmStart.warm_start_increment_jax_safe(mechanicalEnergy, UuGuess, p[0]) - Uu, solverSuccess = EquationSolver.nonlinear_equation_solve(mechanicalEnergy, + Uu, solverSuccess = EquationSolver_Immersed_2.nonlinear_equation_solve(mechanicalEnergy, UuGuess, p, settings, @@ -65,7 +65,7 @@ def nonlinear_solve_with_state_b(mechanicalEnergy, settings, rdata, v): hess_vec_func = lambda w: mechanicalEnergy.hessian_vec(Uu, w) UuZeros = np.zeros_like(Uu) - results = EquationSolver.solve_trust_region_minimization(UuZeros, + results = EquationSolver_Immersed_2.solve_trust_region_minimization(UuZeros, v, hess_vec_func, mechanicalEnergy.apply_precond, diff --git "a/optimism/inverse/NonlinearSolve.py\357\200\272Zone.Identifier" "b/optimism/inverse/NonlinearSolve.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git "a/optimism/inverse/ShapeOpt.py\357\200\272Zone.Identifier" "b/optimism/inverse/ShapeOpt.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git "a/optimism/inverse/TopOpt.py\357\200\272Zone.Identifier" "b/optimism/inverse/TopOpt.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git "a/optimism/material/Gent.py\357\200\272Zone.Identifier" "b/optimism/material/Gent.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git "a/optimism/material/Hardening.py\357\200\272Zone.Identifier" "b/optimism/material/Hardening.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git a/optimism/material/HyperViscoelastic.py b/optimism/material/HyperViscoelastic.py index dea869f9..1dae91c6 100644 --- a/optimism/material/HyperViscoelastic.py +++ b/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_dissipated_energy(dispGrad, state, dt, props) + return _compute_dissipation(dispGrad, state, dt, props) return MaterialModel(compute_energy_density = energy_density, compute_initial_state = compute_initial_state, @@ -63,11 +63,9 @@ def _energy_density(dispGrad, state, dt, props): Ee = Ee_trial - delta_Ev W_neq = _neq_strain_energy(Ee, props) - - Dv = delta_Ev / dt - Psi = _dissipation_potential(Dv, props) + Psi = _incremental_dissipated_energy(Ee, dt, props) - return W_eq + W_neq + dt * Psi + return W_eq + W_neq + Psi def _eq_strain_energy(dispGrad, props): K, G = props[PROPS_K_eq], props[PROPS_G_eq] @@ -83,19 +81,22 @@ def _neq_strain_energy(elasticStrain, props): G_neq = props[PROPS_G_neq] return G_neq * TensorMath.norm_of_deviator_squared(elasticStrain) -def _dissipation_potential(Dv, props): +def _incremental_dissipated_energy(elasticStrain, dt, props): G_neq = props[PROPS_G_neq] tau = props[PROPS_TAU] eta = G_neq * tau - return eta * TensorMath.norm_of_deviator_squared(Dv) + Me = 2. * G_neq * elasticStrain + M_bar = TensorMath.norm_of_deviator_squared(Me) + + return 0.5 * dt * M_bar / eta -def _compute_dissipated_energy(dispGrad, state, dt, props): +def _compute_dissipation(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 - Dv = delta_Ev / dt - return dt * _dissipation_potential(Dv, props) + return _incremental_dissipated_energy(Ee, dt, props) def _compute_state_new(dispGrad, stateOld, dt, props): Ee_trial = _compute_elastic_logarithmic_strain(dispGrad, stateOld) diff --git "a/optimism/material/HyperViscoelastic.py\357\200\272Zone.Identifier" "b/optimism/material/HyperViscoelastic.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git "a/optimism/material/J2Plastic.py\357\200\272Zone.Identifier" "b/optimism/material/J2Plastic.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git "a/optimism/material/LinearElastic.py\357\200\272Zone.Identifier" "b/optimism/material/LinearElastic.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git "a/optimism/material/MaterialModel.py\357\200\272Zone.Identifier" "b/optimism/material/MaterialModel.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git a/optimism/material/MaterialUniaxialSimulator.py b/optimism/material/MaterialUniaxialSimulator.py index 787b7081..110eb209 100644 --- a/optimism/material/MaterialUniaxialSimulator.py +++ b/optimism/material/MaterialUniaxialSimulator.py @@ -4,7 +4,7 @@ import jax import jax.numpy as np -from optimism import EquationSolver as EqSolver +from . import EquationSolver_Immersed_2 as EqSolver from optimism import Objective from optimism import TensorMath diff --git "a/optimism/material/MaterialUniaxialSimulator.py\357\200\272Zone.Identifier" "b/optimism/material/MaterialUniaxialSimulator.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git a/optimism/material/Neohookean.py b/optimism/material/Neohookean.py index cfb5ec76..4055cd60 100644 --- a/optimism/material/Neohookean.py +++ b/optimism/material/Neohookean.py @@ -1,5 +1,5 @@ import jax.numpy as np - +import jax from optimism.material.MaterialModel import MaterialModel # props @@ -49,6 +49,7 @@ def _neohookean_3D_energy_density(dispGrad, internalVariables, props): F = dispGrad + np.eye(3) J = np.linalg.det(F) + #Wvol = 0.125*props[PROPS_LAMBDA]*(J - 1.0/J)**2 Wvol = 0.5*props[PROPS_LAMBDA]*np.log(J)**2 diff --git "a/optimism/material/Neohookean.py\357\200\272Zone.Identifier" "b/optimism/material/Neohookean.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git "a/optimism/material/__init__.py\357\200\272Zone.Identifier" "b/optimism/material/__init__.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git a/optimism/phasefield/MaterialPointSimulator.py b/optimism/phasefield/MaterialPointSimulator.py index 3bb0e244..c9f78b0c 100644 --- a/optimism/phasefield/MaterialPointSimulator.py +++ b/optimism/phasefield/MaterialPointSimulator.py @@ -3,7 +3,7 @@ from jax import disable_jit from optimism.JaxConfig import * -from optimism import EquationSolver as EqSolver +from . import EquationSolver_Immersed_2 as EqSolver from optimism import AlSolver from optimism import Objective from optimism import ConstrainedObjective diff --git "a/optimism/phasefield/MaterialPointSimulator.py\357\200\272Zone.Identifier" "b/optimism/phasefield/MaterialPointSimulator.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git "a/optimism/phasefield/PhaseField.py\357\200\272Zone.Identifier" "b/optimism/phasefield/PhaseField.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git "a/optimism/phasefield/PhaseFieldClassic.py\357\200\272Zone.Identifier" "b/optimism/phasefield/PhaseFieldClassic.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git "a/optimism/phasefield/PhaseFieldLorentzPlastic.py\357\200\272Zone.Identifier" "b/optimism/phasefield/PhaseFieldLorentzPlastic.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git "a/optimism/phasefield/PhaseFieldMaterialModel.py\357\200\272Zone.Identifier" "b/optimism/phasefield/PhaseFieldMaterialModel.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git "a/optimism/phasefield/PhaseFieldThreshold.py\357\200\272Zone.Identifier" "b/optimism/phasefield/PhaseFieldThreshold.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git "a/optimism/phasefield/PhaseFieldThresholdPlastic.py\357\200\272Zone.Identifier" "b/optimism/phasefield/PhaseFieldThresholdPlastic.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git "a/optimism/phasefield/__init__.py\357\200\272Zone.Identifier" "b/optimism/phasefield/__init__.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git "a/optimism/test/.gitignore\357\200\272Zone.Identifier" "b/optimism/test/.gitignore\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git "a/optimism/test/ConvexityPlot.py\357\200\272Zone.Identifier" "b/optimism/test/ConvexityPlot.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git "a/optimism/test/MeshFixture.py\357\200\272Zone.Identifier" "b/optimism/test/MeshFixture.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git "a/optimism/test/TestFixture.py\357\200\272Zone.Identifier" "b/optimism/test/TestFixture.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git "a/optimism/test/__init__.py\357\200\272Zone.Identifier" "b/optimism/test/__init__.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git a/optimism/test/contact/.gitignore b/optimism/test/contact/.gitignore new file mode 100644 index 00000000..0a680a42 --- /dev/null +++ b/optimism/test/contact/.gitignore @@ -0,0 +1,5 @@ +*~ +*.png +*.pyc +*.vtk +__py*__ diff --git "a/optimism/test/contact/.gitignore\357\200\272Zone.Identifier" "b/optimism/test/contact/.gitignore\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git a/optimism/test/contact/__init__.py b/optimism/test/contact/__init__.py new file mode 100644 index 00000000..e69de29b diff --git "a/optimism/test/contact/__init__.py\357\200\272Zone.Identifier" "b/optimism/test/contact/__init__.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git a/optimism/test/contact/test_Contact.py b/optimism/test/contact/test_Contact.py new file mode 100644 index 00000000..94aae40e --- /dev/null +++ b/optimism/test/contact/test_Contact.py @@ -0,0 +1,90 @@ +import unittest +from optimism.JaxConfig import * +from ..MeshFixture import MeshFixture +from optimism.contact import Contact +from optimism.contact import Friction +from optimism import Mesh +from optimism import QuadratureRule + +import itertools, operator + +def sort_uniq(sequence): + return np.array(list(map( + operator.itemgetter(0), + itertools.groupby(sorted(sequence))))) + + +frictionParams = Friction.Params(0.3, 1e-4) + +class TestContactFrictionData(MeshFixture): + + def setUp(self): + targetDispGrad = np.array([[0.1, -0.2],[0.4, -0.1]]) + #targetDispGrad = np.zeros((2,2)) + + m1 = self.create_mesh_and_disp(5, 3, [0., 1.], [0., 1.], + lambda x : targetDispGrad.dot(x), + setNamePostFix='1') + + m2 = self.create_mesh_and_disp(6, 2, [0., 1.], [1., 2.], + lambda x : targetDispGrad.dot(x), + setNamePostFix='2') + + self.mesh, self.U = Mesh.combine_mesh(m1, m2) + + self.surfM = self.mesh.sideSets['top1'] + self.surfI = self.mesh.sideSets['bottom2'] + + self.quadRule = QuadratureRule.create_quadrature_rule_1D(4) + + maxContactNeighbors = 3 + self.interactionList = Contact.get_potential_interaction_list(self.surfM, self.surfI, self.mesh, + self.U, maxContactNeighbors) + self.closestSides,self.sideWeights = \ + Contact.compute_closest_edges_and_field_weights(self.mesh, self.U, + self.quadRule, + self.interactionList, + self.surfI) + + + def test_friction_search_static(self): + + coordsQ = Contact.compute_q_coordinates(self.mesh, self.U, self.quadRule, self.surfI) + coordsFromWeights = Contact.compute_q_coordinates_from_field_weights(self.mesh, self.U, + self.closestSides, + self.sideWeights) + + self.assertArrayNear(coordsQ, coordsFromWeights, 14) + + self.lam = np.ones_like(coordsQ[:,:,0]) + frictionEnergy = Contact.compute_friction_potential(self.mesh, self.U, + self.lam, frictionParams, + self.quadRule, self.surfI, + self.closestSides, self.sideWeights) + + self.assertNear(frictionEnergy, 0.0, 14) + + + def test_friction_search_after_motion(self): + + coordsQ0 = Contact.compute_q_coordinates(self.mesh, self.U, self.quadRule, self.surfI) + + block1Nodes = sort_uniq(self.mesh.conns[self.mesh.blocks['block1']].ravel()) + block2Nodes = sort_uniq(self.mesh.conns[self.mesh.blocks['block2']].ravel()) + + self.U = self.U.at[block1Nodes,0].add(0.1) + self.U = self.U.at[block2Nodes,0].add(-0.25) + + coordsQ = Contact.compute_q_coordinates(self.mesh, self.U, self.quadRule, self.surfI) + coordsFromWeights = Contact.compute_q_coordinates_from_field_weights(self.mesh, self.U, + self.closestSides, + self.sideWeights) + + self.assertArrayNear(coordsQ0, coordsFromWeights-np.array([0.1, 0.0]), 14) + self.assertArrayNear(coordsQ0, coordsQ+np.array([0.25, 0.0]), 14) + + +if __name__ == '__main__': + unittest.main() + + diff --git "a/optimism/test/contact/test_Contact.py\357\200\272Zone.Identifier" "b/optimism/test/contact/test_Contact.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git a/optimism/test/contact/test_Cpp.py b/optimism/test/contact/test_Cpp.py new file mode 100644 index 00000000..84b25da1 --- /dev/null +++ b/optimism/test/contact/test_Cpp.py @@ -0,0 +1,158 @@ +import unittest +import matplotlib.pyplot as plt + +from optimism.JaxConfig import * + +from optimism.contact import EdgeCpp +from optimism import SmoothFunctions +from ..TestFixture import TestFixture + +doPlotting=False + +def compute_grid_field(N, xmin, xmax, evaluation_function): + x = np.linspace(xmin, xmax, N) + y = np.linspace(xmin, xmax, N) + + X,Y = np.meshgrid(x,y) + V = np.array([vmap(lambda xi: np.array([xi,yi]))(x) for yi in y]) + Z = vmap(lambda vs: vmap(lambda v: evaluation_function(v))(vs))(V) + + return X,Y,Z + + +def edges_from_points(p0, p1, p2): + edge1 = np.array([p0, p1]) + edge2 = np.array([p1, p2]) + return np.array([edge1, edge2]) + + +class TestEdgeIntersection(TestFixture): + + def setUp(self): + + #self.v1 = np.array([0., 0.]) + #self.v2 = np.array([ 1., 1.]) + #self.v3 = np.array([ 1., 0,]) + + self.v1 = np.array([-1., -1.]) + self.v2 = np.array([ 1., 1.]) + self.v3 = np.array([ 1., -1,]) + + self.edges = edges_from_points(self.v3, self.v2, self.v1) + self.edgesR = edges_from_points(self.v1, self.v2, self.v3) + + self.xmin = -1.5 + self.xmax = 1.5 + + + def test_cpp_dist_interior(self): + p = np.array([0.5, -0.5]) + dist = EdgeCpp.cpp_distance(self.edges[1], p) + self.assertNear(dist, -np.sqrt(0.5), 14) + + + def test_cpp_dist_exterior(self): + p = np.array([-0.5, 0.5]) + dist = EdgeCpp.cpp_distance(self.edges[1], p) + self.assertNear(dist, np.sqrt(0.5), 14) + + + def test_cpp_dist_corner1(self): + p = np.array([1.5, 1.5]) + dist = EdgeCpp.cpp_distance(self.edges[1], p) + self.assertNear(dist, np.sqrt(0.5), 14) + + + def test_cpp_dist_corner2(self): + p = np.array([2.0, 1.0]) + dist = EdgeCpp.cpp_distance(self.edges[1], p) + self.assertNear(dist, -1.0, 14) + + + def plot_grid(self, X, Y, Z, edges, outname): + + if doPlotting: + plt.clf() + fig1, ax2 = plt.subplots(constrained_layout=True) + c = ax2.contourf(X, Y, Z, levels=[-0.6+0.1*i for i in range(13)]) + + ax2.set_title('Distance contours') + ax2.set_xlabel('x') + ax2.set_ylabel('y') + + cbar = fig1.colorbar(c) + cbar.ax.set_ylabel('distance') + + edge0 = edges[0] + edge1 = edges[1] + + plt.plot(edge0[:,0], + edge0[:,1], + color='k', linestyle='-', linewidth=2) + + plt.plot(edge1[:,0], + edge1[:,1], + color='k', linestyle='-', linewidth=2) + + plt.axis([self.xmin, self.xmax, self.xmin, self.xmax]) + + ax2.set_aspect('equal', adjustable='box') + if outname is not None: + plt.savefig(outname) + + + def test_smooth_1(self): + N = 301 + tol = 1e-1 + + X,Y,Z = compute_grid_field(N, self.xmin, self.xmax, lambda p: EdgeCpp.smooth_distance(self.edges, p, tol)) + self.plot_grid(X, Y, Z, self.edges, 'smooth1.png') + + + def test_smooth_2(self): + N = 301 + tol = 1e-1 + + X,Y,Z = compute_grid_field(N, self.xmin, self.xmax, lambda p: EdgeCpp.smooth_distance(self.edgesR, p, tol)) + self.plot_grid(X, Y, Z, self.edgesR, 'smooth2.png') + + + def test_limits(self): + N = 201 + tol = 6e-1 + Ys = np.linspace(1.8, 2.2, 21) + + v1 = np.array([0., 0.]) + v2 = np.array([ 1., 1.]) + + for i,y in enumerate(Ys): + v3 = np.array([2.0, y]) + edges = edges_from_points(v3, v2, v1) + smooth_dist = lambda p: EdgeCpp.smooth_distance(edges, p, tol) + dist = lambda p: grad(EdgeCpp.smooth_distance,0)(edges, p, tol)[0,0,0] + #dist = smooth_dist + X,Y,Z = compute_grid_field(N, self.xmin, self.xmax, + dist) + self.plot_grid(X, Y, Z, edges, 'smooth_'+str(i).zfill(2)+'.png') + + + def test_plot_smooth_min(self): + N = 201 + x = np.linspace(0.,1.,N) + y = 0.5 + + func = grad(SmoothFunctions.min,2) + + for eps in np.linspace(1e-5, 1e-1, 3): + m = vmap(func, (0,None,None))(x,y,eps) + plt.plot(x,m) + + m = vmap(func, (0,None,None))(x,y,0.) + plt.plot(x,m, '--k') + + if doPlotting: + plt.savefig('smooth.png') + + +if __name__ == '__main__': + unittest.main() diff --git "a/optimism/test/contact/test_Cpp.py\357\200\272Zone.Identifier" "b/optimism/test/contact/test_Cpp.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git a/optimism/test/contact/test_EdgeIntersection.py b/optimism/test/contact/test_EdgeIntersection.py new file mode 100644 index 00000000..66d18a2a --- /dev/null +++ b/optimism/test/contact/test_EdgeIntersection.py @@ -0,0 +1,144 @@ +import matplotlib.pyplot as plt +import unittest + +from optimism.JaxConfig import * +from optimism.contact import EdgeIntersection +from ..TestFixture import TestFixture + + +class TestEdgeIntersection(TestFixture): + + def setUp(self): + self.v1 = np.array([-1.,-1.]) + self.v2 = np.array([1.,1.]) + self.p = np.array([1.,-1.]) + self.raySmoothing = 1e-5 + + + def get_edge(self): + return np.array([self.v1, self.v2]) + + + def get_ray(self, normal): + return np.array([self.p, normal]) + + + def compute_ray_trace(self, normal): + edge = self.get_edge() + ray = self.get_ray(normal) + return EdgeIntersection.compute_valid_ray_trace_distance_smoothed(edge, ray, self.raySmoothing) + + + def test_valid_intersection(self): + normal = np.array([-1.,1.]) + + rayLength,edgeParCoord = self.compute_ray_trace(normal) + + pointOnEdge = self.v1 + edgeParCoord*(self.v2-self.v1) + pointOnRay = self.p + rayLength*normal + + self.assertArrayEqual(pointOnEdge, pointOnRay) + self.assertArrayEqual(pointOnEdge, np.array([0.,0.])) + self.assertEqual(edgeParCoord, 0.5) + self.assertEqual(rayLength , 1.0) + + + def test_valid_intersection_on_edge(self): + normal = np.array([0,1.]) + + rayLength,edgeParCoord = self.compute_ray_trace(normal) + pointOnEdge = self.v1 + edgeParCoord*(self.v2-self.v1) + pointOnRay = self.p + rayLength*normal + + self.assertArrayEqual(pointOnEdge, pointOnRay) + self.assertArrayEqual(pointOnEdge, np.array([1.,1.])) + self.assertEqual(edgeParCoord, 1.0) + self.assertEqual(rayLength, 2.0) + + + def test_smooth_gradient_on_either_side_of_right_edge(self): + ray_trace = lambda normal: self.compute_ray_trace(normal)[0] + ray_trace_grad = grad(ray_trace) + + eps = 1e-15 + normalIn = np.array([-eps,1.]) + normalOut = np.array([+eps,1.]) + + self.assertNear(ray_trace(normalIn), ray_trace(normalOut), 13) + self.assertArrayNear(ray_trace_grad(normalIn), ray_trace_grad(normalOut), 2) + + + def test_smooth_gradient_on_either_side_of_left_edge(self): + ray_trace = lambda normal: self.compute_ray_trace(normal)[0] + ray_trace_grad = grad(ray_trace) + + eps = 1e-15 + normalIn = np.array([-1.,eps]) + normalOut = np.array([-1.,-eps]) + + self.assertNear(ray_trace(normalIn), ray_trace(normalOut), 13) + self.assertArrayNear(ray_trace_grad(normalIn), ray_trace_grad(normalOut), 2) + + + def test_limit_of_ray_smoothing(self): + ray_trace = lambda normal: self.compute_ray_trace(normal)[0] + + normalRight = np.array([0.99*self.raySmoothing, 1.]) + self.assertTrue(ray_trace(normalRight) > 1.0) + self.assertTrue(ray_trace(normalRight) != np.inf) + + normalMoreRight = np.array([0.99999*self.raySmoothing, 1.]) + self.assertTrue(ray_trace(normalMoreRight) > 1.0) + self.assertTrue(ray_trace(normalMoreRight) != np.inf) + + normalWayRight = np.array([1.01*self.raySmoothing, 1.]) + self.assertTrue(ray_trace(normalWayRight) > 1.0) + + normalLeft = np.array([-1., -0.99*self.raySmoothing]) + self.assertTrue(ray_trace(normalLeft) > 1.0) + self.assertTrue(ray_trace(normalLeft) != np.inf) + + normalMoreLeft = np.array([-1., -0.99999*self.raySmoothing]) + self.assertTrue(ray_trace(normalMoreLeft) > 1.0) + self.assertTrue(ray_trace(normalMoreLeft) != np.inf) + + normalWayLeft = np.array([-1., -1.01*self.raySmoothing]) + self.assertTrue(ray_trace(normalWayLeft) > 1.0) + + + def get_ray_length_arg_x(self, xComp): + sol = EdgeIntersection.compute_intersection(self.get_edge(), self.get_ray([xComp, 1.0])) + return sol[1] + + + def get_ray_length_arg_y(self, yComp): + sol = EdgeIntersection.compute_intersection(self.get_edge(), self.get_ray([-1.0, yComp])) + return sol[1] + + + @unittest.skip + def test_plot(self): + N = 150 + xcomp = np.linspace(-0.5, 0.5, N) + rays = [self.get_ray_length_arg_x(x) for x in xcomp] + plt.plot(xcomp, rays) + plt.show() + + + @unittest.skip + def test_plot2(self): + N = 150 + ycomp = np.linspace(-0.5, 0.5, N) + rays = [self.get_ray_length_arg_y(x) for x in ycomp] + plt.plot(ycomp, rays) + plt.show() + + rayDer = [jit(grad(self.get_ray_length_arg_y))(x) for x in ycomp] + plt.plot(ycomp, rayDer) + plt.show() + + +if __name__ == '__main__': + unittest.main() + + diff --git "a/optimism/test/contact/test_EdgeIntersection.py\357\200\272Zone.Identifier" "b/optimism/test/contact/test_EdgeIntersection.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git a/optimism/test/contact/test_LevelsetConstraint.py b/optimism/test/contact/test_LevelsetConstraint.py new file mode 100644 index 00000000..5d46e349 --- /dev/null +++ b/optimism/test/contact/test_LevelsetConstraint.py @@ -0,0 +1,159 @@ +from optimism.JaxConfig import * + +from optimism import AlSolver +from optimism.ConstrainedObjective import ConstrainedObjective +from optimism.ConstrainedObjective import ConstrainedQuasiObjective +from . import EquationSolver_Immersed_2 as EqSolver +from optimism import Mesh +from optimism.material import LinearElastic +from optimism.contact import Friction +from optimism.contact import Levelset +from optimism.contact import LevelsetConstraint +from ..MeshFixture import MeshFixture, defaultProps +from optimism import Mechanics +from optimism import FunctionSpace +from optimism import QuadratureRule +import unittest + +props = {'elastic modulus': 1.0, + 'poisson ratio': 0.25} + +materialModel = LinearElastic.create_material_model_functions(props) + +materialModels = {'block': materialModel} + +settings = EqSolver.get_settings() +alSettings = AlSolver.get_settings() +frictionParams = Friction.Params(mu = 0.3, sReg = 1e-4) + +class TestLevelsetContactConstraint(MeshFixture): + + def setUp(self): + + self.Nx = 5 + self.Ny = 5 + self.xRange = [0.,1.] + self.yRange = [0.,1.] + + self.mesh, self.U = self.create_mesh_and_disp(self.Nx, self.Ny, + self.xRange, self.yRange, + lambda x : 1e-14*x) + + quadratureRule2d = QuadratureRule.create_quadrature_rule_on_triangle(degree=1) + fs = FunctionSpace.construct_function_space(self.mesh, quadratureRule2d) + + EBCs = [FunctionSpace.EssentialBC(nodeSet='left', component=0), + FunctionSpace.EssentialBC(nodeSet='left', component=1), + FunctionSpace.EssentialBC(nodeSet='right', component=0), + FunctionSpace.EssentialBC(nodeSet='right', component=1)] + + self.dofManager = FunctionSpace.DofManager(fs, dim=2, EssentialBCs=EBCs) + + self.quadRule = QuadratureRule.create_quadrature_rule_1D(2) + self.edges = self.mesh.sideSets['top'] + + self.Uu = self.dofManager.get_unknown_values(self.U) + self.Ubc = self.dofManager.get_bc_values(self.U) + + self.mechFuncs = Mechanics.create_multi_block_mechanics_functions(fs, + 'plane strain', + materialModels) + + stateVars = self.mechFuncs.compute_initial_state() + + self.energy_func = lambda Uu,p : \ + self.mechFuncs.compute_strain_energy(self.dofManager.create_field(Uu, self.Ubc), stateVars) + + self.constraint_func = lambda levelset, Uu: LevelsetConstraint. \ + compute_levelset_constraints(levelset, + self.dofManager.create_field(Uu, self.Ubc), + mesh=self.mesh, + quadRule=self.quadRule, + edges=self.edges).ravel() + + + def test_compute_all_positive_constraints_for_far_away_levelset(self): + levelset = partial(Levelset.sphere, xLoc=0.5, yLoc=2.0, R=0.2) + constraints=self.constraint_func(levelset, self.Uu) + self.assertTrue( np.all( constraints > 0.0 ) ) + + def test_some_positive_some_negative_constraints_for_small_sphere_on_edge(self): + levelset = partial(Levelset.sphere, xLoc=0.5, yLoc=1.0, R=0.2) + constraints=self.constraint_func(levelset, self.Uu) + + self.assertFalse( np.all( constraints > 0.0 ) ) + self.assertFalse( np.all( constraints < 0.0 ) ) + + def test_solve(self): + levelset = partial(Levelset.sphere, xLoc=0.5, yLoc=1.05, R=0.2) + + constraint_func = lambda x,p: self.constraint_func(levelset, x) + tol = 1e-4 + + #MeshPlot.plot_mesh_with_field(self.mesh, self.U, fast=False, direction=1, plotName='test.png') + + p = None + lam0 = 0.0*constraint_func(self.Uu, p) + kappa0 = 5 * np.ones(lam0.shape) + objective = ConstrainedObjective(self.energy_func, + constraint_func, + self.Uu, + p, + lam0, + kappa0) + + self.Uu = AlSolver.augmented_lagrange_solve(objective, self.Uu, p, alSettings, settings, useWarmStart=False) + self.U = self.dofManager.create_field(self.Uu, self.Ubc) + + #MeshPlot.plot_mesh_with_field(self.mesh, self.U, fast=False, direction=1, plotName='test.png') + + def test_friction(self): + + initialSurfacePoints = LevelsetConstraint. \ + compute_contact_point_coordinates(self.dofManager.create_field(self.Uu, self.Ubc), + self.mesh, + self.quadRule, + self.edges) + + initialLevelsetCenter = np.array([0.5, 1.05]) + + levelset = partial(Levelset.sphere, + xLoc=initialLevelsetCenter[0], + yLoc=initialLevelsetCenter[1], + R=0.2) + + levelsetMotion = np.array([0.2, 0.0]) + + energy_func = lambda Uu,lam,p: self.energy_func(Uu,p) + \ + LevelsetConstraint.compute_friction_potential(self.dofManager.create_field(Uu, self.Ubc), + initialSurfacePoints, + levelsetMotion, + lam, + self.mesh, + self.quadRule, + self.edges, + frictionParams) + + + constraint_func = lambda x,p: self.constraint_func(levelset, x) + tol = 1e-4 + + #MeshPlot.plot_mesh_with_field(self.mesh, self.U, fast=False, direction=1, plotName='test.png') + p = None + lam0 = 0.0*constraint_func(self.Uu, p) + kappa0 = 5 * np.ones(lam0.shape) + + objective = ConstrainedQuasiObjective(energy_func, constraint_func, + self.Uu, + p, + lam0, + kappa0) + + self.Uu = AlSolver.augmented_lagrange_solve(objective, self.Uu, p, alSettings, settings, useWarmStart=False) + self.U = self.dofManager.create_field(self.Uu, self.Ubc) + + +if __name__ == '__main__': + unittest.main() + + diff --git "a/optimism/test/contact/test_LevelsetConstraint.py\357\200\272Zone.Identifier" "b/optimism/test/contact/test_LevelsetConstraint.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git a/optimism/test/contact/test_MortarGeom.py b/optimism/test/contact/test_MortarGeom.py new file mode 100644 index 00000000..704bf3f1 --- /dev/null +++ b/optimism/test/contact/test_MortarGeom.py @@ -0,0 +1,124 @@ +from ..TestFixture import * +from optimism import QuadratureRule +from optimism.contact import MortarContact +import jax +import jax.numpy as np +import numpy as onp + +usingTestingFilter = False + + +def spline_ramp(eta): + # return 3*eta*eta - 2*eta*eta*eta # a less smoothing ramp + #K = np.array([[3,4,5],[6,12,20],[1,1,1]]) + #r = np.array([0,0,1]) + #abc = np.linalg.solve(K,r) + a = 10 + b = -15 + c = 6 + eta2 = eta*eta + eta3 = eta2*eta + return np.where(eta >= 1.0, 1.0, a * eta3 + b * eta2*eta2 + c * eta2*eta3) + + +def compute_error(edgeA, edgeB, xiA, xiB, g, f_common_normal): + normal = f_common_normal(edgeA, edgeB) + xA = MortarContact.eval_linear_field_on_edge(edgeA, xiA) + xB = MortarContact.eval_linear_field_on_edge(edgeB, xiB) + return np.linalg.norm(xA - xB + g * normal) + + +class TestMortarGeom(TestFixture): + + def setUp(self): + self.f_common_normal = MortarContact.compute_average_normal #compute_normal_from_a + + @unittest.skipIf(usingTestingFilter, '') + def testOffEdges(self): + edgeA = np.array([[0.1, 0.2], [-0.2, 0.3]]) + edgeB = np.array([[-0.2, 0.28], [0.15, -0.25]]) + xiA,xiB,g = MortarContact.compute_intersection(edgeA, edgeB, self.f_common_normal) + for i in range(len(g)): + err = compute_error(edgeA, edgeB, xiA[i], xiB[i], g[i], self.f_common_normal) + self.assertTrue(err < 1e-15) + + + @unittest.skipIf(usingTestingFilter, '') + def testEdgesWithCommonPoint(self): + edgeA = np.array([[0.1, 0.2], [-0.2, 0.3]]) + edgeB = np.array([[-0.2, 0.3], [0.15, -0.25]]) + xiA,xiB,g = MortarContact.compute_intersection(edgeA, edgeB, self.f_common_normal) + for i in range(len(g)): + err = compute_error(edgeA, edgeB, xiA[i], xiB[i], g[i], self.f_common_normal) + self.assertTrue(err < 1e-15) + + + @unittest.skipIf(usingTestingFilter, '') + def testEdgesWithTwoCommonPoints(self): + edgeA = np.array([[0.1, 0.2], [-0.2, 0.3]]) + edgeB = np.array([[-0.2, 0.3], [0.1, 0.2]]) + xiA,xiB,g = MortarContact.compute_intersection(edgeA, edgeB, self.f_common_normal) + for i in range(len(g)): + err = compute_error(edgeA, edgeB, xiA[i], xiB[i], g[i], self.f_common_normal) + self.assertTrue(err < 1e-15) + + + @unittest.skipIf(usingTestingFilter, '') + def testAreaIntegrals(self): + edgeA = np.array([[0.1, 0.1], [0.2, 0.1]]) + edgeB = np.array([[0.22, 0.0], [0.18, 0.0]]) + # eventually... + # can give options on integrating on common vs different surfaces + # can give options on quadrature rule on common surface + commonArea = MortarContact.integrate_with_mortar(edgeA, edgeB, self.f_common_normal, lambda xiA, xiB, g: g) + self.assertNear(commonArea, 0.002, 9) + + + @unittest.skipIf(True, 'Used to visually explore smoothing') + def testSpline(self): + from matplotlib import pyplot as plt + x = np.linspace(0.0,1.0,100) + y = jax.vmap(spline_ramp)(x) + #y = jax.vmap(smooth_linear,(0,None))(x,0.1) + plt.clf() + plt.plot(x,y) + plt.savefig('tmp.png') + + + @unittest.skipIf(False, 'Used to visually explore contact behavior') + def testMortarIntegralOneSided(self): + from matplotlib import pyplot as plt + + def integrate_multipliers(edgeA, edgeB, lambdaA, lambdaB): + def q(xiA, xiB, g): + lamA = MortarContact.eval_linear_field_on_edge(lambdaA, xiA) + lamB = MortarContact.eval_linear_field_on_edge(lambdaB, xiB) + return g * (lamA + lamB) + return MortarContact.integrate_with_mortar(edgeA, edgeB, self.f_common_normal, q, relativeSmoothingSize = 0.03) + + edgeA = np.array([[1.0, 0.1], [2.0, 0.1]]) + edgeB = np.array([[3.0, 0.0], [2.0, 0.0]]) + + lambdaA = np.array([1.0, 1.0]) + lambdaB = np.zeros(2) + + def gap_of_y(y): + ea = edgeA.at[0,0].set(y) + ea = ea.at[1,0].set(y+1.0) + return integrate_multipliers(ea, edgeB, lambdaA, lambdaB) + + N = 500 + edgeAy = np.linspace(0.9, 3.1, N) + energy = jax.vmap(gap_of_y)(edgeAy) + force = jax.vmap(jax.grad(gap_of_y))(edgeAy) + + plt.clf() + plt.plot(edgeAy, energy) + plt.savefig('energy.png') + plt.clf() + plt.plot(edgeAy, force) + plt.savefig('force.png') + + +if __name__ == '__main__': + unittest.main() \ No newline at end of file diff --git "a/optimism/test/contact/test_MortarGeom.py\357\200\272Zone.Identifier" "b/optimism/test/contact/test_MortarGeom.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git a/optimism/test/contact/test_NewtonGlobalization.py b/optimism/test/contact/test_NewtonGlobalization.py new file mode 100644 index 00000000..8947c31e --- /dev/null +++ b/optimism/test/contact/test_NewtonGlobalization.py @@ -0,0 +1,187 @@ +from numpy.random import rand + +from optimism.JaxConfig import * +from ..TestFixture import * +from optimism.NewtonSolver import * +from optimism import AlSolver +from optimism.Objective import Objective +from . import EquationSolver_Immersed_2 as EqSolver +from optimism.ConstrainedObjective import ConstrainedObjective + + +class TestQuadraticSolver(TestFixture): + + def setUp(self): + self.bounds = np.array([0.05, 0.75]) + self.samplePoints = np.linspace( self.bounds[0], self.bounds[1], 200 ) + + + def check_quadratic(self, ps, theta): + p = construct_quadratic(ps) + self.assertNear(p(0), ps[0], 14) + self.assertNear(p(1), ps[1], 14) + + pmin = min(p(self.samplePoints)) + self.assertTrue( p(theta) <= pmin ) + + + def test_constant(self): + ps = np.array([0.5, 0.5, 0.]) + theta = compute_min_p(ps, self.bounds) + self.assertNear(theta, self.bounds[1], 13) + self.check_quadratic(ps, theta) + + + def test_linear(self): + ps = np.array([0.5, 0.9, 0.4]) + theta = compute_min_p(ps, self.bounds) + self.assertNear(theta, self.bounds[0], 13) + self.check_quadratic(ps, theta) + + + def test_negative_linear(self): + ps = np.array([0.5, 0.1, -0.4]) + theta = compute_min_p(ps, self.bounds) + self.assertNear(theta, self.bounds[1], 13) + self.check_quadratic(ps, theta) + + + def test_negative_curvature(self): + ps = np.array([0.5, 0.1, 0.0]) + theta = compute_min_p(ps, self.bounds) + self.assertNear(theta, self.bounds[1], 13) + self.check_quadratic(ps, theta) + + + def test_positive_curvature(self): + ps = np.array([0.25, 0.25, -1.]) + theta = compute_min_p(ps, self.bounds) + self.assertNear(theta, 0.5, 13) + self.check_quadratic(ps, theta) + + + def test_positive_curvature2(self): + ps = np.array([0.5, 2.5, 1.]) + theta = compute_min_p(ps, self.bounds) + self.assertNear(theta, self.bounds[0], 13) + self.check_quadratic(ps, theta) + + + def test_positive_curvature3(self): + ps = np.array([2.5, 0.5, -2.5]) + theta = compute_min_p(ps, self.bounds) + self.assertNear(theta, self.bounds[1], 13) + self.check_quadratic(ps, theta) + + +def objective(x): + return x[0]*(x[0]+1.0) + 0.3*x[1]*(x[1]-0.2) + 0.2*x[2]*(x[2]-0.5) + x[3]*x[3] + 0.00*(x[0]*x[0]*x[1]*x[1] + x[0]*x[1]) + + +def constraint(x): + return np.array([x@x-1.0, x[0]+x[1]-0.5, 1.0-x[0]]) + + +sizex=4 +sizelam=3 + + +def fischer_burmeister(c, l): + return np.sqrt(c**2+l**2) - c - l + + +dObjective = grad(objective) +dConstraint = grad(lambda x,lam: lam @ constraint(x)) + + +def residual(xlam): + x = xlam[:sizex] + lam = xlam[sizex:] + c = constraint(x) + return np.hstack( (dObjective(x) - dConstraint(x, lam), fischer_burmeister(c, lam)) ) + + +def create_linear_op(residual): + return lambda x,v: jvp( residual , (x,), (v,) )[1] + + +@jit +def linear_op(xlam, v): + dr_func = create_linear_op(residual) + return dr_func(xlam, v) + + +def my_func(x): + return 3.*x**2 + 2.*x + 3. + + +class TestGMRESSolver(TestFixture): + + def setUp(self): + self.x = np.ones(sizex) + self.lam = np.ones(sizelam) + + self.etak=1e-4 + self.t=1e-4 + + tol = 5e-9 + self.subProblemSettings = EqSolver.get_settings() + self.alSettings = AlSolver.get_settings() + + + def test_newton_step(self): + xl = np.hstack( (self.x, self.lam) ) + r0 = np.linalg.norm( residual(xl) ) + + xl += newton_step(residual, lambda v: linear_op(xl,v), xl)[0] + r1 = np.linalg.norm( residual(xl) ) + + self.assertTrue( r1 < 10*r0 ) + + + def test_globalized_newton_step_with_cubic(self): + residual = vmap(lambda x: np.sqrt(my_func(x))) + linear_op = create_linear_op(residual) + x = np.array([0.1]) + r0 = np.linalg.norm( residual(x) ) + x += globalized_newton_step(residual, lambda v: linear_op(x,v), x, self.etak, self.t) + r1 = np.linalg.norm( residual(x) ) + + self.assertTrue( r1 < r0 ) + + + def test_globalized_newton_step_nonconvex(self): + xl = np.hstack( (self.x, self.lam) ) + r0 = np.linalg.norm( residual(xl) ) + + xl += globalized_newton_step(residual, lambda v: linear_op(xl,v), xl, self.etak, self.t) + r1 = np.linalg.norm( residual(xl) ) + + self.assertTrue( r1 < r0 ) + + for count in range(2): + xl += globalized_newton_step(residual, lambda v: linear_op(xl,v), xl, self.etak, self.t) + r1 = np.linalg.norm( residual(xl) ) + + + def test_al_solver(self): + # Guess a bit randomly at a good initial penalty stiffness + p = None + unconstrainedObjective = Objective(lambda x,p: objective(x), self.x, p) + randRhs = np.array(rand(self.x.size)) + randRhs *= 0.25 / np.linalg.norm(randRhs) + + penalty0 = unconstrainedObjective.hessian_vec(self.x, randRhs) @ randRhs + + alObjective = ConstrainedObjective(lambda x,p: objective(x), + lambda x,p: constraint(x), + self.x, + p, + self.lam, + penalty0 * np.ones(self.lam.shape)) + + AlSolver.augmented_lagrange_solve(alObjective, self.x, p, + self.alSettings, self.subProblemSettings, useWarmStart=False) + +if __name__ == '__main__': + unittest.main() diff --git "a/optimism/test/contact/test_NewtonGlobalization.py\357\200\272Zone.Identifier" "b/optimism/test/contact/test_NewtonGlobalization.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git a/optimism/test/contact/test_Search.py b/optimism/test/contact/test_Search.py new file mode 100644 index 00000000..5f63893a --- /dev/null +++ b/optimism/test/contact/test_Search.py @@ -0,0 +1,98 @@ +import unittest +from optimism.JaxConfig import * +from optimism import Mesh +from optimism import Surface +from optimism import QuadratureRule +from optimism.contact import Search +from optimism.contact import EdgeIntersection +from ..MeshFixture import MeshFixture + + +def get_best_overlap_vector(mesh, edge, listOfEdges): + + edgeCoords = Surface.get_coords(mesh, edge) + normal = Surface.compute_normal(edgeCoords) + edgeCenter = np.average(edgeCoords, axis=0) + + integrationRay = np.array([edgeCenter, normal]) + + minDistance = np.inf + + for neighbor in listOfEdges: + if not np.all(edge==neighbor): + neighborCoords = Surface.get_coords(mesh, neighbor) + distance, parCoord = EdgeIntersection.compute_valid_ray_trace_distance(neighborCoords, integrationRay) + minDistance = distance if distance < minDistance else minDistance + + return minDistance * integrationRay[1] + + +class TestDoubleMeshFixture(MeshFixture): + + def setUp(self): + self.tol = 1e-7 + self.targetDispGrad = np.array([[0.1, -0.2],[0.4, -0.1]]) + + m1 = self.create_mesh_and_disp(3, 5, [0., 1.], [0., 1.], + lambda x : self.targetDispGrad.dot(x)) + + m2 = self.create_mesh_and_disp(2, 4, [1.1, 2.], [0.,1.], + lambda x : self.targetDispGrad.dot(x)) + + self.mesh, self.disp = Mesh.combine_mesh(m1,m2) + + + def is_integration_edge(self, x): + tol = self.tol + return np.all( x[:,0] > 1.0 - tol ) and np.all( x[:,0] < 1.0 + tol ) + + + def is_contact_edge(self, x): + tol = self.tol + return np.all( x[:,0] > 1.0 - tol ) and np.all( x[:,0] < 1.1 + tol ) + + + def test_correct_number_of_edges_created_for_contact(self): + edges = Surface.create_edges(self.mesh.coords, self.mesh.conns, self.is_contact_edge) + self.assertEqual(7, edges.shape[0]) + + + def test_surface_integral_of_linears(self): + edges = Surface.create_edges(self.mesh.coords, self.mesh.conns, self.is_integration_edge) + quadratureRule = QuadratureRule.create_quadrature_rule_1D(1) + I = 0.0 + for edge in edges: + coords = Surface.get_coords(self.mesh, edge) + I += Surface.integrate_function(quadratureRule, coords, lambda x: x[1]) + self.assertNear(I, 0.5, 14) + + + def test_surface_integral_of_quadratics(self): + edges = Surface.create_edges(self.mesh.coords, self.mesh.conns, self.is_integration_edge) + quadratureRule = QuadratureRule.create_quadrature_rule_1D(2) + I = 0.0 + for edge in edges: + coords = Surface.get_coords(self.mesh, edge) + I += Surface.integrate_function(quadratureRule, coords, lambda x: 3.*x[1]**2) + self.assertNear(I, 1.0, 14) + + + def test_contact_distance_constraint_evaluation(self): + mesh = self.mesh + edges = Surface.create_edges(self.mesh.coords, self.mesh.conns, self.is_contact_edge) + coords = mesh.coords + + for edge in edges: + minOverlapVector = get_best_overlap_vector(mesh, edge, edges) + + edgeCoords = Surface.get_coords(mesh, edge) + if edgeCoords[0][0] < 1 + self.tol: + self.assertArrayNear(minOverlapVector, np.array([0.1, 0.0]), 12) + elif edgeCoords[0][0] > 1.1 - self.tol: + self.assertArrayNear(minOverlapVector, np.array([-0.1, 0.0]), 12) + else: + print("There should be no other edges\n") + + +if __name__ == '__main__': + unittest.main() diff --git "a/optimism/test/contact/test_Search.py\357\200\272Zone.Identifier" "b/optimism/test/contact/test_Search.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git a/optimism/test/contact/test_TwoBodyContact.py b/optimism/test/contact/test_TwoBodyContact.py new file mode 100644 index 00000000..216a8b36 --- /dev/null +++ b/optimism/test/contact/test_TwoBodyContact.py @@ -0,0 +1,119 @@ +import unittest +from optimism.JaxConfig import * +from optimism import Mesh +from optimism import QuadratureRule +from optimism import Surface +from optimism.contact import EdgeIntersection +from optimism import FunctionSpace +from ..MeshFixture import MeshFixture +from optimism.contact import Contact +from matplotlib import pyplot as plt + + +def get_side_set_segments(mesh, sideset, disp): + + def get_current_edge_segments(edge): + XOnEdge = Mesh.get_edge_coords(mesh, edge) + uOnEdge = Mesh.get_edge_field(mesh, edge, disp) + xOnEdge = XOnEdge + uOnEdge + return vmap(lambda x, y: np.array([x,y]))(xOnEdge[:-1], xOnEdge[1:]) + + segments = vmap(get_current_edge_segments)(sideset) + return segments.reshape((segments.shape[0]*segments.shape[1], segments.shape[2], segments.shape[3])) + + +def get_best_overlap_vector(mesh, edge, listOfEdges): + edgeCoords = Surface.get_coords(mesh, edge) + normal = Surface.compute_normal(edgeCoords) + edgeCenter = np.average(edgeCoords, axis=0) + + integrationRay = np.array([edgeCenter, normal]) + + minDistance = np.inf + for neighbor in listOfEdges: + if not np.all(edge==neighbor): + neighborCoords = Surface.get_coords(mesh, neighbor) + distance, parCoord = EdgeIntersection.compute_valid_ray_trace_distance(neighborCoords, integrationRay) + minDistance = distance if distance < minDistance else minDistance + + return minDistance * integrationRay[1] + + +class TwoBodyContactFixture(MeshFixture): + + def setUp(self): + self.tol = 1e-7 + self.targetDispGrad = np.array([[0.1, -0.2],[0.4, -0.1]]) + + m1 = self.create_mesh_and_disp(3, 5, [0.0, 1.0], [0.0, 1.0], + lambda x : self.targetDispGrad.dot(x), '1') + + m2 = self.create_mesh_and_disp(2, 4, [0.9, 2.0], [0.0, 1.0], + lambda x : self.targetDispGrad.dot(x), '2') + + self.mesh, _ = Mesh.combine_mesh(m1, m2) + order=1 # this quadrature on segment contact only works with 1st order meshes + self.mesh = Mesh.create_higher_order_mesh_from_simplex_mesh(self.mesh, order=order, copyNodeSets=False) + self.disp = np.zeros(self.mesh.coords.shape) + + #edgesTop1 = get_side_set_segments(self.mesh, self.mesh.sideSets['right1'], self.disp) + #edgesTop2 = get_side_set_segments(self.mesh, self.mesh.sideSets['left2'], self.disp) + + #plt.clf() + #plt.plot( edgesTop1[:,:,0].T, edgesTop1[:,:,1].T, linestyle='-', color='y', markerfacecolor='red', marker='o') + #plt.plot( edgesTop2[:,:,0].T, edgesTop2[:,:,1].T, linestyle='-', color='g', markerfacecolor='blue', marker='o') + #plt.show() + + nodeSets = Mesh.create_nodesets_from_sidesets(self.mesh) + self.mesh = Mesh.mesh_with_nodesets(self.mesh, nodeSets) + self.quadRule = QuadratureRule.create_quadrature_rule_1D(3) + + + @unittest.skipIf(False, '') + def test_combining_nodesets(self): + self.assertArrayEqual( self.mesh.nodeSets['top1'], [12,13,14] ) + numNodesMesh1 = 15 + self.assertArrayEqual( self.mesh.nodeSets['top2'], np.array([6,7])+numNodesMesh1 ) + + @unittest.skipIf(False, '') + def test_combining_sidesets(self): + self.assertArrayEqual( self.mesh.sideSets['top1'], np.array([[7,1],[15,1]]) ) + numElemsMesh1 = 2*8 + self.assertArrayEqual( self.mesh.sideSets['top2'], np.array([[5+numElemsMesh1,1]]) ) + + + def plot_solution(self, plotName): + import VTKWriter + + writer = VTKWriter.VTKWriter(self.mesh, baseFileName=plotName) + writer.add_nodal_field(name='displacement', + nodalData=self.disp, + fieldType=VTKWriter.VTKFieldType.VECTORS) + + bcs = np.array(self.dofManager.isBc, dtype=int) + writer.add_nodal_field(name='bcs', + nodalData=bcs, + fieldType=VTKWriter.VTKFieldType.VECTORS, + dataType=VTKWriter.VTKDataType.INT) + writer.write() + + + @unittest.skipIf(False, '') + def test_contact_search(self): + self.disp = 0.0*self.disp + + sideA = self.mesh.sideSets['right1'] + sideB = self.mesh.sideSets['left2'] + + interactionList = Contact.get_potential_interaction_list(sideA, sideB, self.mesh, self.disp, 3) + + minDists = Contact.compute_closest_distance_to_each_side(self.mesh, self.disp, self.quadRule, + interactionList, + sideB) + + self.assertArrayNear(minDists, -0.1*np.ones((3,2)), 9) + + +if __name__ == '__main__': + unittest.main() + diff --git "a/optimism/test/contact/test_TwoBodyContact.py\357\200\272Zone.Identifier" "b/optimism/test/contact/test_TwoBodyContact.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git a/optimism/test/contact/test_TwoBodyMortarContact.py b/optimism/test/contact/test_TwoBodyMortarContact.py new file mode 100644 index 00000000..e4ca226b --- /dev/null +++ b/optimism/test/contact/test_TwoBodyMortarContact.py @@ -0,0 +1,81 @@ +from optimism.JaxConfig import * +from optimism import Mesh +from ..MeshFixture import MeshFixture +from optimism.contact import MortarContact +import unittest + +class TwoBodyContactFixture(MeshFixture): + + def setUp(self): + self.targetDispGrad = np.array([[0.1, -0.2],[0.4, -0.1]]) + + m1 = self.create_mesh_and_disp(3, 5, [0.0, 1.0], [0.0, 1.0], + lambda x : self.targetDispGrad.dot(x), '1') + + m2 = self.create_mesh_and_disp(2, 8, [0.9, 2.0], [0.0, 1.0], + lambda x : self.targetDispGrad.dot(x), '2') + + self.mesh, _ = Mesh.combine_mesh(m1, m2) + order=2 + self.mesh = Mesh.create_higher_order_mesh_from_simplex_mesh(self.mesh, order=order, copyNodeSets=False) + self.disp = np.zeros(self.mesh.coords.shape) + + sideA = self.mesh.sideSets['right1'] + sideB = self.mesh.sideSets['left2'] + + self.segmentConnsA = MortarContact.get_facet_connectivities(self.mesh, sideA) + self.segmentConnsB = MortarContact.get_facet_connectivities(self.mesh, sideB) + + + def plot_solution(self, plotName): + from optimism import VTKWriter + writer = VTKWriter.VTKWriter(self.mesh, baseFileName=plotName) + writer.add_nodal_field(name='displacement', + nodalData=self.disp, + fieldType=VTKWriter.VTKFieldType.VECTORS) + writer.add_contact_edges(self.segmentConnsA) + writer.add_contact_edges(self.segmentConnsB) + writer.write() + + + @unittest.skipIf(False, '') + def test_contact_search(self): + #self.plot_solution('mesh') + + neighborList = MortarContact.get_closest_neighbors(self.segmentConnsA, self.segmentConnsB, self.mesh, self.disp, 5) + + def compute_overlap_area(segB, neighborSegsA): + def compute_area_for_segment_pair(segB, indexA): + segA = self.segmentConnsA[indexA] + coordsSegB = self.mesh.coords[segB] + self.disp[segB] + coordsSegA = self.mesh.coords[segA] + self.disp[segA] + return MortarContact.integrate_with_mortar(coordsSegB, coordsSegA, MortarContact.compute_average_normal, lambda xiA, xiB, gap: 1.0, 1e-9) + + return np.sum(vmap(compute_area_for_segment_pair, (None,0))(segB, neighborSegsA)) + + totalSum = np.sum(vmap(compute_overlap_area)(self.segmentConnsB, neighborList)) + self.assertNear(totalSum, 1.0, 5) + + + def test_contact_constraints(self): + #self.plot_solution('mesh') + + neighborList = MortarContact.get_closest_neighbors(self.segmentConnsA, self.segmentConnsB, self.mesh, self.disp, 5) + + nodalGapField = MortarContact.assemble_area_weighted_gaps(self.mesh.coords, self.disp, self.segmentConnsA, self.segmentConnsB, neighborList, MortarContact.compute_average_normal) + nodalAreaField = MortarContact.assemble_nodal_areas(self.mesh.coords, self.disp, self.segmentConnsA, self.segmentConnsB, neighborList, MortarContact.compute_average_normal) + + nodalGapField = vmap(lambda x, d : np.where( d > 0, x / d, 0.0))(nodalGapField, nodalAreaField) + + mask = np.ones(len(nodalAreaField), dtype=np.int8) + nodesB = np.unique(np.concatenate(self.segmentConnsB)) + mask = mask.at[nodesB].set(0) + + self.assertNear(np.sum(nodalAreaField), 1.0, 8) + self.assertNear(np.sum(nodalAreaField), np.sum(nodalAreaField[nodesB]), 14) + self.assertNear(np.sum(nodalAreaField[mask]), 0.0, 14) + + +if __name__ == '__main__': + unittest.main() + diff --git "a/optimism/test/contact/test_TwoBodyMortarContact.py\357\200\272Zone.Identifier" "b/optimism/test/contact/test_TwoBodyMortarContact.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git "a/optimism/test/ill_conditioned_dot_product_data.npz\357\200\272Zone.Identifier" "b/optimism/test/ill_conditioned_dot_product_data.npz\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git "a/optimism/test/ill_conditioned_sum_data.npz\357\200\272Zone.Identifier" "b/optimism/test/ill_conditioned_sum_data.npz\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git a/optimism/test/inverse/FiniteDifferenceFixture.py b/optimism/test/inverse/FiniteDifferenceFixture.py new file mode 100644 index 00000000..04a2d2a8 --- /dev/null +++ b/optimism/test/inverse/FiniteDifferenceFixture.py @@ -0,0 +1,67 @@ +from ..MeshFixture import MeshFixture +from collections import namedtuple +import numpy as onp + +class FiniteDifferenceFixture(MeshFixture): + def assertFiniteDifferenceCheckHasVShape(self, errors, tolerance=1e-6): + minError = min(errors) + self.assertLess(minError, tolerance, "Smallest finite difference error not less than tolerance.") + self.assertLess(minError, errors[0], "Finite difference error does not decrease from initial step size.") + self.assertLess(minError, errors[-1], "Finite difference error does not increase after reaching minimum. Try more finite difference steps.") + + def build_direction_vector(self, numDesignVars, seed=123): + + onp.random.seed(seed) + directionVector = onp.random.uniform(-1.0, 1.0, numDesignVars) + normVector = directionVector / onp.linalg.norm(directionVector) + + return onp.array(normVector) + + def compute_finite_difference_error(self, stepSize, initialParameters): + storedState = self.forward_solve(initialParameters) + originalObjective = self.compute_objective_function(storedState, initialParameters) + gradient = self.compute_gradient(storedState, initialParameters) + + directionVector = self.build_direction_vector(initialParameters.shape[0]) + directionalDerivative = onp.tensordot(directionVector, gradient, axes=1) + + perturbedParameters = initialParameters + stepSize * directionVector + storedState = self.forward_solve(perturbedParameters) + perturbedObjective = self.compute_objective_function(storedState, perturbedParameters) + + fd_value = (perturbedObjective - originalObjective) / stepSize + error = abs(directionalDerivative - fd_value) + + return error + + def compute_finite_difference_errors(self, stepSize, steps, initialParameters, printOutput=True): + storedState = self.forward_solve(initialParameters) + originalObjective = self.compute_objective_function(storedState, initialParameters) + gradient = self.compute_gradient(storedState, initialParameters) + + directionVector = self.build_direction_vector(initialParameters.shape[0]) + directionalDerivative = onp.tensordot(directionVector, gradient, axes=1) + + fd_values = [] + errors = [] + for i in range(0, steps): + perturbedParameters = initialParameters + stepSize * directionVector + storedState = self.forward_solve(perturbedParameters) + perturbedObjective = self.compute_objective_function(storedState, perturbedParameters) + + fd_value = (perturbedObjective - originalObjective) / stepSize + fd_values.append(fd_value) + + error = abs(directionalDerivative - fd_value) + errors.append(error) + + stepSize *= 1e-1 + + if printOutput: + print("\n grad'*dir | FD approx | abs error") + print("--------------------------------------------------------------------------------") + for i in range(0, steps): + print(f" {directionalDerivative} | {fd_values[i]} | {errors[i]}") + + return errors + diff --git "a/optimism/test/inverse/FiniteDifferenceFixture.py\357\200\272Zone.Identifier" "b/optimism/test/inverse/FiniteDifferenceFixture.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git a/optimism/test/inverse/__init__.py b/optimism/test/inverse/__init__.py new file mode 100644 index 00000000..e69de29b diff --git "a/optimism/test/inverse/__init__.py\357\200\272Zone.Identifier" "b/optimism/test/inverse/__init__.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git a/optimism/test/inverse/test_Hyperelastic_gradient_checks.py b/optimism/test/inverse/test_Hyperelastic_gradient_checks.py new file mode 100644 index 00000000..7e0e7041 --- /dev/null +++ b/optimism/test/inverse/test_Hyperelastic_gradient_checks.py @@ -0,0 +1,212 @@ +import unittest + +import jax +import jax.numpy as np +import numpy as onp +from scipy.sparse import linalg + +from . import EquationSolver_Immersed_2 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 +from optimism.material import Neohookean + +from .FiniteDifferenceFixture import FiniteDifferenceFixture + +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): + dispGrad0 = np.array([[0.4, -0.2], + [-0.04, 0.68]]) + self.initialMesh, self.U = self.create_mesh_and_disp(4,4,[0.,1.],[0.,1.], + lambda x: dispGrad0@x) + + shearModulus = 0.855 # MPa + bulkModulus = 100*shearModulus # MPa + youngModulus = 9.0*bulkModulus*shearModulus / (3.0*bulkModulus + shearModulus) + poissonRatio = (3.0*bulkModulus - 2.0*shearModulus) / 2.0 / (3.0*bulkModulus + shearModulus) + props = { + 'elastic modulus': youngModulus, + 'poisson ratio': poissonRatio, + 'version': 'coupled' + } + self.materialModel = Neohookean.create_material_model_functions(props) + + self.quadRule = QuadratureRule.create_quadrature_rule_on_triangle(degree=1) + + self.EBCs = [FunctionSpace.EssentialBC(nodeSet='all_boundary', component=0), + FunctionSpace.EssentialBC(nodeSet='all_boundary', component=1)] + + self.steps = 2 + + def forward_solve(self, parameters): + + self.mesh = Mesh.construct_mesh_from_basic_data(parameters.reshape(self.initialMesh.coords.shape),\ + self.initialMesh.conns, self.initialMesh.blocks,\ + self.initialMesh.nodeSets, self.initialMesh.sideSets) + + functionSpace = FunctionSpace.construct_function_space(self.mesh, self.quadRule) + mechFuncs = Mechanics.create_mechanics_functions(functionSpace, + "plane strain", + self.materialModel) + self.dofManager = FunctionSpace.DofManager(functionSpace, 2, self.EBCs) + Ubc = self.dofManager.get_bc_values(self.U) + + Ubc_inc = Ubc / self.steps + ivs = mechFuncs.compute_initial_state() + p = Objective.Params(bc_data=np.zeros(Ubc.shape), state_data=ivs) + + def compute_energy(Uu, p): + U = self.dofManager.create_field(Uu, p.bc_data) + internalVariables = p.state_data + return mechFuncs.compute_strain_energy(U, internalVariables) + + Uu = 0.0*self.dofManager.get_unknown_values(self.U) + self.objective = Objective.Objective(compute_energy, Uu, p) + + storedState = [] + storedState.append((Uu, p)) + + for step in range(1, self.steps+1): + p = Objective.param_index_update(p, 0, step*Ubc_inc) + Uu, solverSuccess = EqSolver.nonlinear_equation_solve(self.objective, Uu, p, EqSolver.get_settings(), useWarmStart=False) + storedState.append((Uu, p)) + + return storedState + + 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) + 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 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) + + def total_work_increment(self, Uu, Uu_prev, p, p_prev, coords, nodal_forces): + index = (self.mesh.nodeSets['left'], 1) # arbitrarily choosing left side nodeset for reaction force + + U = self.dofManager.create_field(Uu, p.bc_data) + force = np.array(nodal_forces(U, p, coords).at[index].get()) + disp = U.at[index].get() + + U_prev = self.dofManager.create_field(Uu_prev, p_prev.bc_data) + force_prev = np.array(nodal_forces(U_prev, p_prev, coords).at[index].get()) + disp_prev = U_prev.at[index].get() + + return 0.5*np.tensordot((force + force_prev),(disp - disp_prev), axes=1) + + def total_work_objective(self, storedState, parameters): + parameters = parameters.reshape(self.mesh.coords.shape) + 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, 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): + 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)) + + parameters = parameters.reshape(self.mesh.coords.shape) + gradient = np.zeros(parameters.shape) + adjointLoad = np.zeros(storedState[0][0].shape) + + for step in reversed(range(1, self.steps+1)): + Uu = storedState[step][0] + p = storedState[step][1] + + 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, energyFuncs.nodal_forces) + + adjointLoad -= df_du + + n = self.dofManager.get_unknown_size() + self.objective.p = p # have to update parameters to get precond to work + self.objective.update_precond(Uu) # update preconditioner for use in cg (will converge in 1 iteration as long as the preconditioner is not approximate) + dRdu = linalg.LinearOperator((n, n), lambda V: onp.asarray(self.objective.hessian_vec(Uu, V))) + dRdu_decomp = linalg.LinearOperator((n, n), lambda V: onp.asarray(self.objective.apply_precond(V))) + adjointVector = linalg.cg(dRdu, onp.array(adjointLoad, copy=False), rtol=1e-10, atol=0.0, M=dRdu_decomp)[0] + + gradient += df_dx + gradient += residualInverseFuncs.residual_jac_coords_vjp(Uu, p, parameters, adjointVector) + + adjointLoad = -df_dun + # The action dRdU * lambda (the same as lambda^T * dRdU) - For Hyperelastic the residual is not dependent on Uu_prev, so don't actually need this term + + return gradient.ravel() + + + + def test_self_adjoint_gradient(self): + self.compute_objective_function = self.strain_energy_objective + self.compute_gradient = self.strain_energy_gradient + + stepSize = 1e-7 + + error = self.compute_finite_difference_error(stepSize, self.initialMesh.coords.ravel()) + self.assertLessEqual(error, 1e-7) + + @unittest.expectedFailure + def test_non_self_adjoint_gradient_without_adjoint_solve(self): + self.compute_objective_function = self.total_work_objective + self.compute_gradient = self.total_work_gradient_just_jax + + stepSize = 1e-7 + + error = self.compute_finite_difference_error(stepSize, self.initialMesh.coords.ravel()) + self.assertLessEqual(error, 1e-7) + + def test_non_self_adjoint_gradient_with_adjoint_solve(self): + self.compute_objective_function = self.total_work_objective + self.compute_gradient = self.total_work_gradient_with_adjoint + + stepSize = 1e-7 + + error = self.compute_finite_difference_error(stepSize, self.initialMesh.coords.ravel()) + self.assertLessEqual(error, 1e-7) + + + +if __name__ == '__main__': + unittest.main() \ No newline at end of file diff --git "a/optimism/test/inverse/test_Hyperelastic_gradient_checks.py\357\200\272Zone.Identifier" "b/optimism/test/inverse/test_Hyperelastic_gradient_checks.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git a/optimism/test/inverse/test_J2Plastic_gradient_checks.py b/optimism/test/inverse/test_J2Plastic_gradient_checks.py new file mode 100644 index 00000000..c6a6f5b3 --- /dev/null +++ b/optimism/test/inverse/test_J2Plastic_gradient_checks.py @@ -0,0 +1,294 @@ +import unittest + +import jax +import jax.numpy as np +import numpy as onp +from scipy.sparse import linalg + +from . import EquationSolver_Immersed_2 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 +from optimism.material import J2Plastic as J2 + +from .FiniteDifferenceFixture import FiniteDifferenceFixture + +from optimism.inverse import MechanicsInverse +from optimism.inverse import AdjointFunctionSpace +from collections import namedtuple + +EnergyFunctions = namedtuple('EnergyFunctions', + ['energy_function_coords', + 'nodal_forces']) + +class J2GlobalMeshAdjointSolveFixture(FiniteDifferenceFixture): + def setUp(self): + dispGrad0 = np.array([[0.4, -0.2], + [-0.04, 0.68]]) + self.initialMesh, self.U = self.create_mesh_and_disp(4,4,[0.,1.],[0.,1.], + lambda x: dispGrad0@x) + + E = 100.0 + poisson = 0.321 + H = 1e-2 * E + Y0 = 0.3 * E + props = { + 'elastic modulus': E, + 'poisson ratio': poisson, + 'yield strength': Y0, + 'kinematics': 'small deformations', + 'hardening model': 'linear', + 'hardening modulus': H + } + self.materialModel = J2.create_material_model_functions(props) + + self.quadRule = QuadratureRule.create_quadrature_rule_on_triangle(degree=1) + + self.EBCs = [FunctionSpace.EssentialBC(nodeSet='all_boundary', component=0), + FunctionSpace.EssentialBC(nodeSet='all_boundary', component=1)] + + self.steps = 2 + + def forward_solve(self, parameters): + self.mesh = Mesh.construct_mesh_from_basic_data(parameters.reshape(self.initialMesh.coords.shape),\ + self.initialMesh.conns, self.initialMesh.blocks,\ + self.initialMesh.nodeSets, self.initialMesh.sideSets) + + functionSpace = FunctionSpace.construct_function_space(self.mesh, self.quadRule) + mechFuncs = Mechanics.create_mechanics_functions(functionSpace, + "plane strain", + self.materialModel) + self.dofManager = FunctionSpace.DofManager(functionSpace, 2, self.EBCs) + Ubc = self.dofManager.get_bc_values(self.U) + + Ubc_inc = Ubc / self.steps + ivs = mechFuncs.compute_initial_state() + p = Objective.Params(bc_data=np.zeros(Ubc.shape), state_data=ivs) + + def compute_energy(Uu, p): + U = self.dofManager.create_field(Uu, p.bc_data) + internalVariables = p.state_data + return mechFuncs.compute_strain_energy(U, internalVariables) + + Uu = 0.0*self.dofManager.get_unknown_values(self.U) + self.objective = Objective.Objective(compute_energy, Uu, p) + + storedState = [] + storedState.append((Uu, p)) + + for step in range(1, self.steps+1): + p = Objective.param_index_update(p, 0, step*Ubc_inc) + Uu, solverSuccess = EqSolver.nonlinear_equation_solve(self.objective, Uu, p, EqSolver.get_settings(), useWarmStart=False) + U = self.dofManager.create_field(Uu, p.bc_data) + ivs = mechFuncs.compute_updated_internal_variables(U, p.state_data) + p = Objective.param_index_update(p, 1, ivs) + storedState.append((Uu, p)) + + return storedState + + def setup_energy_functions(self): + shapeOnRef = Interpolants.compute_shapes(self.mesh.parentElement, self.quadRule.xigauss) + + def energy_function_all_dofs(U, ivs, 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) + return mech_funcs.compute_strain_energy(U, ivs) + + def energy_function_coords(Uu, p, ivs, coords): + U = self.dofManager.create_field(Uu, p.bc_data) + return energy_function_all_dofs(U, ivs, coords) + + nodal_forces = jax.grad(energy_function_all_dofs, argnums=0) + + return EnergyFunctions(energy_function_coords, jax.jit(nodal_forces)) + + def total_work_increment(self, Uu, Uu_prev, ivs, ivs_prev, p, p_prev, coordinates, nodal_forces): + index = (self.mesh.nodeSets['left'], 1) # arbitrarily choosing left side nodeset for reaction force + + U = self.dofManager.create_field(Uu, p.bc_data) + force = np.array(nodal_forces(U, ivs, coordinates).at[index].get()) + disp = U.at[index].get() + + U_prev = self.dofManager.create_field(Uu_prev, p_prev.bc_data) + force_prev = np.array(nodal_forces(U_prev, ivs_prev, coordinates).at[index].get()) + disp_prev = U_prev.at[index].get() + + return 0.5*np.tensordot((force + force_prev),(disp - disp_prev), axes=1) + + def total_work_objective(self, storedState, parameters): + parameters = parameters.reshape(self.mesh.coords.shape) + 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] + ivs = p.state_data + ivs_prev = p_prev.state_data + + val += self.total_work_increment(Uu, Uu_prev, ivs, ivs_prev, p, p_prev, parameters, energyFuncs.nodal_forces) + + return val + + def total_work_gradient(self, storedState, parameters): + functionSpace = FunctionSpace.construct_function_space(self.mesh, self.quadRule) + energyFuncs = self.setup_energy_functions() + ivsUpdateInverseFuncs = MechanicsInverse.create_ivs_update_inverse_functions(functionSpace, + "plane strain", + self.materialModel) + + residualInverseFuncs = MechanicsInverse.create_path_dependent_residual_inverse_functions(energyFuncs.energy_function_coords) + + compute_df = jax.grad(self.total_work_increment, (0, 1, 2, 3, 6)) + + parameters = parameters.reshape(self.mesh.coords.shape) + gradient = np.zeros(parameters.shape) + mu = np.zeros(storedState[0][1].state_data.shape) + adjointLoad = np.zeros(storedState[0][0].shape) + + for step in reversed(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] + ivs = p.state_data + ivs_prev = p_prev.state_data + + dc_dcn = ivsUpdateInverseFuncs.ivs_update_jac_ivs_prev(self.dofManager.create_field(Uu, p.bc_data), ivs_prev) + + df_du, df_dun, df_dc, df_dcn, df_dx = compute_df(Uu, Uu_prev, ivs, ivs_prev, p, p_prev, parameters, energyFuncs.nodal_forces) + + mu += df_dc + adjointLoad -= df_du + adjointLoad -= self.dofManager.get_unknown_values(ivsUpdateInverseFuncs.ivs_update_jac_disp_vjp(self.dofManager.create_field(Uu, p.bc_data), ivs_prev, mu)) + + n = self.dofManager.get_unknown_size() + p_objective = Objective.Params(bc_data=p.bc_data, state_data=p_prev.state_data) # remember R is a function of ivs_prev + self.objective.p = p_objective + self.objective.update_precond(Uu) # update preconditioner for use in cg (will converge in 1 iteration as long as the preconditioner is not approximate) + dRdu = linalg.LinearOperator((n, n), lambda V: onp.asarray(self.objective.hessian_vec(Uu, V))) + dRdu_decomp = linalg.LinearOperator((n, n), lambda V: onp.asarray(self.objective.apply_precond(V))) + adjointVector = linalg.cg(dRdu, onp.array(adjointLoad, copy=False), rtol=1e-10, atol=0.0, M=dRdu_decomp)[0] + + gradient += df_dx + gradient += residualInverseFuncs.residual_jac_coords_vjp(Uu, p, ivs_prev, parameters, adjointVector) + gradient += ivsUpdateInverseFuncs.ivs_update_jac_coords_vjp(self.dofManager.create_field(Uu, p.bc_data), ivs_prev, parameters, mu) + + mu = np.einsum('ijk,ijkn->ijn', mu, dc_dcn) + mu += residualInverseFuncs.residual_jac_ivs_prev_vjp(Uu, p, ivs_prev, parameters, adjointVector) + mu += df_dcn + + adjointLoad = -df_dun + + return gradient.ravel() + + def compute_L2_norm_difference(self, uSteps, ivsSteps, bcsSteps, coordinates, nodal_forces): + index = (self.mesh.nodeSets['left'], 1) # arbitrarily choosing left side nodeset for reaction force + + numerator = 0.0 + denominator= 0.0 + for i in range(0, len(self.targetSteps)): + step = self.targetSteps[i] + Uu = uSteps[step] + bc_data = bcsSteps[step] + ivs = ivsSteps[step] + + U = self.dofManager.create_field(Uu, bc_data) + force = np.sum(np.array(nodal_forces(U, ivs, coordinates).at[index].get())) + + diff = force - self.targetForces[i] + numerator += diff*diff + denominator += self.targetForces[i]*self.targetForces[i] + + return np.sqrt(numerator/denominator) + + def target_curve_objective(self, storedState, parameters): + parameters = parameters.reshape(self.mesh.coords.shape) + energyFuncs = self.setup_energy_functions() + + uSteps = np.stack([storedState[i][0] for i in range(0, self.steps+1)], axis=0) + ivsSteps = np.stack([storedState[i][1].state_data for i in range(0, self.steps+1)], axis=0) + bcsSteps = np.stack([storedState[i][1].bc_data for i in range(0, self.steps+1)], axis=0) + + return self.compute_L2_norm_difference(uSteps, ivsSteps, bcsSteps, parameters, energyFuncs.nodal_forces) + + def target_curve_gradient(self, storedState, parameters): + functionSpace = FunctionSpace.construct_function_space(self.mesh, self.quadRule) + energyFuncs = self.setup_energy_functions() + ivsUpdateInverseFuncs = MechanicsInverse.create_ivs_update_inverse_functions(functionSpace, + "plane strain", + self.materialModel) + + residualInverseFuncs = MechanicsInverse.create_path_dependent_residual_inverse_functions(energyFuncs.energy_function_coords) + + parameters = parameters.reshape(self.mesh.coords.shape) + + # derivatives of F + uSteps = np.stack([storedState[i][0] for i in range(0, self.steps+1)], axis=0) + ivsSteps = np.stack([storedState[i][1].state_data for i in range(0, self.steps+1)], axis=0) + bcsSteps = np.stack([storedState[i][1].bc_data for i in range(0, self.steps+1)], axis=0) + df_du, df_dc, gradient = jax.grad(self.compute_L2_norm_difference, (0, 1, 3))(uSteps, ivsSteps, bcsSteps, parameters, energyFuncs.nodal_forces) + + mu = np.zeros(ivsSteps[0].shape) + adjointLoad = np.zeros(uSteps[0].shape) + + for step in reversed(range(1, self.steps+1)): + Uu = uSteps[step] + p = storedState[step][1] + p_prev = storedState[step-1][1] + ivs_prev = ivsSteps[step-1] + + dc_dcn = ivsUpdateInverseFuncs.ivs_update_jac_ivs_prev(self.dofManager.create_field(Uu, p.bc_data), ivs_prev) + + mu += df_dc[step] + adjointLoad -= df_du[step] + adjointLoad -= self.dofManager.get_unknown_values(ivsUpdateInverseFuncs.ivs_update_jac_disp_vjp(self.dofManager.create_field(Uu, p.bc_data), ivs_prev, mu)) + + n = self.dofManager.get_unknown_size() + p_objective = Objective.Params(bc_data=p.bc_data, state_data=p_prev.state_data) # remember R is a function of ivs_prev + self.objective.p = p_objective + self.objective.update_precond(Uu) # update preconditioner for use in cg (will converge in 1 iteration as long as the preconditioner is not approximate) + dRdu = linalg.LinearOperator((n, n), lambda V: onp.asarray(self.objective.hessian_vec(Uu, V))) + dRdu_decomp = linalg.LinearOperator((n, n), lambda V: onp.asarray(self.objective.apply_precond(V))) + adjointVector = linalg.cg(dRdu, onp.array(adjointLoad, copy=False), rtol=1e-10, atol=0.0, M=dRdu_decomp)[0] + + gradient += residualInverseFuncs.residual_jac_coords_vjp(Uu, p, ivs_prev, parameters, adjointVector) + gradient += ivsUpdateInverseFuncs.ivs_update_jac_coords_vjp(self.dofManager.create_field(Uu, p.bc_data), ivs_prev, parameters, mu) + + mu = np.einsum('ijk,ijkn->ijn', mu, dc_dcn) + mu += residualInverseFuncs.residual_jac_ivs_prev_vjp(Uu, p, ivs_prev, parameters, adjointVector) + + adjointLoad = np.zeros(storedState[0][0].shape) + + return gradient.ravel() + + def test_total_work_gradient_with_adjoint_solve(self): + self.compute_objective_function = self.total_work_objective + self.compute_gradient = self.total_work_gradient + + stepSize = 1e-7 + + error = self.compute_finite_difference_error(stepSize, self.initialMesh.coords.ravel()) + self.assertLessEqual(error, 1e-7) + + def test_target_curve_gradient_with_adjoint_solve(self): + self.compute_objective_function = self.target_curve_objective + self.compute_gradient = self.target_curve_gradient + + self.targetSteps = [1, 2] + self.targetForces = [4.5, 5.5] # [4.542013626078756, 5.7673988583067555] actual forces + + stepSize = 1e-8 + + error = self.compute_finite_difference_error(stepSize, self.initialMesh.coords.ravel()) + self.assertLessEqual(error, 1e-6) + + + +if __name__ == '__main__': + unittest.main() diff --git "a/optimism/test/inverse/test_J2Plastic_gradient_checks.py\357\200\272Zone.Identifier" "b/optimism/test/inverse/test_J2Plastic_gradient_checks.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git a/optimism/test/inverse/test_J2Plastic_inverse.py b/optimism/test/inverse/test_J2Plastic_inverse.py new file mode 100644 index 00000000..c38fd641 --- /dev/null +++ b/optimism/test/inverse/test_J2Plastic_inverse.py @@ -0,0 +1,277 @@ +import unittest + +import jax +import jax.numpy as np + +from . import EquationSolver_Immersed_2 as EqSolver +from optimism import FunctionSpace +from optimism.material import J2Plastic as J2 +from optimism import Mechanics +from optimism import Mesh +from optimism import Objective +from optimism import QuadratureRule +from optimism import TensorMath + +from ..TestFixture import TestFixture +from ..MeshFixture import MeshFixture + +from optimism.inverse import MechanicsInverse + +def make_disp_grad_from_strain(strain): + return linalg.expm(strain) - np.identity(3) + +def small_strain_linear_hardening_analytic_gradients(disp_grad, state, state_previous, props): + sqrt32 = np.sqrt(3./2.) + mu = 0.5*props['elastic modulus']/(1.0 + props['poisson ratio']) + c0 = 3.*mu + props['hardening modulus'] + + vol_projection = np.zeros((9,9)).at[(0,0,0,4,4,4,8,8,8),(0,4,8,0,4,8,0,4,8)].set(1.) + dev_projection = np.eye(9,9) - vol_projection/3. + sym_projection = np.zeros((9,9)).at[(1,1,2,2,3,3,5,5,6,6,7,7),(1,3,2,6,1,3,5,7,2,6,5,7)].set(0.5).at[(0,4,8),(0,4,8)].set(1.) + + eqps = state[0] + deqps = eqps - state_previous[0] + + eps_pn = state_previous[1:] + eps = TensorMath.sym(disp_grad) + eps_e_tr = eps - eps_pn.reshape((3,3)) + + dev_eps_e_tr = TensorMath.dev(eps_e_tr) + dev_eps_e_tr_squared = np.tensordot(dev_eps_e_tr, dev_eps_e_tr) + n_tr = 1./np.sqrt(dev_eps_e_tr_squared) * dev_eps_e_tr + + deqps_deps = (2.*mu * sqrt32/ c0)*n_tr + deps_p_deps = (sqrt32 * deqps / np.sqrt(dev_eps_e_tr_squared)) * (np.tensordot(dev_projection,sym_projection,axes=1) - np.kron(n_tr.ravel(),n_tr.ravel()).reshape(9,9)) + deps_p_deps += sqrt32 * np.kron((2.*mu * sqrt32/ c0) * n_tr.ravel(), n_tr.ravel()).reshape(9,9) + + deqps_deqps_n = 3.*mu / c0 + deqps_deps_p_n = -(2.*mu * sqrt32/ c0)*n_tr + deqps_dc_n = np.append(deqps_deqps_n, deqps_deps_p_n.ravel()) + + deps_p_deqps_n = (3.*mu / c0 - 1.) * sqrt32 * n_tr + deps_p_deps_p_n = np.eye(9,9) - (sqrt32 * deqps / np.sqrt(dev_eps_e_tr_squared)) * (dev_projection - np.kron(n_tr.ravel(),n_tr.ravel()).reshape(9,9)) + deps_p_deps_p_n -= sqrt32 * np.kron((2.*mu * sqrt32/ c0) * n_tr.ravel(), n_tr.ravel()).reshape(9,9) + deps_p_dc_n = np.hstack((deps_p_deqps_n.ravel()[:,None],deps_p_deps_p_n)) + + return np.vstack((deqps_deps.ravel(), deps_p_deps)), np.vstack((deqps_dc_n, deps_p_dc_n)) + + +class J2MaterialPointUpdateGradsFixture(TestFixture): + def setUp(self): + E = 100.0 + poisson = 0.321 + Y0 = 0.2*E + H = 1.0e-2*E + + self.props = {'elastic modulus': E, + 'poisson ratio': poisson, + 'yield strength': Y0, + 'kinematics': 'small deformations', + 'hardening model': 'linear', + 'hardening modulus': H} + + materialModel = J2.create_material_model_functions(self.props) + + self.compute_state_new = jax.jit(materialModel.compute_state_new) + self.compute_initial_state = materialModel.compute_initial_state + + self.compute_state_new_derivs = jax.jit(jax.jacfwd(self.compute_state_new, (0, 1))) + + def test_jax_computation_of_state_derivs_at_elastic_step(self): + dt = 1.0 + + initial_state = self.compute_initial_state() + dc_dugrad, dc_dc_n = self.compute_state_new_derivs(np.zeros((3,3)), initial_state, dt) + + # Check derivative sizes + self.assertEqual(dc_dugrad.shape, (10,3,3)) + self.assertEqual(dc_dc_n.shape, (10,10)) + + # check derivatives w.r.t. displacement grad + self.assertArrayNear(dc_dugrad.ravel(), np.zeros((10,3,3)).ravel(), 12) + + # check derivatives w.r.t. previous step internal variables + self.assertArrayNear(dc_dc_n.ravel(), np.eye(10).ravel(), 12) + + def test_jax_computation_of_state_derivs_at_plastic_step(self): + dt = 1.0 + # get random displacement gradient + key = jax.random.PRNGKey(0) + dispGrad = jax.random.uniform(key, (3, 3)) + initial_state = self.compute_initial_state() + + state = self.compute_state_new(dispGrad, initial_state, dt) + dc_dugrad, dc_dc_n = self.compute_state_new_derivs(dispGrad, initial_state, dt) + + # Compute data for gold values (assuming small strain and linear kinematics) + dc_dugrad_gold, dc_dc_n_gold = small_strain_linear_hardening_analytic_gradients(dispGrad, state, initial_state, self.props) + + # Check derivative sizes + self.assertEqual(dc_dugrad.shape, (10,3,3)) + self.assertEqual(dc_dc_n.shape, (10,10)) + + # check derivatives w.r.t. displacement grad + self.assertArrayNear(dc_dugrad[0,:,:].ravel(), dc_dugrad_gold[0,:].ravel(), 12) + self.assertArrayNear(dc_dugrad[1:,:,:].ravel(), dc_dugrad_gold[1:,:].ravel() , 12) + + # check derivatives w.r.t. previous step internal variables + self.assertNear(dc_dc_n[0,0], dc_dc_n_gold[0,0], 12) + self.assertArrayNear(dc_dc_n[0,1:], dc_dc_n_gold[0,1:], 12) + self.assertArrayNear(dc_dc_n[1:,0], dc_dc_n_gold[1:,0], 12) + self.assertArrayNear(dc_dc_n[1:,1:].ravel(), dc_dc_n_gold[1:,1:].ravel() , 12) + + +class J2GlobalMeshUpdateGradsFixture(MeshFixture): + @classmethod + def setUpClass(cls): + dispGrad0 = np.array([[0.4, -0.2], + [-0.04, 0.68]]) + mf = MeshFixture() + cls.mesh, cls.U = mf.create_mesh_and_disp(4,4,[0.,1.],[0.,1.], + lambda x: dispGrad0@x) + + E = 100.0 + poisson = 0.321 + H = 1e-2 * E + Y0 = 0.3 * E + + cls.props = {'elastic modulus': E, + 'poisson ratio': poisson, + 'yield strength': Y0, + 'kinematics': 'small deformations', + 'hardening model': 'linear', + 'hardening modulus': H} + + cls.materialModel = J2.create_material_model_functions(cls.props) + cls.quadRule = QuadratureRule.create_quadrature_rule_on_triangle(degree=1) + cls.fs = FunctionSpace.construct_function_space(cls.mesh, cls.quadRule) + cls.mechFuncs = Mechanics.create_mechanics_functions(cls.fs, + "plane strain", + cls.materialModel) + cls.ivs_prev = cls.mechFuncs.compute_initial_state() + + EBCs = [FunctionSpace.EssentialBC(nodeSet='all_boundary', component=0), + FunctionSpace.EssentialBC(nodeSet='all_boundary', component=1)] + cls.dofManager = FunctionSpace.DofManager(cls.fs, 2, EBCs) + cls.Ubc = cls.dofManager.get_bc_values(cls.U) + + p = Objective.Params(None, cls.ivs_prev, None, None, None) + UuGuess = 0.0*cls.dofManager.get_unknown_values(cls.U) + + def compute_energy(Uu, p): + U = cls.dofManager.create_field(Uu, cls.Ubc) + internalVariables = p[1] + return cls.mechFuncs.compute_strain_energy(U, internalVariables) + + objective = Objective.Objective(compute_energy, UuGuess, p) + cls.Uu, solverSuccess = EqSolver.nonlinear_equation_solve(objective, UuGuess, p, EqSolver.get_settings(), useWarmStart=False) + U = cls.dofManager.create_field(cls.Uu, cls.Ubc) + cls.ivs = cls.mechFuncs.compute_updated_internal_variables(U, cls.ivs_prev) + + def test_state_derivs_at_elastic_step(self): + + def update_internal_vars_test(Uu, ivs_prev): + U = self.dofManager.create_field(Uu) + ivs = self.mechFuncs.compute_updated_internal_variables(U, ivs_prev) + return ivs + + Uu = 0.0*self.dofManager.get_unknown_values(self.U) + + update_internal_variables_derivs = jax.jacfwd(update_internal_vars_test, (0,1)) + dc_du, dc_dc_n = update_internal_variables_derivs(Uu, self.ivs_prev) + + nElems = Mesh.num_elements(self.mesh) + nQpsPerElem = QuadratureRule.len(self.quadRule) + nIntVars = 10 + nFreeDofs = Uu.shape[0] + + self.assertEqual(dc_du.shape, (nElems,nQpsPerElem,nIntVars,nFreeDofs)) + self.assertEqual(dc_dc_n.shape, (nElems,nQpsPerElem,nIntVars,nElems,nQpsPerElem,nIntVars)) + + for i in range(0,nElems): + for j in range(0,nQpsPerElem): + self.assertArrayNear(dc_du[i,j,:,:].ravel(), np.zeros((nIntVars,nFreeDofs)).ravel(), 12) + self.assertArrayNear(dc_dc_n[i,j,:,i,j,:].ravel(), np.eye(nIntVars).ravel(), 12) + + def test_state_derivs_at_plastic_step(self): + + def update_internal_vars_test(U, ivs_prev): + ivs = self.mechFuncs.compute_updated_internal_variables(U, ivs_prev) + return ivs + + update_internal_variables_derivs = jax.jacfwd(update_internal_vars_test, (0,1)) + + U = self.dofManager.create_field(self.Uu, self.Ubc) + dc_du, dc_dc_n = update_internal_variables_derivs(U, self.ivs_prev) + + nElems = Mesh.num_elements(self.mesh) + nQpsPerElem = QuadratureRule.len(self.quadRule) + nIntVars = 10 + nDims = 2 + nNodes = Mesh.num_nodes(self.mesh) + + self.assertEqual(dc_du.shape, (nElems,nQpsPerElem,nIntVars,nNodes,nDims)) + self.assertEqual(dc_dc_n.shape, (nElems,nQpsPerElem,nIntVars,nElems,nQpsPerElem,nIntVars)) + + for i in range(0,nElems): + for j in range(0,nQpsPerElem): + state = self.ivs[i,j,:] + initial_state = self.ivs_prev[i,j,:] + + conn = self.mesh.conns[i] + Uele = U[conn] + shapeGrads = self.fs.shapeGrads[i,j,:,:] + dispGrad = TensorMath.tensor_2D_to_3D(np.tensordot(Uele,shapeGrads,axes=[0,0])) + Be_mat = np.zeros((9,6)).\ + at[(0,3),(0,1)].set(shapeGrads[0,0]).at[(1,4),(0,1)].set(shapeGrads[0,1]).\ + at[(0,3),(2,3)].set(shapeGrads[1,0]).at[(1,4),(2,3)].set(shapeGrads[1,1]).\ + at[(0,3),(4,5)].set(shapeGrads[2,0]).at[(1,4),(4,5)].set(shapeGrads[2,1]) + + dc_dugrad_gold, dc_dc_n_gold = small_strain_linear_hardening_analytic_gradients(dispGrad, state, initial_state, self.props) + + dc_duele_gold = np.tensordot(dc_dugrad_gold, Be_mat, axes=1) + dc_er_du = dc_du[i,j,:,:,:] + + self.assertArrayNear(dc_er_du[:,conn,:].ravel(), dc_duele_gold.reshape(10,3,2).ravel(), 10) + self.assertArrayNear(np.delete(dc_er_du, conn, axis=1).ravel(), np.zeros((10,nNodes-conn.shape[0],2)).ravel(), 10) + + for p in range(0,nElems): + for q in range(0,nQpsPerElem): + if(i == p and j == q): + self.assertArrayNear(dc_dc_n[i,j,:,i,j,:].ravel(), dc_dc_n_gold.ravel(), 10) + else: + self.assertArrayNear(dc_dc_n[i,j,:,p,q,:].ravel(), np.zeros((nIntVars,nIntVars)).ravel(), 10) + + def test_state_derivs_computed_locally_at_plastic_step(self): + + ivsUpdateInverseFuncs = MechanicsInverse.create_ivs_update_inverse_functions(self.fs, + "plane strain", + self.materialModel) + + U = self.dofManager.create_field(self.Uu, self.Ubc) + dc_dc_n = ivsUpdateInverseFuncs.ivs_update_jac_ivs_prev(U, self.ivs_prev) + + nElems = Mesh.num_elements(self.mesh) + nQpsPerElem = QuadratureRule.len(self.quadRule) + nIntVars = 10 + + self.assertEqual(dc_dc_n.shape, (nElems,nQpsPerElem,nIntVars,nIntVars)) + + for i in range(0,nElems): + for j in range(0,nQpsPerElem): + state = self.ivs[i,j,:] + initial_state = self.ivs_prev[i,j,:] + + conn = self.mesh.conns[i] + Uele = U[conn] + shapeGrads = self.fs.shapeGrads[i,j,:,:] + dispGrad = TensorMath.tensor_2D_to_3D(np.tensordot(Uele,shapeGrads,axes=[0,0])) + + _, dc_dc_n_gold = small_strain_linear_hardening_analytic_gradients(dispGrad, state, initial_state, self.props) + + self.assertArrayNear(dc_dc_n[i,j,:,:].ravel(), dc_dc_n_gold.ravel(), 10) + + + +if __name__ == '__main__': + unittest.main() diff --git "a/optimism/test/inverse/test_J2Plastic_inverse.py\357\200\272Zone.Identifier" "b/optimism/test/inverse/test_J2Plastic_inverse.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git a/optimism/test/material/__init__.py b/optimism/test/material/__init__.py new file mode 100644 index 00000000..e69de29b diff --git "a/optimism/test/material/__init__.py\357\200\272Zone.Identifier" "b/optimism/test/material/__init__.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git a/optimism/test/material/plotJ2Plastic.py b/optimism/test/material/plotJ2Plastic.py new file mode 100644 index 00000000..89067944 --- /dev/null +++ b/optimism/test/material/plotJ2Plastic.py @@ -0,0 +1,31 @@ +from matplotlib import pyplot as plt + +from optimism.JaxConfig import * +from optimism.material.MaterialPointUniaxialSimulator import MaterialPointUniaxialSimulator +from optimism.material import J2Plastic as J2 + + +def plot_J2_uniaxial(): + E = 200.0e3 + nu = 0.25 + Y0 = 350.0 + H = 1e-1*E + Ysat = 800.0 + eps0 = 0.01 + props = {'elastic modulus': E, + 'poisson ratio': nu, + 'yield strength': Y0, + 'hardening model': 'voce', + 'hardening modulus': H, + 'saturation strength': Ysat, + 'reference plastic strain': eps0} + materialModel = J2.create_material_model_functions(props) + maxStrain = 20.0*Y0/E + strainRate = 1.0 + simulator = MaterialPointUniaxialSimulator(materialModel, maxStrain, strainRate, steps=20) + out = simulator.run() + plt.plot(out.strainHistory, out.stressHistory, marker='o') + plt.show() + +if __name__ == "__main__": + plot_J2_uniaxial() diff --git "a/optimism/test/material/plotJ2Plastic.py\357\200\272Zone.Identifier" "b/optimism/test/material/plotJ2Plastic.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git a/optimism/test/material/test_Gent.py b/optimism/test/material/test_Gent.py new file mode 100644 index 00000000..e7aaa347 --- /dev/null +++ b/optimism/test/material/test_Gent.py @@ -0,0 +1,88 @@ +import numpy as onp +from scipy.spatial.transform import Rotation +import unittest + +import jax +import jax.numpy as np + +from optimism.material import Gent +from .. import TestFixture + + +class TestGentMaterial(TestFixture.TestFixture): + def setUp(self): + self.kappa = 100.0 + self.mu = 10.0 + self.Jm = 4.0 # pick large enough to avoid singularity in tests + + properties = {"bulk modulus": self.kappa, + "shear modulus": self.mu, + "Jm parameter": self.Jm} + + self.material = Gent.create_material_functions(properties) + + self.internalVariables = self.material.compute_initial_state() + + self.dt = 0.0 + + + def test_zero_point(self): + dispGrad = np.zeros((3, 3)) + W = self.material.compute_energy_density(dispGrad, self.internalVariables, self.dt) + self.assertLessEqual(W, np.linalg.norm(dispGrad)*1e-10) + + + def test_frame_indifference(self): + # generate a random displacement gradient + key = jax.random.PRNGKey(1) + dispGrad = jax.random.uniform(key, (3, 3)) + + + W = self.material.compute_energy_density(dispGrad, self.internalVariables, self.dt) + for i in range(10): + Q = Rotation.random(random_state=i).as_matrix() + dispGradTransformed = Q@(dispGrad + np.identity(3)) - np.identity(3) + WStar = self.material.compute_energy_density(dispGradTransformed, self.internalVariables, self.dt) + self.assertAlmostEqual(W, WStar, 12) + + + def test_correspondence_with_linear_elasticity(self): + zero = np.zeros((3, 3)) + C = jax.hessian(self.material.compute_energy_density)(zero, self.internalVariables, self.dt) + + lam = self.kappa - 2/3*self.mu + + CLinear = onp.zeros((3, 3, 3, 3)) + for i in range(3): + for j in range(3): + for k in range(3): + for l in range(3): + CLinear[i, j, k, l] += self.mu if i == k and j == l else 0 + CLinear[i, j, k, l] += self.mu if i == l and j == k else 0 + CLinear[i, j, k, l] += lam if i == j and k == l else 0 + + self.assertArrayNear(C, CLinear, 12) + + + def test_finite_extensibility(self): + # incompressible uniaxial extension + # find stretch such that the strain energy just reaches the singularity. + lockStretchCandidates = onp.roots([1.0, 0.0, -(self.Jm + 3), 2.0]) + lockStretch = onp.max(lockStretchCandidates) + stretch = lockStretch*(1 + 1e-6) # account for finite precision of root finder + I1 = stretch**2 + 2/stretch + self.assertGreater(I1 - 3, self.Jm) + + # Check that energy is indeed infinite + # (actually nan, since it produces a negative argument to a logarithm) + F = np.diag(np.array([stretch, 1/np.sqrt(stretch), 1/np.sqrt(stretch)])) + W = self.material.compute_energy_density(F - np.identity(3), self.internalVariables, self.dt) + self.assertTrue(np.isnan(W)) + + stretch = lockStretch*(1 - 1e-6) + F = np.diag(np.array([stretch, 1/np.sqrt(stretch), 1/np.sqrt(stretch)])) + W = self.material.compute_energy_density(F - np.identity(3), self.internalVariables, self.dt) + self.assertFalse(np.isnan(W)) + +if __name__ == "__main__": + unittest.main() \ No newline at end of file diff --git "a/optimism/test/material/test_Gent.py\357\200\272Zone.Identifier" "b/optimism/test/material/test_Gent.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git a/optimism/test/material/test_Hardening.py b/optimism/test/material/test_Hardening.py new file mode 100644 index 00000000..3ff6412a --- /dev/null +++ b/optimism/test/material/test_Hardening.py @@ -0,0 +1,83 @@ +import unittest + +from optimism.JaxConfig import * +from optimism.material import J2Plastic +from optimism.material import Hardening +from ..TestFixture import TestFixture + + +E = 100.0 +nu = 0.321 +Y0 = E*0.003 + +class VoceHardeningTestFixture(TestFixture): + def setUp(self): + self.Ysat = 4.0*Y0 + self.eps0 = 0.05 + self.props = {'hardening model': 'voce', + 'elastic modulus': E, + 'poisson ratio': nu, + 'yield strength': Y0, + 'saturation strength': self.Ysat, + 'reference plastic strain': self.eps0} + self.plasticEnergy, self.flowStrength = Hardening.create_hardening_model(self.props) + + + def test_voce_hardening_zero_point(self): + eqps = 0.0 + Wp = self.plasticEnergy(eqps, eqpsOld=0.0, dt=0.0) + self.assertLess(Wp, 1e-14) + + + def test_voce_hardening_yield_strength(self): + eqps = 0.0 + Y = self.flowStrength(eqps, eqpsOld=0.0, dt=0.0) + self.assertAlmostEqual(Y, Y0, 12) + + + def test_voce_hardening_saturates_to_correct_value(self): + eqps = 15.0*self.eps0 + Y = self.flowStrength(eqps, eqpsOld=0.0, dt=0.0) + self.assertAlmostEqual(Y, self.Ysat, 5) + + + +class PowerLawHardeningTestFixture(TestFixture): + def setUp(self): + self.n = 4.0 + self.eps0 = 0.05 + self.props = {'hardening model': 'power law', + 'elastic modulus': E, + 'poisson ratio': nu, + 'yield strength': Y0, + 'hardening exponent': self.n, + 'reference plastic strain': self.eps0} + self.plasticEnergy, self.flowStrength = Hardening.create_hardening_model(self.props) + + + def test_power_law_hardening_zero_point(self): + Wp = self.plasticEnergy(eqps=0.0, eqpsOld=0.0, dt=0.0) + self.assertLess(Wp, 1e-14) + + + def test_power_law_hardening_yield_strength(self): + eqps = 0.0 + Y = self.flowStrength(eqps, eqpsOld=0.0, dt=0.0) + self.assertAlmostEqual(Y, Y0, 12) + + + def test_power_law_strength_increases(self): + eqps = 5.0*self.eps0 + Y = self.flowStrength(eqps, eqpsOld=0.0, dt=0.0) + self.assertGreater(Y, Y0) + + + def test_power_law_hardening_slope_is_finite_at_origin(self): + hardeningRate = jacfwd(self.flowStrength) + eqps = 0.0 + H = hardeningRate(eqps, eqpsOld=0.0, dt=0.0) + self.assertTrue(np.isfinite(H)) + + +if __name__ == '__main__': + unittest.main() diff --git "a/optimism/test/material/test_Hardening.py\357\200\272Zone.Identifier" "b/optimism/test/material/test_Hardening.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git a/optimism/test/material/test_HyperVisco.py b/optimism/test/material/test_HyperVisco.py new file mode 100644 index 00000000..19c3a547 --- /dev/null +++ b/optimism/test/material/test_HyperVisco.py @@ -0,0 +1,139 @@ +import unittest + +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 ..TestFixture import TestFixture +from optimism.material import MaterialUniaxialSimulator +from optimism.TensorMath import deviator +from optimism.TensorMath import log_symm + +plotting = False + +class HyperViscoModelFixture(TestFixture): + def setUp(self): + + G_eq = 0.855 # MPa + K_eq = 1000*G_eq # MPa + 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.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 = jax.jit(materialModel.compute_state_new) + self.compute_initial_state = materialModel.compute_initial_state + self.compute_material_qoi = jax.jit(materialModel.compute_material_qoi) + +class HyperViscoUniaxialStrain(HyperViscoModelFixture): + def test_loading_only(self): + strain_rate = 1.0e-2 + total_time = 100.0 + n_steps = 100 + dt = total_time / n_steps + times = np.linspace(0.0, total_time, n_steps) + Fs = jax.vmap( + lambda t: np.array( + [[np.exp(strain_rate * t), 0.0, 0.0], + [0.0, 1.0, 0.0], + [0.0, 0.0, 1.0]] + ) + )(times) + state_old = self.compute_initial_state() + energies = np.zeros(n_steps) + states = np.zeros((n_steps, state_old.shape[0])) + + # 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) + 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)) + + e_e_11 = strain_rate * times - e_v_11 + e_e_22 = 0.5 * e_v_11 + + Ee_analytic = jax.vmap( + lambda e_11, e_22: np.array( + [[e_11, 0., 0.], + [0., e_22, 0.], + [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) + + # 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) + + 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.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') + plt.plot(times, e_e_22, linestyle='--', color='red') + plt.xlabel('Time (s)') + plt.ylabel('Viscous Elastic Strain') + 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.legend(["Numerical", "Analytic"], loc=0, frameon=True) + plt.savefig('times_dissipation.png') + + + +if __name__ == '__main__': + unittest.main() diff --git "a/optimism/test/material/test_HyperVisco.py\357\200\272Zone.Identifier" "b/optimism/test/material/test_HyperVisco.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git a/optimism/test/material/test_J2Plastic.py b/optimism/test/material/test_J2Plastic.py new file mode 100644 index 00000000..e50e8e10 --- /dev/null +++ b/optimism/test/material/test_J2Plastic.py @@ -0,0 +1,332 @@ +import unittest + +import jax +import jax.numpy as np +from jax.scipy import linalg +from matplotlib import pyplot as plt + +from . import EquationSolver_Immersed_2 as EqSolver +from optimism import FunctionSpace +from optimism.material import J2Plastic as J2 +from optimism.material import MaterialUniaxialSimulator +from optimism import Mechanics +from optimism import Mesh +from optimism import Objective +from optimism import QuadratureRule +from ..TestFixture import TestFixture +from ..MeshFixture import MeshFixture + + +plotting=False + + +def make_disp_grad_from_strain(strain): + return linalg.expm(strain) - np.identity(3) + + +class GradOfPlasticityModelFixture(TestFixture): + def setUp(self): + + E = 100.0 + poisson = 0.321 + Y0 = 0.3*E + H = 1.0e-2*E + + self.props = {'elastic modulus': E, + 'poisson ratio': poisson, + 'yield strength': Y0, + 'hardening model': 'linear', + 'hardening modulus': H} + + materialModel = J2.create_material_model_functions(self.props) + + self.energy_density = jax.jit(materialModel.compute_energy_density) + self.stress_func = jax.jit(jax.grad(self.energy_density, 0)) + self.compute_state_new = materialModel.compute_state_new + self.tangents_func = jax.hessian(self.energy_density) + self.compute_initial_state = materialModel.compute_initial_state + + def test_zero_point(self): + dispGrad = np.zeros((3,3)) + state = self.compute_initial_state() + dt = 1.0 + + energy = self.energy_density(dispGrad, state, dt) + self.assertNear(energy, 0.0, 12) + + stress = self.stress_func(dispGrad, state, dt) + self.assertArrayNear(stress, np.zeros((3,3)), 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) + + + def test_elastic_strain_path(self): + strain = np.zeros((3,3)) + strain_inc = 1e-6 + stateOld = self.compute_initial_state() + dt = 1.0 + + strainHistory = [] + stressHistory = [] + eqpsHistory = [] + energyHistory = [] + for i in range(10): + dispGrad = make_disp_grad_from_strain(strain) + F = dispGrad + np.identity(3) + energy = self.energy_density(dispGrad, stateOld, dt) + strainHistory.append(strain[0,0]) + stateNew = self.compute_state_new(dispGrad, stateOld, dt) + stressNew = self.stress_func(dispGrad, stateNew, dt) @ F.T + stressHistory.append(stressNew[0,0]) + eqpsHistory.append(stateNew[J2.EQPS]) + energyHistory.append(energy) + + stateOld = stateNew + strain = dispGrad.at[0,0].add(strain_inc) + + E = self.props['elastic modulus'] + nu = self.props['poisson ratio'] + stressNewExp = E * (1.0 - nu) / ((1.0 + nu) * (1.0 - 2.0*nu)) * strainHistory[-1] + self.assertNear(stressHistory[-1], stressNewExp, 12) + self.assertArrayNear(eqpsHistory, np.zeros(len(eqpsHistory)), 12) + + + def test_plastic_strain_path(self): + strain = np.zeros((3,3)) + strain_inc = 0.2 + stateOld = self.compute_initial_state() + dt = 1.0 + + strainHistory = [] + stressHistory = [] + tangentsHistory = [] + eqpsHistory = [] + energyHistory = [] + for i in range(2): + strain = strain.at[0,0].add(strain_inc) + dispGrad = make_disp_grad_from_strain(strain) + F = dispGrad + np.identity(3) + energy = self.energy_density(dispGrad, stateOld, dt) + tangentsNew = self.tangents_func(dispGrad, stateOld, dt) + stateNew = self.compute_state_new(dispGrad, stateOld, dt) + stressNew = self.stress_func(dispGrad, stateOld, dt) @ F.T + strainHistory.append(strain[0,0]) + stressHistory.append(stressNew[0,0]) + tangentsHistory.append(tangentsNew[0,0,0,0]) + eqpsHistory.append(stateNew[J2.EQPS]) + energyHistory.append(energy) + + stateOld = stateNew + + if plotting: + plt.figure() + plt.plot(strainHistory, energyHistory, marker='o', fillstyle='none') + plt.xlabel('strain') + plt.ylabel('potential density') + + plt.figure() + plt.plot(strainHistory, stressHistory, marker='o', fillstyle='none') + plt.xlabel('strain') + plt.ylabel('stress') + + plt.figure() + plt.plot(strainHistory, eqpsHistory, marker='o') + plt.xlabel('strain') + plt.ylabel('Eq Pl Strain') + + plt.figure() + plt.plot(strainHistory, tangentsHistory, marker='o') + plt.xlabel('strain') + plt.ylabel('Tangent modulus') + + plt.show() + + E = self.props['elastic modulus'] + nu = self.props['poisson ratio'] + mu = 0.5*E/(1.0 + nu) + lam = E*nu/(1.0 + nu)/(1.0 - 2.0*nu) + Y0 = self.props['yield strength'] + H = self.props['hardening modulus'] + strainEnd = strainHistory[-1] + # for solutions, refer to jax-fem/papers/plane_strain_unit_test.pdf + eqpsExp = [max(0.0, (2.0 * i * mu - Y0) / (3.0*mu + H) ) for i in strainHistory] + stressNewExp = (2.0 * mu * Y0 + \ + 2.0 * strainEnd * pow(mu,2) + \ + strainEnd * H * lam + \ + 2.0 * strainEnd * H * mu + \ + 3.0 * strainEnd * lam * mu) \ + / (3.0 * mu + H) + + self.assertNear( stressHistory[-1], stressNewExp, 11 ) + self.assertArrayNear(eqpsHistory, eqpsExp, 12) + + +class J2UpdateFixture(TestFixture): + def setUp(self): + E = 100.0 + poisson = 0.321 + Y0 = 0.3*E + eps0 = Y0/E + n = 4 + + self.props = {'elastic modulus': E, + 'poisson ratio': poisson, + 'yield strength': Y0, + 'hardening model': 'power law', + 'hardening exponent': n, + 'reference plastic strain': eps0} + + materialModel = J2.create_material_model_functions(self.props) + + self.compute_state_new = jax.jit(materialModel.compute_state_new) + self.compute_initial_state = materialModel.compute_initial_state + + def test_update_only_happens_once(self): + dt = 1.0 + + # get 10,000 random displacement gradients + key = jax.random.PRNGKey(0) + dispGrads = jax.random.uniform(key, (10000, 3, 3)) + + # get 10,000 copies of the initial state + states = np.tile(self.compute_initial_state(), (10000, 1)) + + states_new = jax.vmap(self.compute_state_new, (0, 0, None))(dispGrads, states, dt) + + states_should_be_unchanged = jax.vmap(self.compute_state_new, (0, 0, None))(dispGrads, states_new, dt) + self.assertArrayEqual(states_new, states_should_be_unchanged) + + +class J2PlasticUniaxial(TestFixture): + + def setUp(self): + E = 100.0e3 + nu = 0.25 + Y0 = 30.0 + H = E/200 + + properties = {'elastic modulus': E, + 'poisson ratio': nu, + 'yield strength': Y0, + 'kinematics': 'large deformations', + 'hardening model': 'linear', + 'hardening modulus': H} + self.E = E + self.Y0 = Y0 + self.H = H + self.mat = J2.create_material_model_functions(properties) + + + def test_uniaxial(self): + strainRate = 1e-3 + + def constant_true_strain_rate(t): + return np.expm1(strainRate*t) + + maxTime = 20.0 + + uniaxial = MaterialUniaxialSimulator.run(self.mat, constant_true_strain_rate, maxTime, steps=10) + + logStrainHistory = np.log(1.0 + uniaxial.strainHistory[:,0,0]) + yieldStrain = self.Y0/self.E + exact = np.where(logStrainHistory < yieldStrain, + self.E*logStrainHistory, + self.E/(self.E + self.H)*(self.H*logStrainHistory + self.Y0)) + + # convert Piola stress output to Kirchhoff stress + I = np.identity(3) + kirchhoffStressHistory = jax.vmap(lambda H, P: P@(H + I).T)(uniaxial.strainHistory, uniaxial.stressHistory) + + self.assertArrayNear(kirchhoffStressHistory[:,0,0], exact, 2) + + +class PlasticityOnMesh(MeshFixture): + + def test_plasticity_with_mesh(self): + + dispGrad0 = np.array([[0.4, -0.2], + [-0.04, 0.68]]) + mesh, U = self.create_mesh_and_disp(4,4,[0.,1.],[0.,1.], + lambda x: dispGrad0@x) + + E = 100.0 + poisson = 0.321 + H = 1e-2 * E + Y0 = 0.3 * E + + props = {'elastic modulus': E, + 'poisson ratio': poisson, + 'yield strength': Y0, + 'hardening model': 'linear', + 'hardening modulus': H} + + materialModel = J2.create_material_model_functions(props) + + quadRule = QuadratureRule.create_quadrature_rule_on_triangle(degree=1) + fs = FunctionSpace.construct_function_space(mesh, quadRule) + + mechFuncs = Mechanics.create_mechanics_functions(fs, + "plane strain", + materialModel) + + EBCs = [FunctionSpace.EssentialBC(nodeSet='all_boundary', component=0), + FunctionSpace.EssentialBC(nodeSet='all_boundary', component=1)] + dofManager = FunctionSpace.DofManager(fs, 2, EBCs) + + Ubc = dofManager.get_bc_values(U) + + nElems = Mesh.num_elements(mesh) + nQpsPerElem = QuadratureRule.len(quadRule) + internalVariables = mechFuncs.compute_initial_state() + + tOld = 0.0 + t = 1.0 + p = Objective.Params(None, internalVariables, None, None, np.array([t, tOld])) + + + def compute_energy(Uu, p): + U = dofManager.create_field(Uu, Ubc) + internalVariables = p[1] + dt = p.time[0] - p.time[1] + return mechFuncs.compute_strain_energy(U, internalVariables, dt) + + UuGuess = 0.0*dofManager.get_unknown_values(U) + + objective = Objective.Objective(compute_energy, UuGuess, p) + + Uu, solverSuccess = EqSolver.nonlinear_equation_solve(objective, UuGuess, p, EqSolver.get_settings(), useWarmStart=False) + self.assertTrue(solverSuccess) + + U = dofManager.create_field(Uu, Ubc) + + dispGrads = FunctionSpace.compute_field_gradient(fs, U) + for dg in dispGrads.reshape(nElems*nQpsPerElem, 2, 2): + self.assertArrayNear(dg, dispGrad0, 11) + + internalVariablesNew = mechFuncs.compute_updated_internal_variables(U, p[1], t - tOld) + + # check to make sure plastic strain evolved + # if this fails, make the applied displacement grad bigger + eqpsField = internalVariablesNew[:,:,J2.EQPS] + self.assertTrue(eqpsField[0] > 1e-8) + + +if __name__ == '__main__': + unittest.main() diff --git "a/optimism/test/material/test_J2Plastic.py\357\200\272Zone.Identifier" "b/optimism/test/material/test_J2Plastic.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git a/optimism/test/material/test_LinearElastic.py b/optimism/test/material/test_LinearElastic.py new file mode 100644 index 00000000..982319e4 --- /dev/null +++ b/optimism/test/material/test_LinearElastic.py @@ -0,0 +1,55 @@ +from scipy.spatial.transform import Rotation +import unittest + +import jax +import jax.numpy as np + +from optimism.material import LinearElastic +from .. import TestFixture + +class TestLinearElasticMaterial(TestFixture.TestFixture): + def setUp(self): + self.E = 10.0 + self.nu = 0.25 + + properties = {"elastic modulus": self.E, + "poisson ratio": self.nu, + "strain measure": "logarithmic"} + + self.material = LinearElastic.create_material_model_functions(properties) + + self.internalVariables = self.material.compute_initial_state() + self.dt = 0.0 + + + def test_zero_point(self): + dispGrad = np.zeros((3, 3)) + W = self.material.compute_energy_density(dispGrad, self.internalVariables, self.dt) + self.assertLessEqual(W, np.linalg.norm(dispGrad)*1e-10) + + + def test_finite_deformation_frame_indifference(self): + # generate a random displacement gradient + key = jax.random.PRNGKey(1) + dispGrad = jax.random.uniform(key, (3, 3)) + + W = self.material.compute_energy_density(dispGrad, self.internalVariables, self.dt) + for i in range(10): + Q = Rotation.random(random_state=i).as_matrix() + dispGradTransformed = Q@(dispGrad + np.identity(3)) - np.identity(3) + WStar = self.material.compute_energy_density(dispGradTransformed, self.internalVariables, self.dt) + self.assertAlmostEqual(W, WStar, 12) + + + def test_internal_state_update(self): + # generate a random displacement gradient + key = jax.random.PRNGKey(1) + dispGrad = jax.random.uniform(key, (3, 3)) + + internalVariablesNew = self.material.compute_state_new(dispGrad, self.internalVariables, self.dt) + # linear elastic has no internal variables + self.assertEqual(internalVariablesNew.size, 0) + + +if __name__ == "__main__": + unittest.main() \ No newline at end of file diff --git "a/optimism/test/material/test_LinearElastic.py\357\200\272Zone.Identifier" "b/optimism/test/material/test_LinearElastic.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git a/optimism/test/material/test_MaterialUniaxialSimulator.py b/optimism/test/material/test_MaterialUniaxialSimulator.py new file mode 100644 index 00000000..6717d690 --- /dev/null +++ b/optimism/test/material/test_MaterialUniaxialSimulator.py @@ -0,0 +1,34 @@ +import unittest + +import jax.numpy as np + +from optimism.material import MaterialUniaxialSimulator +from optimism.material import LinearElastic +from ..TestFixture import TestFixture + + +class MaterialUniaxialSimulatorFixture(TestFixture): + + def test_uniaxial_state_achieved(self): + E = 100.0 + nu = 0.25 + + properties = {"elastic modulus": E, + "poisson ratio": nu, + "strain measure": "logarithmic"} + material = LinearElastic.create_material_model_functions(properties) + engineering_strain_rate = 1e-3 + + def strain_history(t): + return engineering_strain_rate*t + + maxTime = 1000.0 + response = MaterialUniaxialSimulator.run(material, strain_history, maxTime, steps=20) + + for stress in response.stressHistory[1:]: + self.assertGreater(stress[0,0], 0.0) # axial stress is nonnegative + self.assertTrue(np.abs(np.all(stress.ravel()[1:]) < 1e-8)) # all other components near zero + + +if __name__ == '__main__': + unittest.main() diff --git "a/optimism/test/material/test_MaterialUniaxialSimulator.py\357\200\272Zone.Identifier" "b/optimism/test/material/test_MaterialUniaxialSimulator.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git a/optimism/test/material/test_RateSensitivity.py b/optimism/test/material/test_RateSensitivity.py new file mode 100644 index 00000000..b47aa395 --- /dev/null +++ b/optimism/test/material/test_RateSensitivity.py @@ -0,0 +1,93 @@ +import unittest + +import jax + +from optimism.material import Hardening +from ..TestFixture import TestFixture + +class RateSensitivityFixture(TestFixture): + def test_power_law_scales_correctly(self): + S = 10.0 + m = 2.0 + epsDot0 = 1e-3 + + compute_overstress = jax.grad(Hardening.power_law_rate_sensitivity) + + eqpsDot1 = 2.5 + dt = 1.5e-2 + eqpsOld = 0.37 + eqps1 = eqpsOld + eqpsDot1*dt + stress1 = compute_overstress(eqps1, eqpsOld, dt, S, m, epsDot0) + + eqpsDot2 = eqpsDot1*2.0 + eqps2 = eqpsOld + eqpsDot2*dt + stress2 = compute_overstress(eqps2, eqpsOld, dt, S, m, epsDot0) + + self.assertNear(stress2/stress1, (eqpsDot2/eqpsDot1)**(1/m), 10) + + + def test_property_parsing(self): + S = 10.0 + m = 1.0 + epsDot0 = 0.1 + Y0 = 3.0 + props = {"rate sensitivity": "power law", + "rate sensitivity exponent": m, + "rate sensitivity stress": S, + "reference plastic strain rate": epsDot0, + "hardening model": "linear", + "hardening modulus": 0.0, + "yield strength": Y0} + + hardening = Hardening.create_hardening_model(props) + + # viscous overstress will increase flow stress above initial yield + eqps = 2.0 + eqpsOld = 1.0 + dt = 1.0e1 + flowStress = hardening.compute_flow_stress(eqps, eqpsOld, dt) + self.assertGreater(flowStress, Y0) + + # if the time period is huge, the strain rate is so small that the + # overstress should be lost in the floating point truncation error + dt = 1e25 + flowStress = hardening.compute_flow_stress(eqps, eqpsOld, dt) + self.assertEqual(flowStress, Y0) + + +from optimism.material import J2Plastic + +class RateSentivityInsideJ2(TestFixture): + def test_kinetic_potential_works_inside_J2(self): + E = 1.0 + nu = 0.3 + + Y0 = 1e-4 + H = E / 50 + + S = Y0 + m = 2.0 + epsDot0 = 0.1 + + props = {"elastic modulus": E, + "poisson ratio": nu, + "yield strength": Y0, + "hardening model": "linear", + "hardening modulus": H, + "rate sensitivity": "power law", + "rate sensitivity stress": S, + "rate sensitivity exponent": m, + "reference plastic strain rate": epsDot0} + + material = J2Plastic.create_material_model_functions(props) + + key = jax.random.PRNGKey(0) + dispGrad = jax.random.uniform(key, (3, 3)) + internalState = material.compute_initial_state() + dt = 1.0 + W = material.compute_energy_density(dispGrad, internalState, dt) + self.assertGreater(W, 0.0) + + +if __name__ == "__main__": + unittest.main() \ No newline at end of file diff --git "a/optimism/test/material/test_RateSensitivity.py\357\200\272Zone.Identifier" "b/optimism/test/material/test_RateSensitivity.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git "a/optimism/test/patch.json\357\200\272Zone.Identifier" "b/optimism/test/patch.json\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git "a/optimism/test/patch_2_blocks.g\357\200\272Zone.Identifier" "b/optimism/test/patch_2_blocks.g\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git a/optimism/test/phasefield/.gitignore b/optimism/test/phasefield/.gitignore new file mode 100644 index 00000000..a1363379 --- /dev/null +++ b/optimism/test/phasefield/.gitignore @@ -0,0 +1 @@ +*.pdf diff --git "a/optimism/test/phasefield/.gitignore\357\200\272Zone.Identifier" "b/optimism/test/phasefield/.gitignore\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git a/optimism/test/phasefield/__init__.py b/optimism/test/phasefield/__init__.py new file mode 100644 index 00000000..e69de29b diff --git "a/optimism/test/phasefield/__init__.py\357\200\272Zone.Identifier" "b/optimism/test/phasefield/__init__.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git a/optimism/test/phasefield/plotPhaseFieldThresholdModel.py b/optimism/test/phasefield/plotPhaseFieldThresholdModel.py new file mode 100644 index 00000000..8a865e9d --- /dev/null +++ b/optimism/test/phasefield/plotPhaseFieldThresholdModel.py @@ -0,0 +1,54 @@ +from matplotlib import pyplot as plt +import unittest + +from optimism.JaxConfig import * +from optimism.phasefield import PhaseFieldThreshold as Model +from optimism import TensorMath +from optimism.phasefield.MaterialPointSimulator import MaterialPointSimulator + +plotting=True + +E = 10.0 +nu = 0.25 +l = 1.0 +maxStrain = 0.05 +critStrain = 0.025 +psiC = 0.5*E*critStrain**2 +G0 = psiC/(3.0/16/l) * 1.1 +props = {'elastic modulus': E, + 'poisson ratio': nu, + 'critical energy release rate': G0, + 'critical strain energy density': psiC, + 'regularization length': l, + 'kinematics': 'small deformations'} +strainRate = 1.0e-3 + + +class PhaseFieldThresholdUniaxialFixture(unittest.TestCase): + def setUp(self): + materialModel = Model.create_material_model_functions(props) + self.simulator = MaterialPointSimulator(materialModel, maxStrain, strainRate, steps=20) + + + def testUniaxial(self): + output = self.simulator.run() + + if plotting: + fig, axs = plt.subplots(2,2) + axs[0,0].plot(output.strainHistory, output.stressHistory, marker='o') + axs[0,0].set(xlabel='strain', ylabel='stress') + + axs[1,0].plot(output.strainHistory, output.phaseHistory, marker='o') + axs[1,0].set(xlabel='strain', ylabel='phase') + + axs[0,1].plot(output.strainHistory, output.energyHistory, marker='o') + axs[0,1].plot(output.strainHistory, psiC*np.ones_like(output.strainHistory), '--') + axs[0,1].set(xlabel='strain', ylabel='strain energy density') + + fig.set_size_inches(12.0, 10.0) + plt.tight_layout() + plt.savefig('phaseUniaxial.pdf') + + +if __name__ == '__main__': + unittest.main() diff --git "a/optimism/test/phasefield/plotPhaseFieldThresholdModel.py\357\200\272Zone.Identifier" "b/optimism/test/phasefield/plotPhaseFieldThresholdModel.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git a/optimism/test/phasefield/plotSandiaModelUniaxial.py b/optimism/test/phasefield/plotSandiaModelUniaxial.py new file mode 100644 index 00000000..142f10d6 --- /dev/null +++ b/optimism/test/phasefield/plotSandiaModelUniaxial.py @@ -0,0 +1,76 @@ +from matplotlib import pyplot as plt +import unittest + +from optimism.JaxConfig import * +from optimism.phasefield import SandiaModel as Model +from optimism import TensorMath +from optimism.phasefield.MaterialPointSimulator import MaterialPointSimulator + +plotting=True + +E = 10.0 +nu = 0.25 +l = 1.0 +H = 1e-1 * E +Y0 = 0.01*E +maxStrain = 0.05 +critStrain = 0.025 +psiC = 0.5*E*((H*critStrain + Y0)/(H + E))**2 +G0 = psiC/(3.0/16/l) * 1.1 +void0 = 0.0 +props = {'elastic modulus': E, + 'poisson ratio': nu, + 'yield strength': Y0, + 'hardening model': 'linear', + 'hardening modulus': H, + 'critical energy release rate': G0, + 'critical strain energy density': psiC, + 'regularization length': l, + 'kinematics': 'small deformations', + 'void growth prefactor': 1.0, + 'void growth exponent': 1.0, + 'initial void fraction': 0.0} +strainRate = 1.0e-3 + + +class PhaseFieldUniaxialFixture(unittest.TestCase): + def setUp(self): + materialModel = Model.create_material_model_functions(props) + self.simulator = MaterialPointSimulator(materialModel, maxStrain, strainRate, steps=20) + + + def testUniaxial(self): + output = self.simulator.run() + + eqpsHistory = vmap(lambda Q : Q[Model.STATE_EQPS])(output.internalVariableHistory) + voidHistory = vmap(lambda Q : Q[Model.STATE_VOID_FRACTION])(output.internalVariableHistory) + psiCHistory = Model.update_critical_energy(psiC, voidHistory) + + if plotting: + fig, axs = plt.subplots(2,3) + axs[0,0].plot(output.strainHistory, output.stressHistory, marker='o') + axs[0,0].set(xlabel='strain', ylabel='stress') + + axs[1,0].plot(output.strainHistory, output.phaseHistory, marker='o') + axs[1,0].set(xlabel='strain', ylabel='phase') + + axs[0,1].plot(output.strainHistory, eqpsHistory, marker='o') + axs[0,1].set(xlabel='strain', ylabel='equivalent plastic strain') + +# axs[1,1].plot(output.strainHistory, output.energyHistory, marker='o') +# axs[1,1].set(xlabel='strain', ylabel='Potential density') + axs[1,1].plot(output.strainHistory, voidHistory, marker='o') + axs[1,1].set(xlabel='strain', ylabel='Void fraction') + + axs[0,2].plot(output.strainHistory, psiCHistory, marker='o') + axs[0,2].plot(output.strainHistory, np.ones_like(psiCHistory)*psiC, linestyle='dashed') + axs[0,2].set(xlabel='strain', ylabel=r'$\Psi$') + axs[0,2].legend((r'$\Psi$',r'$\Psi_c$',r'$\Psi_{c0}$')) + + fig.set_size_inches(12.0, 10.0) + plt.tight_layout() + plt.savefig('phaseUniaxial.pdf') + + +if __name__ == '__main__': + unittest.main() diff --git "a/optimism/test/phasefield/plotSandiaModelUniaxial.py\357\200\272Zone.Identifier" "b/optimism/test/phasefield/plotSandiaModelUniaxial.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git a/optimism/test/phasefield/test_PhaseFieldLorentzPlastic.py b/optimism/test/phasefield/test_PhaseFieldLorentzPlastic.py new file mode 100644 index 00000000..4790f3db --- /dev/null +++ b/optimism/test/phasefield/test_PhaseFieldLorentzPlastic.py @@ -0,0 +1,115 @@ +from matplotlib import pyplot as plt +from jax import random +from scipy.spatial.transform import Rotation as R +import unittest + +from optimism.JaxConfig import * +from . import EquationSolver_Immersed_2 as EqSolver +from optimism import Objective +from .. import TestFixture +from ..MeshFixture import MeshFixture +from optimism.phasefield import PhaseFieldLorentzPlastic as Model +from optimism import SparseMatrixAssembler +from optimism import TensorMath +from optimism import Mesh + + +plotting=False + + +class GradOfPlasticPhaseFieldModelFixture(TestFixture.TestFixture): + + def setUp(self): + self.E = 100.0 + self.nu = 0.321 + self.Gc = 40.0 + self.psiC = 0.5*self.E + self.l = 1.0 + self.Y0 = 0.3*self.E + self.H = 1.0e-2*self.E + props = {'elastic modulus': self.E, + 'poisson ratio': self.nu, + 'critical energy release rate': self.Gc, + 'critical strain energy density': self.psiC, + 'regularization length': self.l, + 'yield strength': self.Y0, + 'hardening model': 'linear', + 'hardening modulus': self.H, + 'kinematics': 'large deformations'} + self.model = Model.create_material_model_functions(props) + self.flux_func = grad(self.model.compute_energy_density, (0,1,2)) + self.internalVariables = self.model.compute_initial_state() + self.dt = 1.0 + + + def test_zero_point(self): + dispGrad = np.zeros((3,3)) + phase = 0. + phaseGrad = np.zeros(3) + + energy = self.model.compute_energy_density(dispGrad, phase, phaseGrad, self.internalVariables, self.dt) + self.assertNear(energy, 0.0, 12) + + stress, phaseForce, phaseGradForce = self.flux_func(dispGrad, phase, phaseGrad, self.internalVariables, self.dt) + self.assertArrayNear(stress, np.zeros((3,3)), 12) + self.assertNear(phaseForce, 3.0/8.0*self.Gc/self.l, 12) + self.assertArrayNear(phaseGradForce, np.zeros(3), 12) + + + def test_rotation_invariance(self): + key = random.PRNGKey(0) + dispGrad = random.uniform(key, (3,3)) + key, subkey = random.split(key) + phase = random.uniform(subkey) + key,subkey = random.split(key) + phaseGrad = random.uniform(subkey, (3,)) + energy = self.model.compute_energy_density(dispGrad, phase, phaseGrad, self.internalVariables, self.dt) + + Q = R.random(random_state=1234).as_matrix() + dispGradStar = Q@(dispGrad + np.identity(3)) - np.identity(3) + phaseStar = phase + phaseGradStar = Q@phaseGrad + internalVariablesStar = self.internalVariables + energyStar = self.model.compute_energy_density(dispGradStar, phaseStar, phaseGradStar, internalVariablesStar, self.dt) + self.assertNear(energy, energyStar, 12) + + + def test_elastic_energy(self): + strainBelowYield = 0.5*self.Y0/self.E # engineering strain + dispGrad = np.diag(np.exp(strainBelowYield*np.array([1.0, -self.nu, -self.nu])))-np.identity(3) + phase = 0.0 + phaseGrad = np.zeros(3) + + energy = self.model.compute_energy_density(dispGrad, phase, phaseGrad, self.internalVariables, self.dt) + energyExact = 0.5*self.E*strainBelowYield**2 + self.assertNear(energy, energyExact, 12) + + piolaStress,_,_ = self.flux_func(dispGrad, phase, phaseGrad, self.internalVariables, self.dt) + mandelStress = piolaStress@(dispGrad + np.identity(3)).T + stressExact = np.zeros((3,3)).at[0,0].set(self.E*strainBelowYield) + self.assertArrayNear(mandelStress, stressExact, 12) + + + def test_plastic_stress(self): + + strain11 = 1.1*self.Y0/self.E + eqps = (self.E*strain11 - self.Y0)/(self.H + self.E) + elasticStrain11 = strain11 - eqps + lateralStrain = -self.nu*elasticStrain11 - 0.5*eqps + strains = np.array([strain11, lateralStrain, lateralStrain]) + dispGrad = np.diag(np.exp(strains)) - np.identity(3) + phase = 0.0 + phaseGrad = np.zeros(3) + + energyExact = 0.5*self.E*elasticStrain11**2 + self.Y0*eqps + 0.5*self.H*eqps**2 + energy = self.model.compute_energy_density(dispGrad, phase, phaseGrad, self.internalVariables, self.dt) + self.assertNear(energy, energyExact, 12) + + stress,_,_ = self.flux_func(dispGrad, phase, phaseGrad, self.internalVariables, self.dt) + mandelStress = stress@(dispGrad + np.identity(3)).T + mandelStress11Exact = self.E*(strain11 - eqps) + self.assertNear(mandelStress[0,0], mandelStress11Exact, 12) + + +if __name__ == '__main__': + unittest.main() diff --git "a/optimism/test/phasefield/test_PhaseFieldLorentzPlastic.py\357\200\272Zone.Identifier" "b/optimism/test/phasefield/test_PhaseFieldLorentzPlastic.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git a/optimism/test/phasefield/test_PhaseFieldLorentzPlasticPatch.py b/optimism/test/phasefield/test_PhaseFieldLorentzPlasticPatch.py new file mode 100644 index 00000000..f1ffb7dd --- /dev/null +++ b/optimism/test/phasefield/test_PhaseFieldLorentzPlasticPatch.py @@ -0,0 +1,240 @@ +from scipy.sparse import diags +from matplotlib import pyplot as plt + +from optimism.JaxConfig import * +from optimism.phasefield import PhaseFieldLorentzPlastic as material +from optimism.phasefield import PhaseField +from optimism.phasefield import MaterialPointSimulator +from .. import MeshFixture +from optimism import Mesh +from optimism import AlSolver +from . import EquationSolver_Immersed_2 as EqSolver +from optimism import FunctionSpace +from optimism import Interpolants +from optimism import Objective +from optimism.ConstrainedObjective import ConstrainedObjective +from optimism import QuadratureRule +from optimism import SparseMatrixAssembler +from optimism import TensorMath + + +E = 1.0 +nu = 0.25 +Y0 = 1.0 +H = 1e-2 * E +ell = 0.25 +rpOverL = 3.0 +Gc = 3.0*np.pi*Y0**2/(E/(1.0-nu**2)) * ell * rpOverL +psiC = 3/16 * Gc/ell + +props = {'elastic modulus': E, + 'poisson ratio': nu, + 'yield strength': Y0, + 'hardening model': 'linear', + 'hardening modulus': H, + 'critical energy release rate': Gc, + 'critical strain energy density': psiC, + 'regularization length': ell} + +materialModel = material.create_material_model_functions(props) + +subProblemSettings = EqSolver.get_settings() +alSettings = AlSolver.get_settings() + +class TestSingleMeshFixture(MeshFixture.MeshFixture): + + def setUp(self): + self.Nx = 7 + self.Ny = 7 + xRange = [0.,1.] + yRange = [0.,1.] + + self.targetFieldGrad = np.array([[0.1, -0.2],[0.4, -0.1], [0,0]]) + self.targetFieldOffset = np.array([0.5, -0.5, 0.0]) + + self.mesh, self.U = self.create_mesh_and_disp(self.Nx, self.Ny, xRange, yRange, + lambda x: self.targetFieldGrad@x + self.targetFieldOffset) + + quadRule = QuadratureRule.create_quadrature_rule_on_triangle(degree=2) + self.fs = FunctionSpace.construct_function_space(self.mesh, quadRule) + + dim = 3 # 2 displacements and scalar phase + ebc = [FunctionSpace.EssentialBC(nodeSet='all_boundary', component=0), + FunctionSpace.EssentialBC(nodeSet='all_boundary', component=1)] + self.dofManager = FunctionSpace.DofManager(self.fs, dim, ebc) + self.Ubc = self.dofManager.get_bc_values(self.U) + + self.bvpFunctions = PhaseField.create_phasefield_functions(self.fs, + "plane strain", + materialModel) + + self.nElements = Mesh.num_elements(self.mesh) + self.nQuadPtsPerElem = QuadratureRule.len(quadRule) + self.stateVars = self.bvpFunctions.compute_initial_state() + + dofToUnknown = self.dofManager.dofToUnknown.reshape(self.U.shape) + self.phaseIds = dofToUnknown[self.dofManager.isUnknown[:,2],2] + + + def test_sparse_hessian_at_zero_phase(self): + elementStiffnesses = self.bvpFunctions.\ + compute_element_stiffnesses(self.U, self.stateVars) + + KSparse = SparseMatrixAssembler.assemble_sparse_stiffness_matrix(elementStiffnesses, + self.mesh.conns, + self.dofManager) + KSparseDensified = np.array(KSparse.todense()) + + Ubc = self.dofManager.get_bc_values(self.U) + + def objective_func(Uu): + U = self.dofManager.create_field(Uu, Ubc) + return self.bvpFunctions.compute_internal_energy(U, self.stateVars) + + compute_dense_hessian = hessian(objective_func) + KDense = compute_dense_hessian(self.dofManager.get_unknown_values(self.U)) + + self.assertArrayNear(KSparseDensified, KDense, 11) + + + def test_sparse_hessian_at_nonzero_phase(self): + U = self.U.at[:,2].set(np.linspace(0.1, 0.9, self.Nx*self.Ny)) + + elementStiffnesses = self.bvpFunctions.\ + compute_element_stiffnesses(U, self.stateVars) + + KSparse = SparseMatrixAssembler.assemble_sparse_stiffness_matrix(elementStiffnesses, + self.mesh.conns, + self.dofManager) + KSparseDensified = np.array(KSparse.todense()) + + Ubc = self.dofManager.get_bc_values(U) + + def objective_func(Uu): + U = self.dofManager.create_field(Uu, Ubc) + return self.bvpFunctions.compute_internal_energy(U, self.stateVars) + + compute_dense_hessian = hessian(objective_func) + KDense = compute_dense_hessian(self.dofManager.get_unknown_values(U)) + + self.assertArrayNear(KSparseDensified, KDense, 11) + + + def test_constrained_hessian(self): + U = self.U + + c = U[self.dofManager.isUnknown[:,2],2] + lam = np.zeros_like(c) + kappa = np.ones_like(lam) + + elementStiffnesses =self.bvpFunctions.\ + compute_element_stiffnesses(U, self.stateVars) + + + KSparse = SparseMatrixAssembler.assemble_sparse_stiffness_matrix(elementStiffnesses, + self.mesh.conns, + self.dofManager) + constraintHessian = self.bvpFunctions.compute_constraint_hessian(lam, + kappa, + c, + self.dofManager) + + KSparse += diags(constraintHessian, format='csc') + KDensified = np.array(KSparse.todense()) + + # construct constrained obj + Ubc = self.dofManager.get_bc_values(U) + Uu = self.dofManager.get_unknown_values(U) + + def objective_function(Uu, p): + U = self.dofManager.create_field(Uu, Ubc) + return self.bvpFunctions.compute_internal_energy(U, self.stateVars) + + def constraint_function(Uu, p): + return Uu[self.phaseIds] + + objective = ConstrainedObjective(objective_function, + constraint_function, + Uu, + None, + lam, + kappa) + + KDense = objective.hessian(Uu) + + self.assertArrayNear(KDensified, KDense, 11) + + + def no_test_uniaxial(self): + """Restore this test to see the uniaxial stress response + """ + psiC = 3.0/16.0*Gc/ell + strainC = np.sqrt(2.0*psiC/E) + maxStrain = 3.0*strainC + strainRate = 1.0e-3 + simulator = MaterialPointSimulator.MaterialPointSimulator(compute_energy_density, + compute_energy_density, + compute_state_new, + maxStrain, + strainRate, + self.stateVars, + steps=20) + output = simulator.run() + + fig, axs = plt.subplots(2) + axs[0].plot(output.strainHistory, output.stressHistory, marker='o') + axs[0].set(xlabel='strain', ylabel='stress') + + axs[1].plot(output.strainHistory, output.phaseHistory, marker='o') + axs[1].set(xlabel='strain', ylabel='phase') + + fig.set_size_inches(6.0, 10.0) + plt.tight_layout() + plt.show() + #plt.savefig('phaseUniaxial.pdf') + + + def test_patch_test(self): + + def objective_function(Uu, p): + U = self.dofManager.create_field(Uu, self.Ubc) + return self.bvpFunctions.compute_internal_energy(U, self.stateVars) + + def constraint_function(Uu, p): + return Uu[self.phaseIds] + + Uu_guess = self.dofManager.get_unknown_values(self.U) + + p0 = Objective.Params(None, self.stateVars) + lam0 = grad(objective_function, 0)(Uu_guess, p0)[self.phaseIds] + kappaGuess = 10.0 #tbd + kappa0 = kappaGuess * np.ones(lam0.shape) + + objective = ConstrainedObjective(objective_function, + constraint_function, + Uu_guess, + p0, + lam0, + kappa0) + + + Uu = AlSolver.augmented_lagrange_solve(objective, Uu_guess, p0, + alSettings, subProblemSettings, + useWarmStart=False) + + U = self.dofManager.create_field(Uu, self.Ubc) + + fieldGrads = FunctionSpace.compute_field_gradient(self.fs, U) + fieldGrads = fieldGrads.reshape((self.nElements*self.nQuadPtsPerElem, + fieldGrads.shape[2], fieldGrads.shape[3])) + for fg in fieldGrads: + self.assertArrayNear(fg, self.targetFieldGrad, 10) + + lagrangian = lambda x : objective_function(x, 0.0) - constraint_function(x, 0.0) @ objective.lam + grad_lagrangian = grad(lagrangian) + self.assertArrayNear( grad_lagrangian(Uu), np.zeros(Uu.shape), 10) + + +if __name__ == '__main__': + MeshFixture.unittest.main() + diff --git "a/optimism/test/phasefield/test_PhaseFieldLorentzPlasticPatch.py\357\200\272Zone.Identifier" "b/optimism/test/phasefield/test_PhaseFieldLorentzPlasticPatch.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git a/optimism/test/phasefield/test_PhaseFieldThreshold.py b/optimism/test/phasefield/test_PhaseFieldThreshold.py new file mode 100644 index 00000000..d5414419 --- /dev/null +++ b/optimism/test/phasefield/test_PhaseFieldThreshold.py @@ -0,0 +1,90 @@ +from matplotlib import pyplot as plt +from jax import random +from scipy.spatial.transform import Rotation as R + +from optimism.JaxConfig import * +from . import EquationSolver_Immersed_2 as EqSolver +from optimism import Objective +from .. import TestFixture +from ..MeshFixture import MeshFixture +from optimism.phasefield import PhaseFieldThreshold as Model +from optimism import SparseMatrixAssembler +from optimism import TensorMath +from optimism import Mesh + + +plotting=False + + +class PhaseFieldThresholdModelFixture(TestFixture.TestFixture): + + def setUp(self): + self.E = 100.0 + self.nu = 0.321 + self.Gc = 40.0 + self.l = 1.0 + props = {'elastic modulus': self.E, + 'poisson ratio': self.nu, + 'critical energy release rate': self.Gc, + 'regularization length': self.l, + 'kinematics': 'large deformations'} + self.model = Model.create_material_model_functions(props) + self.flux_func = grad(self.model.compute_energy_density, (0,1,2)) + self.internalVariables = self.model.compute_initial_state() + self.dt = 0.0 # unused - material has no rate dependence + + + def test_zero_point(self): + dispGrad = np.zeros((3,3)) + phase = 0. + phaseGrad = np.zeros(3) + + energy = self.model.compute_energy_density(dispGrad, phase, phaseGrad, self.internalVariables, self.dt) + self.assertNear(energy, 0.0, 12) + + stress, phaseForce, phaseGradForce = self.flux_func(dispGrad, phase, phaseGrad, self.internalVariables, self.dt) + self.assertArrayNear(stress, np.zeros((3,3)), 12) + self.assertNear(phaseForce, 3.0/8.0*self.Gc/self.l, 12) + self.assertArrayNear(phaseGradForce, np.zeros(3), 12) + + + def test_rotation_invariance(self): + key = random.PRNGKey(0) + dispGrad = random.uniform(key, (3,3)) + key, subkey = random.split(key) + phase = random.uniform(subkey) + key,subkey = random.split(key) + phaseGrad = random.uniform(subkey, (3,)) + dt = 0.0 + energy = self.model.compute_energy_density(dispGrad, phase, phaseGrad, self.internalVariables, self.dt) + + Q = R.random(random_state=1234).as_matrix() + dispGradStar = Q@(dispGrad + np.identity(3)) - np.identity(3) + phaseStar = phase + phaseGradStar = Q@phaseGrad + internalVariablesStar = self.internalVariables + energyStar = self.model.compute_energy_density(dispGradStar, phaseStar, phaseGradStar, internalVariablesStar, self.dt) + self.assertNear(energy, energyStar, 12) + + + def test_uniaxial_energy(self): + strain = 0.1 # engineering strain + F = np.diag(np.exp(strain*np.array([1.0, -self.nu, -self.nu]))) + dispGrad = F - np.identity(3) + phase = 0.15 + phaseGrad = np.zeros(3) + + energy = self.model.compute_energy_density(dispGrad, phase, phaseGrad, self.internalVariables, self.dt) + + g = (1.0 - phase)**2 + energyExact = g*0.5*self.E*strain**2 + 3/8*self.Gc/self.l*phase + self.assertNear(energy, energyExact, 12) + + piolaStress,_,_ = self.flux_func(dispGrad, phase, phaseGrad, self.internalVariables, self.dt) + kStress = piolaStress@(dispGrad + np.identity(3)).T + kStressExact = np.zeros((3,3)).at[0,0].set(g*self.E*strain) + self.assertArrayNear(kStress, kStressExact, 12) + + +if __name__ == '__main__': + TestFixture.unittest.main() diff --git "a/optimism/test/phasefield/test_PhaseFieldThreshold.py\357\200\272Zone.Identifier" "b/optimism/test/phasefield/test_PhaseFieldThreshold.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git a/optimism/test/phasefield/test_PhaseFieldThresholdPatch.py b/optimism/test/phasefield/test_PhaseFieldThresholdPatch.py new file mode 100644 index 00000000..dd1376f6 --- /dev/null +++ b/optimism/test/phasefield/test_PhaseFieldThresholdPatch.py @@ -0,0 +1,202 @@ +from scipy.sparse import diags +from matplotlib import pyplot as plt + +from optimism.JaxConfig import * +from optimism.phasefield import PhaseFieldThreshold as material +from optimism.phasefield import PhaseField +from optimism.phasefield import MaterialPointSimulator +from .. import MeshFixture +from optimism import Mesh +from optimism import AlSolver +from . import EquationSolver_Immersed_2 as EqSolver +from optimism import FunctionSpace +from optimism import Interpolants +from optimism import Objective +from optimism.ConstrainedObjective import ConstrainedObjective +from optimism import QuadratureRule +from optimism import SparseMatrixAssembler +from optimism import TensorMath + + +E = 1.0 +nu = 0.25 +Gc = 0.2 +ell = 0.25 + +props = {'elastic modulus': E, + 'poisson ratio': nu, + 'critical energy release rate': Gc, + 'regularization length': ell, + 'kinematics': 'large deformations'} + +materialModel = material.create_material_model_functions(props) + +subProblemSettings = EqSolver.get_settings() +alSettings = AlSolver.get_settings() + +class TestSingleMeshFixture(MeshFixture.MeshFixture): + + def setUp(self): + self.Nx = 7 + self.Ny = 7 + xRange = [0.,1.] + yRange = [0.,1.] + + self.targetFieldGrad = np.array([[0.1, -0.2],[0.4, -0.1], [0,0]]) + self.targetFieldOffset = np.array([0.5, -0.5, 0.0]) + + self.mesh, self.U = self.create_mesh_and_disp(self.Nx, self.Ny, xRange, yRange, + lambda x: self.targetFieldGrad@x + self.targetFieldOffset) + + quadRule = QuadratureRule.create_quadrature_rule_on_triangle(degree=2) + self.fs = FunctionSpace.construct_function_space(self.mesh, quadRule) + + ebc = [FunctionSpace.EssentialBC(nodeSet='all_boundary', component=0), + FunctionSpace.EssentialBC(nodeSet='all_boundary', component=1)] + dim = 3 # 2 displacements and scalar phase + self.dofManager = FunctionSpace.DofManager(self.fs, dim, ebc) + self.Ubc = self.dofManager.get_bc_values(self.U) + + self.bvpFunctions = PhaseField.create_phasefield_functions(self.fs, + "plane strain", + materialModel) + + self.nElements = Mesh.num_elements(self.mesh) + self.nQuadPtsPerElem = QuadratureRule.len(quadRule) + self.stateVars = self.bvpFunctions.compute_initial_state() + + dofToUnknown = self.dofManager.dofToUnknown.reshape(self.U.shape) + self.phaseIds = dofToUnknown[self.dofManager.isUnknown[:,2],2] + + + def test_sparse_hessian_at_zero_phase(self): + elementStiffnesses = self.bvpFunctions.\ + compute_element_stiffnesses(self.U, self.stateVars) + + KSparse = SparseMatrixAssembler.assemble_sparse_stiffness_matrix(elementStiffnesses, + self.mesh.conns, + self.dofManager) + KSparseDensified = np.array(KSparse.todense()) + + Ubc = self.dofManager.get_bc_values(self.U) + + def objective_func(Uu): + U = self.dofManager.create_field(Uu, Ubc) + return self.bvpFunctions.compute_internal_energy(U, self.stateVars) + + compute_dense_hessian = hessian(objective_func) + KDense = compute_dense_hessian(self.dofManager.get_unknown_values(self.U)) + + self.assertArrayNear(KSparseDensified, KDense, 11) + + + def test_sparse_hessian_at_nonzero_phase(self): + U = self.U.at[:,2].set(np.linspace(0.1, 0.9, self.Nx*self.Ny)) + + elementStiffnesses = self.bvpFunctions.\ + compute_element_stiffnesses(U, self.stateVars) + + KSparse = SparseMatrixAssembler.assemble_sparse_stiffness_matrix(elementStiffnesses, + self.mesh.conns, + self.dofManager) + KSparseDensified = np.array(KSparse.todense()) + + Ubc = self.dofManager.get_bc_values(U) + + def objective_func(Uu): + U = self.dofManager.create_field(Uu, Ubc) + return self.bvpFunctions.compute_internal_energy(U, self.stateVars) + + compute_dense_hessian = hessian(objective_func) + KDense = compute_dense_hessian(self.dofManager.get_unknown_values(U)) + + self.assertArrayNear(KSparseDensified, KDense, 11) + + + def test_constrained_hessian(self): + U = self.U + + c = U[self.dofManager.isUnknown[:,2],2] + lam = np.zeros_like(c) + kappa = np.ones_like(lam) + + elementStiffnesses =self.bvpFunctions.\ + compute_element_stiffnesses(U, self.stateVars) + + KSparse = SparseMatrixAssembler.assemble_sparse_stiffness_matrix(elementStiffnesses, + self.mesh.conns, + self.dofManager) + constraintHessian = self.bvpFunctions.compute_constraint_hessian(lam, + kappa, + c, + self.dofManager) + + KSparse += diags(constraintHessian, 0, format='csc') + KDensified = np.array(KSparse.todense()) + + # construct constrained obj + Ubc = self.dofManager.get_bc_values(U) + Uu = self.dofManager.get_unknown_values(U) + + def objective_function(Uu, p): + U = self.dofManager.create_field(Uu, Ubc) + return self.bvpFunctions.compute_internal_energy(U, self.stateVars) + + def constraint_function(Uu, p): + U = self.dofManager.create_field(Uu, Ubc) + return U[self.dofManager.isUnknown[:,2],2] + + objective = ConstrainedObjective(objective_function, + constraint_function, + Uu, + None, + lam, + kappa) + + KDense = objective.hessian(Uu) + + self.assertArrayNear(KDensified, KDense, 11) + + + def test_patch_test(self): + + def objective_function(Uu, p): + U = self.dofManager.create_field(Uu, self.Ubc) + return self.bvpFunctions.compute_internal_energy(U, self.stateVars) + + def constraint_function(Uu, p): + return Uu[self.phaseIds] + + Uu_guess = self.dofManager.get_unknown_values(self.U) + + p0 = Objective.Params(None, self.stateVars) + lam0 = grad(objective_function, 0)(Uu_guess, p0)[self.phaseIds] + kappa0 = 10.0*np.amax(lam0) * np.ones(lam0.shape) + + objective = ConstrainedObjective(objective_function, + constraint_function, + Uu_guess, + p0, + lam0, + kappa0) + + Uu = AlSolver.augmented_lagrange_solve(objective, Uu_guess, p0, + alSettings, subProblemSettings, + useWarmStart=False) + + U = self.dofManager.create_field(Uu, self.Ubc) + + fieldGrads = FunctionSpace.compute_field_gradient(self.fs, U) + fieldGrads = fieldGrads.reshape((self.nElements*self.nQuadPtsPerElem, + fieldGrads.shape[2], fieldGrads.shape[3])) + for fg in fieldGrads: + self.assertArrayNear(fg, self.targetFieldGrad, 10) + + lagrangian = lambda x : objective_function(x, 0.0) - constraint_function(x, 0.0) @ objective.lam + grad_lagrangian = grad(lagrangian) + self.assertArrayNear( grad_lagrangian(Uu), np.zeros(Uu.shape), 10) + + +if __name__ == '__main__': + MeshFixture.unittest.main() + diff --git "a/optimism/test/phasefield/test_PhaseFieldThresholdPatch.py\357\200\272Zone.Identifier" "b/optimism/test/phasefield/test_PhaseFieldThresholdPatch.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git a/optimism/test/phasefield/test_PhaseFieldThresholdPlastic.py b/optimism/test/phasefield/test_PhaseFieldThresholdPlastic.py new file mode 100644 index 00000000..5c4f676d --- /dev/null +++ b/optimism/test/phasefield/test_PhaseFieldThresholdPlastic.py @@ -0,0 +1,186 @@ +import unittest +from matplotlib import pyplot as plt +from jax import random +from scipy.spatial.transform import Rotation as R + +from optimism.JaxConfig import * +from . import EquationSolver_Immersed_2 as EqSolver +from optimism import Objective +from ..TestFixture import TestFixture +from ..MeshFixture import MeshFixture +from optimism.phasefield import PhaseFieldThresholdPlastic as Model +from optimism import SparseMatrixAssembler +from optimism import TensorMath +from optimism import Mesh + + +plotting=True + +flux_func = jacfwd(Model.energy_density, (0,1,2)) + + +class GradOfPlasticPhaseFieldModelFixture(TestFixture): + + def setUp(self): + E = 100.0 + poisson = 0.321 + Gc = 40.0 + l = 1.0 + Y0 = 0.3*E + H = 1.0e-2*E + self.props = Model.make_properties(E, poisson, Gc, l, Y0, H) + + + def test_zero_point(self): + dispGrad = np.zeros((3,3)) + phase = 0. + phaseGrad = np.zeros(3) + state = Model.make_initial_state() + + energy = Model.energy_density(dispGrad, phase, phaseGrad, state, self.props) + self.assertNear(energy, 0.0, 12) + + stress, phaseForce, phaseGradForce = flux_func(dispGrad, phase, phaseGrad, state, self.props) + self.assertArrayNear(stress, np.zeros((3,3)), 12) + self.assertNear(phaseForce, 3.0/8.0*self.props.Gc/self.props.l, 12) + self.assertArrayNear(phaseGradForce, np.zeros(3), 12) + + + def test_rotation_invariance(self): + key = random.PRNGKey(0) + dispGrad = random.uniform(key, (3,3)) + key, subkey = random.split(key) + phase = random.uniform(subkey) + key,subkey = random.split(key) + phaseGrad = random.uniform(subkey, (3,)) + internalVariables = Model.make_initial_state() + energy = Model.energy_density(dispGrad, phase, phaseGrad, internalVariables, self.props) + + Q = R.random(random_state=1234).as_matrix() + dispGradStar = Q@dispGrad@Q.T + phaseStar = phase + phaseGradStar = Q@phaseGrad + internalVariablesStar = internalVariables + energyStar = Model.energy_density(dispGradStar, phaseStar, phaseGradStar, internalVariablesStar, self.props) + self.assertNear(energy, energyStar, 12) + + + def test_elastic_energy(self): + strainBelowYield = 0.5*self.props.Y0/self.props.E + dispGrad = strainBelowYield*np.diag(np.array([1.0, -self.props.nu, -self.props.nu])) + phase = 0.0 + phaseGrad = np.zeros(3) + state = Model.make_initial_state() + + energy = Model.energy_density(dispGrad, phase, phaseGrad, state, self.props) + WExact = 0.5*self.props.E*strainBelowYield**2 + self.assertNear(energy, WExact, 12) + + stress,_,_ = flux_func(dispGrad, phase, phaseGrad, state, self.props) + stressExact = np.zeros((3,3)).at[0,0].set(self.props.E*strainBelowYield) + self.assertArrayNear(stress, stressExact, 12) + + + def test_plastic_stress(self): + huge = 1e10 + + E = 100.0 + nu = 0.321 + Gc = huge + l = 1.0 + Y0 = 0.3*E + H = 1.0e-2*E + props = Model.make_properties(E, nu, Gc, l, Y0, H) + + strainAboveYield = 1.25*self.props.Y0/self.props.E + dispGrad = strainAboveYield*np.diag(np.array([1.0, -self.props.nu, -self.props.nu])) + phase = 0.0 + phaseGrad = np.zeros(3) + state = Model.make_initial_state() + stress,_,_ = flux_func(dispGrad, phase, phaseGrad, state, self.props) + print("stress=",stress) + + + def no_test_plastic_strain_path(self): + + strain = np.zeros((3,3)) + dispGrad = np.zeros((3,3)) + strain_inc = 0.1 + stateOld = np.zeros((10,)) + + energy_density = jit(J2.energy_density) + stress_func = jit(grad(lambda elStrain, state, props: J2.energy_density(elStrain, state, props, doUpdate=False))) + tangents_func = jit(hessian(J2.energy_density)) + compute_state_new = jit(J2.compute_state_new) + + strainHistory = [] + stressHistory = [] + tangentsHistory = [] + eqpsHistory = [] + energyHistory = [] + for i in range(10): + energy = energy_density(strain, stateOld, self.props) + tangentsNew = tangents_func(strain, stateOld, self.props) + stateNew = compute_state_new(strain, stateOld, self.props) + stressNew = stress_func(strain, stateNew, self.props) + strainHistory.append(strain[0,0]) + stressHistory.append(stressNew[0,0]) + tangentsHistory.append(tangentsNew[0,0,0,0]) + eqpsHistory.append(stateNew[J2.EQPS]) + energyHistory.append(energy) + + stateOld = stateNew + strain = strain.at[0,0].add(strain_inc) + dispGrad = strain + + if plotting: + plt.figure() + plt.plot(strainHistory, energyHistory, marker='o', fillstyle='none') + plt.xlabel('strain') + plt.ylabel('potential density') + + plt.savefig('energy.png') + + plt.figure() + plt.plot(strainHistory, stressHistory, marker='o', fillstyle='none') + plt.xlabel('strain') + plt.ylabel('stress') + + plt.savefig('stress.png') + + plt.figure() + plt.plot(strainHistory, eqpsHistory, marker='o') + plt.xlabel('strain') + plt.ylabel('Eq Pl Strain') + + plt.savefig('eqps.png') + + plt.figure() + plt.plot(strainHistory, tangentsHistory, marker='o') + plt.xlabel('strain') + plt.ylabel('Tangent modulus') + + plt.savefig('tangent.png') + + E = self.props[J2.PROPS_E] + mu = self.props[J2.PROPS_MU] + nu = self.props[J2.PROPS_NU] + lam = E*nu/(1.0 + nu)/(1.0 - 2.0*nu) + Y0 = self.props[J2.PROPS_Y0] + H = self.props[J2.PROPS_H] + strainEnd = 0.9 + # for solutions, refer to jax-fem/papers/plane_strain_unit_test.pdf + eqpsExp = [max(0.0, (2.0 * i * mu - Y0) / (3.0*mu + H) ) for i in strainHistory] + stressNewExp = (2.0 * mu * Y0 + \ + 2.0 * strainEnd * pow(mu,2) + \ + strainEnd * H * lam + \ + 2.0 * strainEnd * H * mu + \ + 3.0 * strainEnd * lam * mu) \ + / (3.0 * mu + H) + + self.assertNear( stressHistory[-1], stressNewExp, 11) + for i in range(0, len(strainHistory)): + self.assertNear( eqpsHistory[i], eqpsExp[i], 12) + +if __name__ == '__main__': + unittest.main() diff --git "a/optimism/test/phasefield/test_PhaseFieldThresholdPlastic.py\357\200\272Zone.Identifier" "b/optimism/test/phasefield/test_PhaseFieldThresholdPlastic.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git a/optimism/test/phasefield/test_PhaseFieldUniaxial.py b/optimism/test/phasefield/test_PhaseFieldUniaxial.py new file mode 100644 index 00000000..fe7e7d6f --- /dev/null +++ b/optimism/test/phasefield/test_PhaseFieldUniaxial.py @@ -0,0 +1,66 @@ +from matplotlib import pyplot as plt +import unittest +from optimism.JaxConfig import * +from optimism.phasefield import PhaseFieldThresholdPlastic +from optimism import TensorMath +from optimism.phasefield.MaterialPointSimulator import MaterialPointSimulator + +plotting=False + +E = 10.0 +nu = 0.25 +# Gc goes here in order +l = 1.0 +H = 1e-1 * E +Y0 = 0.01*E + +maxStrain = 0.05 +critStrain = 0.025 +psiC = 0.5*E*((H*critStrain + Y0)/(H + E))**2 +Gc = psiC/(3.0/16/l) + +props = PhaseFieldThresholdPlastic.make_properties(E, nu, Gc, l, Y0, H) + +def energy_density(dispGrad, phase, phaseGrad, internalVariables, doUpdate=True): + return PhaseFieldThresholdPlastic.energy_density(dispGrad, phase, phaseGrad, + internalVariables, props=props, doUpdate=doUpdate) + +def update(dispGrad, phase, phaseGrad, internalVariables): + return PhaseFieldThresholdPlastic.compute_state_new(dispGrad, phase, phaseGrad, + internalVariables, props=props) + + +class PhaseFieldUniaxialFixture(unittest.TestCase): + + @unittest.skip + def testUniaxial(self): + + internalVariables = PhaseFieldThresholdPlastic.make_initial_state() + + strainRate = 1.0e-3 + simulator = MaterialPointSimulator(energy_density, energy_density, update, maxStrain, strainRate, internalVariables, steps=20) + output = simulator.run() + + eqpsHistory = [Q[PhaseFieldThresholdPlastic.STATE_EQPS] for Q in output.internalVariableHistory] + + if plotting: + fig, axs = plt.subplots(2,2) + axs[0,0].plot(output.strainHistory, output.stressHistory, marker='o') + axs[0,0].set(xlabel='strain', ylabel='stress') + + axs[1,0].plot(output.strainHistory, output.phaseHistory, marker='o') + axs[1,0].set(xlabel='strain', ylabel='phase') + + axs[0,1].plot(output.strainHistory, eqpsHistory, marker='o') + axs[0,1].set(xlabel='strain', ylabel='equivalent plastic strain') + + axs[1,1].plot(output.strainHistory, output.energyHistory, marker='o') + axs[1,1].set(xlabel='strain', ylabel='Potential density') + + fig.set_size_inches(12.0, 10.0) + plt.tight_layout() + plt.savefig('phaseUniaxial.pdf') + + +if __name__ == '__main__': + unittest.main() diff --git "a/optimism/test/phasefield/test_PhaseFieldUniaxial.py\357\200\272Zone.Identifier" "b/optimism/test/phasefield/test_PhaseFieldUniaxial.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git "a/optimism/test/read_material_property_test.exo\357\200\272Zone.Identifier" "b/optimism/test/read_material_property_test.exo\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git a/optimism/test/test_AxisymmPatchTest.py b/optimism/test/test_AxisymmPatchTest.py index 4a2cbdfe..909bac1f 100644 --- a/optimism/test/test_AxisymmPatchTest.py +++ b/optimism/test/test_AxisymmPatchTest.py @@ -2,7 +2,7 @@ from optimism import FunctionSpace from optimism.material import LinearElastic as MatModel from optimism import Mechanics -from optimism.EquationSolver import newton_solve +from EquationSolver_Immersed_2 import newton_solve from optimism import QuadratureRule from . import MeshFixture diff --git "a/optimism/test/test_AxisymmPatchTest.py\357\200\272Zone.Identifier" "b/optimism/test/test_AxisymmPatchTest.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git "a/optimism/test/test_DofManager.py\357\200\272Zone.Identifier" "b/optimism/test/test_DofManager.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git a/optimism/test/test_EquationSolver.py b/optimism/test/test_EquationSolver.py index 8e66e1c0..41629fcc 100644 --- a/optimism/test/test_EquationSolver.py +++ b/optimism/test/test_EquationSolver.py @@ -1,5 +1,5 @@ from optimism.JaxConfig import * -from optimism import EquationSolver as es +from . import EquationSolver_Immersed_2 as es from optimism import Objective from .TestFixture import TestFixture import unittest diff --git "a/optimism/test/test_EquationSolver.py\357\200\272Zone.Identifier" "b/optimism/test/test_EquationSolver.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git "a/optimism/test/test_FunctionSpace.py\357\200\272Zone.Identifier" "b/optimism/test/test_FunctionSpace.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git "a/optimism/test/test_Interpolants.py\357\200\272Zone.Identifier" "b/optimism/test/test_Interpolants.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git "a/optimism/test/test_JaxConfig.py\357\200\272Zone.Identifier" "b/optimism/test/test_JaxConfig.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git "a/optimism/test/test_LinAlg.py\357\200\272Zone.Identifier" "b/optimism/test/test_LinAlg.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git a/optimism/test/test_Math.py b/optimism/test/test_Math.py index 04790807..4431c1d4 100644 --- a/optimism/test/test_Math.py +++ b/optimism/test/test_Math.py @@ -5,7 +5,7 @@ import jax.numpy as np from optimism import Math -from optimism.test import TestFixture +from . import TestFixture SUM_TESTDATA_FILENAME = os.path.join(os.path.dirname(__file__), 'ill_conditioned_sum_data.npz') DOTPROD_TESTDATA_FILENAME = os.path.join(os.path.dirname(__file__), 'ill_conditioned_dot_product_data.npz') diff --git "a/optimism/test/test_Math.py\357\200\272Zone.Identifier" "b/optimism/test/test_Math.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git "a/optimism/test/test_Mechanics.py\357\200\272Zone.Identifier" "b/optimism/test/test_Mechanics.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git "a/optimism/test/test_Mesh.py\357\200\272Zone.Identifier" "b/optimism/test/test_Mesh.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git "a/optimism/test/test_MinimizeScalar.py\357\200\272Zone.Identifier" "b/optimism/test/test_MinimizeScalar.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git a/optimism/test/test_Newmark.py b/optimism/test/test_Newmark.py index cd0b6b62..c69fb3dd 100644 --- a/optimism/test/test_Newmark.py +++ b/optimism/test/test_Newmark.py @@ -3,7 +3,7 @@ import jax import jax.numpy as np -from optimism import EquationSolver +from . import EquationSolver_Immersed_2 from optimism import FunctionSpace from optimism import Mechanics from optimism import Mesh @@ -18,7 +18,7 @@ nu = 0.0 rho = 1.0 -trSettings = EquationSolver.get_settings(max_cg_iters=50, +trSettings = EquationSolver_Immersed_2.get_settings(max_cg_iters=50, max_trust_iters=500, min_tr_size=1e-13, tol=4e-12, @@ -184,7 +184,7 @@ def time_step(self, Uu, Vu, Au, objective, dt): UuPredicted, Vu = self.dynamicsFunctions.predict(Uu, Vu, Au, dt) objective.p = Objective.param_index_update(objective.p, 5, UuPredicted) - Uu, solverSuccess = EquationSolver.nonlinear_equation_solve(objective, + Uu, solverSuccess = EquationSolver_Immersed_2.nonlinear_equation_solve(objective, UuPredicted, objective.p, trSettings, @@ -299,7 +299,7 @@ def energy(Uu, p): # The predictor value is exact. # Choose a bad initial guess (zero) to force an actual solve. - Uu, solverSuccess = EquationSolver.nonlinear_equation_solve(objective, + Uu, solverSuccess = EquationSolver_Immersed_2.nonlinear_equation_solve(objective, dofManager.get_unknown_values(UPrediction), objective.p, trSettings, @@ -365,7 +365,7 @@ def traction(X, N): # The predictor value is exact. # Choose a bad initial guess (zero) to force an actual solve. - Uu, solverSuccess = EquationSolver.nonlinear_equation_solve(objective, + Uu, solverSuccess = EquationSolver_Immersed_2.nonlinear_equation_solve(objective, dofManager.get_unknown_values(UPrediction), objective.p, trSettings, diff --git "a/optimism/test/test_Newmark.py\357\200\272Zone.Identifier" "b/optimism/test/test_Newmark.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git "a/optimism/test/test_Objective.py\357\200\272Zone.Identifier" "b/optimism/test/test_Objective.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git a/optimism/test/test_PatchTest.py b/optimism/test/test_PatchTest.py index 9bd797d7..0cd45ca5 100644 --- a/optimism/test/test_PatchTest.py +++ b/optimism/test/test_PatchTest.py @@ -9,7 +9,7 @@ from optimism import Mesh from optimism import Mechanics from optimism.Timer import Timer -from optimism.EquationSolver import newton_solve +from EquationSolver_Immersed_2 import newton_solve from optimism import QuadratureRule from . import MeshFixture diff --git "a/optimism/test/test_PatchTest.py\357\200\272Zone.Identifier" "b/optimism/test/test_PatchTest.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git a/optimism/test/test_PatchTestPou.py b/optimism/test/test_PatchTestPou.py index be7d3bcb..6e939214 100644 --- a/optimism/test/test_PatchTestPou.py +++ b/optimism/test/test_PatchTestPou.py @@ -7,7 +7,7 @@ #from optimism.material import Neohookean as MatModel from optimism import Mesh from optimism import Mechanics -from optimism.EquationSolver import newton_solve +from EquationSolver_Immersed_2 import newton_solve from optimism import QuadratureRule from . import MeshFixture from optimism.TensorMath import tensor_2D_to_3D diff --git "a/optimism/test/test_PatchTestPou.py\357\200\272Zone.Identifier" "b/optimism/test/test_PatchTestPou.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git "a/optimism/test/test_QuadratureRule.py\357\200\272Zone.Identifier" "b/optimism/test/test_QuadratureRule.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git a/optimism/test/test_ReadExodusMesh.py b/optimism/test/test_ReadExodusMesh.py index 6de0f6c0..ef5ed257 100644 --- a/optimism/test/test_ReadExodusMesh.py +++ b/optimism/test/test_ReadExodusMesh.py @@ -12,7 +12,7 @@ from optimism import FunctionSpace from optimism.material import LinearElastic as MaterialModel -from optimism.EquationSolver import newton_solve +from EquationSolver_Immersed_2 import newton_solve from optimism import QuadratureRule from . import TestFixture from optimism import Mechanics diff --git "a/optimism/test/test_ReadExodusMesh.py\357\200\272Zone.Identifier" "b/optimism/test/test_ReadExodusMesh.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git a/optimism/test/test_ReadMesh.py b/optimism/test/test_ReadMesh.py index 93130d60..878abdef 100644 --- a/optimism/test/test_ReadMesh.py +++ b/optimism/test/test_ReadMesh.py @@ -5,7 +5,7 @@ from optimism.material import LinearElastic as MaterialModel from optimism import Mesh from optimism.Timer import Timer -from optimism.EquationSolver import newton_solve +from EquationSolver_Immersed_2 import newton_solve from optimism import QuadratureRule from optimism import Surface from . import TestFixture diff --git "a/optimism/test/test_ReadMesh.py\357\200\272Zone.Identifier" "b/optimism/test/test_ReadMesh.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git "a/optimism/test/test_ScalarRootFinder.py\357\200\272Zone.Identifier" "b/optimism/test/test_ScalarRootFinder.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git "a/optimism/test/test_SmoothFunctions.py\357\200\272Zone.Identifier" "b/optimism/test/test_SmoothFunctions.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git "a/optimism/test/test_SparseMatrix.py\357\200\272Zone.Identifier" "b/optimism/test/test_SparseMatrix.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git "a/optimism/test/test_SparsePreconditioner.py\357\200\272Zone.Identifier" "b/optimism/test/test_SparsePreconditioner.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git "a/optimism/test/test_Surface.py\357\200\272Zone.Identifier" "b/optimism/test/test_Surface.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git "a/optimism/test/test_TensorMath.py\357\200\272Zone.Identifier" "b/optimism/test/test_TensorMath.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git a/optimism/test/test_Traction.py b/optimism/test/test_Traction.py index 75cd3eed..d07407e4 100644 --- a/optimism/test/test_Traction.py +++ b/optimism/test/test_Traction.py @@ -8,7 +8,7 @@ from optimism import Mesh from optimism import Mechanics from optimism.Timer import Timer -from optimism.EquationSolver import newton_solve +from EquationSolver_Immersed_2 import newton_solve from optimism import QuadratureRule from . import MeshFixture diff --git "a/optimism/test/test_Traction.py\357\200\272Zone.Identifier" "b/optimism/test/test_Traction.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git a/optimism/test/test_TrustRegionSPG.py b/optimism/test/test_TrustRegionSPG.py index 1385906a..7a776c56 100644 --- a/optimism/test/test_TrustRegionSPG.py +++ b/optimism/test/test_TrustRegionSPG.py @@ -8,7 +8,7 @@ from . import TestFixture from scipy import sparse -from optimism import EquationSolver +from . import EquationSolver_Immersed_2 def energy(x, params): p = params[0] @@ -272,8 +272,8 @@ def test_trust_region_spg_on_unbounded_problem(self): # self.assertTrue(onp.testing.assert_array_less(-eigvals, np.zeros(self.x.shape))) def no_test_cgunbound(self): - settings = EquationSolver.get_settings() - sol, solverSuccess = EquationSolver.nonlinear_equation_solve(self.obj, self.x, self.p, + settings = EquationSolver_Immersed_2.get_settings() + sol, solverSuccess = EquationSolver_Immersed_2.nonlinear_equation_solve(self.obj, self.x, self.p, settings, useWarmStart=False) self.assertTrue(solverSuccess) self.assertNear(np.linalg.norm(self.obj.gradient(sol)), 0, 7) @@ -314,8 +314,8 @@ def test_spg_on_rosenbrock(self): def no_test_steihaug_on_rosenbrock(self): - settings = EquationSolver.get_settings() - sol, solverSuccess = EquationSolver.nonlinear_equation_solve(self.obj, self.x, self.p, + settings = EquationSolver_Immersed_2.get_settings() + sol, solverSuccess = EquationSolver_Immersed_2.nonlinear_equation_solve(self.obj, self.x, self.p, settings, useWarmStart=False) self.assertTrue(solverSuccess) print('sol', sol) diff --git "a/optimism/test/test_TrustRegionSPG.py\357\200\272Zone.Identifier" "b/optimism/test/test_TrustRegionSPG.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git "a/optimism/test/test_VTKWriter.py\357\200\272Zone.Identifier" "b/optimism/test/test_VTKWriter.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git "a/optimism/test/test_VolumeAverageJ.py\357\200\272Zone.Identifier" "b/optimism/test/test_VolumeAverageJ.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git a/optimism/test/treigen/run_tests.sh b/optimism/test/treigen/run_tests.sh new file mode 100644 index 00000000..1bd55322 --- /dev/null +++ b/optimism/test/treigen/run_tests.sh @@ -0,0 +1 @@ +pytest -s diff --git "a/optimism/test/treigen/run_tests.sh\357\200\272Zone.Identifier" "b/optimism/test/treigen/run_tests.sh\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git a/optimism/test/treigen/test_treigen.py b/optimism/test/treigen/test_treigen.py new file mode 100644 index 00000000..70dbb9e7 --- /dev/null +++ b/optimism/test/treigen/test_treigen.py @@ -0,0 +1,81 @@ +from optimism.JaxConfig import * +from optimism.treigen import treigen +from testutils import * +from jax.numpy.linalg import norm + +import numpy as onp + + + +def test_tr_spd_and_inside(): + A = np.array( [np.array([1.0,2.0,0.3]), + np.array([2.0,4.5,0.0]), + np.array([0.3,0.0,5.0])] ) + b = np.array( onp.random.random(3) ) + + Delta = 100 + s = treigen.solve(A, b, Delta) + + assert_array_approx(A@s, -b) + assert( b@s < 0 ) + assert( norm(s) < Delta ) + + +def test_tr_spd_and_outside(): + A = np.array( [np.array([1.0,2.0,0.3]), + np.array([2.0,4.5,0.0]), + np.array([0.3,0.0,5.0])] ) + b = np.array( onp.random.random(3) ) + Delta = 1e-4 + s = treigen.solve(A, b, Delta) + assert( b@s < 0 ) + assert( norm(s) == approx(Delta) ) + + +def test_tr_not_spd_and_outside(): + A = np.array( [np.array([1.0,2.0,0.3]), + np.array([2.0,4.5,0.0]), + np.array([0.3,0.0,-5.0])] ) + b = np.array( onp.random.random(3) ) + Delta = 1e-4 + s = treigen.solve(A, b, Delta) + assert( b@s < 0 ) + assert( norm(s) == approx(Delta) ) + + +def test_tr_the_not_so_hard_case(): + A = np.array( [np.array([1.0, 2.0, 0.0]), + np.array([2.0, 4.5, 0.0]), + np.array([0.0, 0.0,-1.1])] ) + b = np.array([0.,0.0,1.]) + + Delta = 1e4 + s = treigen.solve(A, b, Delta) + assert( b@s < 0 ) + assert( norm(s) == approx(Delta) ) + + +def test_tr_the_easy_hard_case(): + A = np.array( [np.array([1.0, 2.0, 0.0]), + np.array([2.0, 4.5, 0.0]), + np.array([0.0, 0.0,-1.1])] ) + b = np.array([2.0,1.0,0.0]) + + Delta = 1e-2 + s = treigen.solve(A, b, Delta) + + assert( b@s < 0 ) + assert( norm(s) == approx(Delta) ) + + +def test_tr_the_hard_case(): + A = np.array( [np.array([1.0, 2.0, 0.0]), + np.array([2.0, 4.5, 0.0]), + np.array([0.0, 0.0,-1.1])] ) + b = np.array([2.0,1.0,0.0]) + + Delta = 1e3 + s = treigen.solve(A, b, Delta) + + assert( norm(s) == approx(Delta) ) + diff --git "a/optimism/test/treigen/test_treigen.py\357\200\272Zone.Identifier" "b/optimism/test/treigen/test_treigen.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git a/optimism/test/treigen/testutils.py b/optimism/test/treigen/testutils.py new file mode 100644 index 00000000..6f51bfe7 --- /dev/null +++ b/optimism/test/treigen/testutils.py @@ -0,0 +1,6 @@ +from pytest import approx + +def assert_array_approx(a,b): + assert(len(a)==len(b)) + for i in range(len(a)): + assert a[i] == approx(b[i]) diff --git "a/optimism/test/treigen/testutils.py\357\200\272Zone.Identifier" "b/optimism/test/treigen/testutils.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b diff --git "a/optimism/treigen/treigen.py\357\200\272Zone.Identifier" "b/optimism/treigen/treigen.py\357\200\272Zone.Identifier" new file mode 100644 index 00000000..e69de29b