diff --git a/Project.toml b/Project.toml index ce5721a..8325383 100644 --- a/Project.toml +++ b/Project.toml @@ -44,7 +44,7 @@ Distributions = "0.25" DocStringExtensions = "0.9" Enzyme = "0.13" Exodus = "0.14" -FiniteElementContainers = "0.10" +FiniteElementContainers = "0.11" ForwardDiff = "1" KernelAbstractions = "0.9" Krylov = "0.9" diff --git a/examples/cube/script_explicit_dynamics.jl b/examples/cube/script_explicit_dynamics.jl new file mode 100644 index 0000000..eaf3a2d --- /dev/null +++ b/examples/cube/script_explicit_dynamics.jl @@ -0,0 +1,48 @@ +using ConstitutiveModels +using Cthonios +using FiniteElementContainers + +# file management +mesh_file = Base.source_dir() * "/cube.g" +output_file = splitext(mesh_file)[1] * "-output.exo" + +# Times +times = TimeStepper(0., 1., 11) + +# Physics +physics = (; + cube = SolidMechanics( + ThreeDimensional(), NeoHookean() + ) +) +props = (; + cube = Dict{String, Any}( + "Young's modulus" => 1., + "Poisson's ratio" => 0.45 + ) +) + +# Boundary Conditions +func(_) = 1.0 + +vel_ics = InitialCondition[ + InitialCondition("displ_z", "cube", func) +] + +# Simulation +sim = SingleDomainSimulation( + mesh_file, output_file, + times, physics, props +) + +objective = Cthonios.ExplicitDynamicsObjective() +objective_cache, U, p = Cthonios.setup_caches(objective, sim, 0.1) + +# small hack we should remove +mesh = UnstructuredMesh(mesh_file) +vel_ics = InitialConditions(mesh, objective_cache.assembler.dof, vel_ics) +# small hack above we should removes + +Cthonios.initialize!(objective_cache, U, p; vel_ics = vel_ics) +solver = Cthonios.ExplicitSolver(objective_cache, p) +Cthonios.run!(solver, objective_cache, U, p, sim) diff --git a/examples/cube/script_quasistatic.jl b/examples/cube/script_quasistatic.jl index 5eb4b86..27c795c 100644 --- a/examples/cube/script_quasistatic.jl +++ b/examples/cube/script_quasistatic.jl @@ -26,10 +26,10 @@ func_1(x, t) = 1.0 * t func_2(x, t) = 0.0 dirichlet_bcs = [ - DirichletBC("displ_x", "ssx-", func_2), - DirichletBC("displ_y", "ssy-", func_2), - DirichletBC("displ_z", "ssz-", func_2), - DirichletBC("displ_z", "ssz+", func_1) + DirichletBC("displ_x", func_2; sideset_name = "ssx-"), + DirichletBC("displ_y", func_2; sideset_name = "ssy-"), + DirichletBC("displ_z", func_2; sideset_name = "ssz-"), + DirichletBC("displ_z", func_1; sideset_name = "ssz+") ] # Simulation @@ -38,7 +38,7 @@ sim = SingleDomainSimulation( times, physics, props; dirichlet_bcs=dirichlet_bcs ) -objective = Cthonios.QuasiStaticObjective() -objective_cache = Cthonios.setup_cache(objective, sim) -solver = Cthonios.TrustRegionSolver(objective_cache; use_warm_start=true) -Cthonios.run!(objective_cache, solver, sim) # eventually remove sim from call +objective = QuasiStaticObjective() +objective_cache, U, p = Cthonios.setup_caches(objective, sim) +solver = TrustRegionSolver(objective_cache, p; use_warm_start=true) +Cthonios.run!(solver, objective_cache, U, p, sim) # eventually remove sim from call diff --git a/examples/wip/clamped/script_v0.3_dynamics.jl b/examples/hole_array/script_implicit_dynamics.jl similarity index 76% rename from examples/wip/clamped/script_v0.3_dynamics.jl rename to examples/hole_array/script_implicit_dynamics.jl index 71bc830..e6e95c0 100644 --- a/examples/wip/clamped/script_v0.3_dynamics.jl +++ b/examples/hole_array/script_implicit_dynamics.jl @@ -3,22 +3,22 @@ using Cthonios using FiniteElementContainers # Mesh -mesh_file = Base.source_dir() * "/clamped.g" +mesh_file = Base.source_dir() * "/hole_array.exo" # Times -times = TimeStepper(0., 1.e-5, 11) +times = TimeStepper(0., 1., 100) # Physics physics = (; - clamped = SolidMechanics( + Block1 = SolidMechanics( PlaneStrain(), NeoHookean() # PlaneStrain(), LinearElastic() ) ) props = (; - clamped = Dict{String, Any}( + Block1 = Dict{String, Any}( "Young's modulus" => 1., - "Poisson's ratio" => 0.0 + "Poisson's ratio" => 0.3 ) ) props = Cthonios.create_properties(physics, props) @@ -26,12 +26,13 @@ props = Cthonios.create_properties(physics, props) # Boundary Conditions func_1(x, t) = -1. * t#(0., -5. * t) func_2(x, t) = 0.0 +func_3(x, t) = @SVector [0., -0.025 * t] dirichlet_bcs = DirichletBC[ DirichletBC("displ_x", "yminus_sideset", func_2), DirichletBC("displ_y", "yminus_sideset", func_2), DirichletBC("displ_x", "yplus_sideset", func_2), - DirichletBC("displ_y", "yplus_sideset", func_2) + DirichletBC("displ_y", "yplus_sideset", func_1) ] # Simulation setup @@ -40,7 +41,7 @@ sim = SingleDomainSimulation( mesh_file, times, physics, props; dirichlet_bcs=dirichlet_bcs ) -objective_cache = Cthonios.ImplicitDynamicsObjectiveCacheNew(sim) +objective_cache = Cthonios.ImplicitDynamicsObjectiveCache(sim) solver = Cthonios.NewtonSolver(objective_cache) # fill!(objective_cache.solution_rate, 10.) diff --git a/examples/hole_array/script_quasistatic.jl b/examples/hole_array/script_quasistatic.jl index 85e5f60..856418a 100644 --- a/examples/hole_array/script_quasistatic.jl +++ b/examples/hole_array/script_quasistatic.jl @@ -1,6 +1,8 @@ using ConstitutiveModels using Cthonios +using Enzyme using FiniteElementContainers +Enzyme.Compiler.VERBOSE_ERRORS[] = true # function sim_test() # file management @@ -8,30 +10,50 @@ mesh_file = Base.source_dir() * "/mesh/hole_array.exo" output_file = splitext(mesh_file)[1] * "-output.exo" # Times -times = TimeStepper(0., 1., 40) +# times = TimeStepper(0., 1., 4) +times = TimeStepper(0., 1., 20) # Physics physics = (; Block1 = SolidMechanics( PlaneStrain(), NeoHookean() + # PlaneStrain(), LinearElastic() + # PlaneStrain(), Gent() + # PlaneStrain(), LinearElastoPlasticity(VonMisesYieldSurface(LinearIsotropicHardening())) ) ) +E = 70.e3 +ν = 0.3 +σ_y = 200. +H = 180. props = (; Block1 = Dict{String, Any}( + "density" => 1., "Young's modulus" => 1., - "Poisson's ratio" => 0.45 + "Poisson's ratio" => 0.48, + # "Poisson's ratio" => 0.3, + "Jm" => 3. ) + # Block1 = Dict( + # "Young's modulus" => E, + # "Poisson's ratio" => ν, + # "yield surface" => "VonMisesYieldSurface", + # "isotropic hardening model" => "LinearIsotropicHardening", + # "yield stress" => σ_y, + # "hardening modulus" => H + # ) ) # Boundary Conditions +# func_1(x, t) = -0.075 * t func_1(x, t) = -7.5 * t func_2(x, t) = 0.0 dirichlet_bcs = [ - DirichletBC("displ_x", "yminus_sideset", func_2), - DirichletBC("displ_y", "yminus_sideset", func_2), - DirichletBC("displ_x", "yplus_sideset", func_2), - DirichletBC("displ_y", "yplus_sideset", func_1) + DirichletBC("displ_x", func_2; sideset_name = "yminus_sideset"), + DirichletBC("displ_y", func_2; sideset_name = "yminus_sideset"), + DirichletBC("displ_x", func_2; sideset_name = "yplus_sideset"), + DirichletBC("displ_y", func_1; sideset_name = "yplus_sideset") ] # Simulation @@ -43,48 +65,57 @@ sim = SingleDomainSimulation( objective = QuasiStaticObjective() objective_cache, U, p = setup_caches(objective, sim) -qoi = QOIExtractor( +qoi1 = QOIExtractor( objective_cache, helmholtz_free_energy, +, FiniteElementContainers.L2QuadratureField, Float64; + is_material_qoi = true, reduction_2 = + ) +qoi2 = QOIExtractor( + objective_cache, residual, identity, + H1Field, Float64; + is_field_qoi = true, + reduction_2 = identity +) solver = TrustRegionSolver(objective_cache, p; use_warm_start=true) -Cthonios.run!(solver, objective_cache, U, p, sim) # eventually remove sim from call - -function design_objective!(f, qois, Us, ps, params) - U, p = Us[end], ps[end] - # example for setting coords as params - copyto!(p.h1_coords, params) +# solver = Cthonios.NewtonSolver(objective_cache) - # example for setting properties as params - # for (n, val) in enumerate(values(params)) - # copyto!(values(p.properties)[n], val) - # end - Cthonios._value!(f, qois[1], U, p) - return nothing -end +Cthonios.run!(solver, objective_cache, U, p, sim) -obj = Cthonios.DesignObjective(design_objective!, [qoi], U, p) +# function design_objective!(f, qoi_storages, qois, Us, ps, params) +# U, p = Us[end], ps[end] -coord_params = deepcopy(p.h1_coords) -props_params = deepcopy(p.properties) -display(props_params) -push!(obj.stored_solutions, deepcopy(U)) -push!(obj.stored_parameters, deepcopy(p)) +# # example for setting coords as params +# copyto!(p.h1_coords, params) +# # example for setting properties as params +# # for (n, val) in enumerate(values(params)) +# # copyto!(values(p.properties)[n], val) +# # end +# Cthonios.value!(f, qoi_storages[end], qois[1], U, p) +# return nothing +# end +# obj = Cthonios.DesignObjective(design_objective!, [qoi1, qoi2], U, p) +# Cthonios.forward_problem!(obj, solver, objective_cache, U, p) -# obj = Cthonios.DesignObjective(design_objective!, sens) +# # pick out parameters for design objective function +# coord_params = deepcopy(p.h1_coords) +# # props_params = deepcopy(p.properties) -f = Cthonios.value(obj, coord_params) -f, dfdX = Cthonios.gradient_and_value(obj, coord_params) -# f, dfdp = Cthonios.gradient_and_value(obj, props_params) +# f = Cthonios.value(obj, coord_params) +# # f = zeros(1) +# f = create_field(objective_cache) +# Cthonios.value!(f, qoi2, U, p) +# f +# # f, dfdX = Cthonios.gradient_and_value(obj, coord_params) +# # # f, dfdp = Cthonios.gradient_and_value(obj, props_params) -display(f) -display(dfdX) -# # display(dfdU) -# # end -# # sim_test() \ No newline at end of file +# # display(f) +# # display(dfdX) +# # # # display(dfdU) +# # # # end +# # # # sim_test() \ No newline at end of file diff --git a/examples/hole_array/script_quasistatic_to_eigen.jl b/examples/hole_array/script_quasistatic_to_eigen.jl new file mode 100644 index 0000000..269471f --- /dev/null +++ b/examples/hole_array/script_quasistatic_to_eigen.jl @@ -0,0 +1,49 @@ +using ConstitutiveModels +using Cthonios +using FiniteElementContainers + +# file management +# mesh_file = Base.source_dir() * "/hole_array_tri6.exo" +mesh_file = Base.source_dir() * "/mesh/hole_array.exo" +output_file = splitext(mesh_file)[1] * "-output.exo" + +# Times +times = TimeStepper(0., 1., 40) + +# Physics +physics = (; + Block1 = SolidMechanics( + PlaneStrain(), NeoHookean() + ) +) +props = (; + Block1 = Dict{String, Any}( + "Young's modulus" => 1., + "Poisson's ratio" => 0.45 + ) +) + +# Boundary Conditions +func_1(x, t) = -0.0 * t +func_2(x, t) = 0.0 + +dirichlet_bcs = [ + DirichletBC("displ_x", "yminus_sideset", func_2), + DirichletBC("displ_y", "yminus_sideset", func_2), + DirichletBC("displ_x", "yplus_sideset", func_2), + DirichletBC("displ_y", "yplus_sideset", func_1) +] + +# Simulation +sim = SingleDomainSimulation( + mesh_file, output_file, + times, physics, props; + dirichlet_bcs=dirichlet_bcs +) +solver = x -> Cthonios.TrustRegionSolverGPU(x; use_warm_start=true) +objective_cache = Cthonios.run!(sim, QuasiStaticObjectiveCache, solver) + +solver = Cthonios.EigenSolver(objective_cache, 10) +Cthonios.solve!(solver) + +# post-processing \ No newline at end of file diff --git a/examples/wip/block_on_block_2d/script.jl b/examples/wip/block_on_block_2d/script.jl index 4eb69ee..7299df2 100644 --- a/examples/wip/block_on_block_2d/script.jl +++ b/examples/wip/block_on_block_2d/script.jl @@ -3,8 +3,9 @@ using Cthonios using FiniteElementContainers mesh_file = Base.source_dir() * "/mesh_tri3.exo" +output_file = splitext(mesh_file)[1] * "-output.exo" -times = TimeStepper(0., 1., 20) +times = TimeStepper(0., 1., 400) physics = (; bottom = SolidMechanics(PlaneStrain(), NeoHookean()), @@ -16,20 +17,20 @@ props = (; "Poisson's ratio" => 0.3 ), top = Dict{String, Any}( - "Young's modulus" => 0.25, + "Young's modulus" => 10., "Poisson's ratio" => 0.3 ) ) -props = Cthonios.create_properties(physics, props) +# props = Cthonios.create_properties(physics, props) func_1(x, t) = 0. func_2(x, t) = -0.5 * t dirichlet_bcs = [ - DirichletBC("displ_x", "bottom", func_1), - DirichletBC("displ_y", "bottom", func_1), - DirichletBC("displ_x", "top", func_1), - DirichletBC("displ_y", "top", func_2) + DirichletBC("displ_x", func_1; sideset_name = "bottom"), + DirichletBC("displ_y", func_1; sideset_name = "bottom"), + DirichletBC("displ_x", func_1; sideset_name = "top"), + DirichletBC("displ_y", func_2; sideset_name = "top") ] contact_pairs = [ @@ -37,48 +38,81 @@ contact_pairs = [ ] sim = SingleDomainSimulation( - mesh_file, times, physics, props; + mesh_file, output_file, + times, physics, props; dirichlet_bcs=dirichlet_bcs ) -objective_cache = QuasiStaticObjectiveCache(sim) -# solver = TrustRegionSolver(objective_cache, Cthonios.parameters(objective_cache), TimerOutput()) -# mesh = UnstructuredMesh(mesh_file) -# pp = PostProcessor(mesh, "output.e", objective_cache.assembler.dof.var) -# Cthonios.run!(solver, pp) -# run!(solver) -# +objective = QuasiStaticObjective() +objective_cache, U, p = setup_caches(objective, sim) +# solver = TrustRegionSolver(objective_cache, p) +# Cthonios.run!(solver, objective_cache, U, p, sim) + mesh = UnstructuredMesh(mesh_file) -X = mesh.nodal_coords -contact_enforcement = Cthonios.MortarContact(0.05, 0.1) -contact_formulation = Cthonios.PenaltyContact() -contact_caches = Cthonios.create_contact_pair_caches(mesh, contact_pairs) -# @time Cthonios.update_interactions!(contact_caches, X, U) -# @time Cthonios.update_closest_edges_and_weights!(contact_caches, X, U) -FiniteElementContainers.update_time!(objective_cache.parameters) -FiniteElementContainers.update_bc_values!(objective_cache.parameters) -U = objective_cache.solution - -block_top_conns = mesh.element_conns.top.data |> unique |> sort -U[:, block_top_conns] .= -0.05 -@time Cthonios.update_nearest_neighbors!(contact_caches, X, U) -# @time Cthonios.contact_potential(contact, contact_caches, X, U, 0.0005) -@time @show Cthonios.assemble_contact_scalar!( - contact_enforcement, contact_formulation, contact_caches, X, U +dof = objective_cache.assembler.dof +# X = mesh.nodal_coords +X = p.h1_coords +contact_caches = Cthonios.create_contact_pair_caches( + mesh, dof, contact_pairs; + max_neighbors = 4 ) -@time Cthonios.assemble_contact_vector!( - contact_enforcement, contact_formulation, contact_caches, X, U + +mortar = Cthonios.MortarContact(0.05, 0.1) +contact_cache = Cthonios.PenaltyContactObjectiveCache( + mortar, contact_caches, objective_cache ) -# @time Cthonios.assemble_contact_matrix!( -# contact_enforcement, contact_formulation, contact_caches, X, U -# ) +solver = TrustRegionSolver( + contact_cache, p; + max_trust_iters = 500, + preconditioner = Cthonios.NoPreconditioner, + use_warm_start = false +) +# solver = Cthonios.NewtonSolver(contact_cache) +# Cthonios.update_nearest_neighbors!(contact_caches, X, U) +# Cthonios.step!(solver, contact_cache, U, p) +# Cthonios.step!(solver, contact_cache, U, p) + +# Cthonios.value(contact_cache, U, p) +# Cthonios.gradient(contact_cache, U, p) +Cthonios.run!(solver, contact_cache, U, p, sim) + +# # # let's move the top block a little bit to make it penetrate +# # top_nodes = unique(getfield(mesh.element_conns, :top)) +# # # @show top_nodes +# # X[2, top_nodes] .-= 0.01 -# @show "Edge 1" -# edge_1_close = values(contact_caches)[1].interactions[:, end] -# @show X[:, values(contact_caches[1]).driver_surf.side_nodes[:, end]] +# FiniteElementContainers.update_time!(p) +# FiniteElementContainers.update_bc_values!(p) +# FiniteElementContainers.update_field_dirichlet_bcs!(U, p.dirichlet_bcs) -# for edge in edge_1_close -# @show edge -# @show X[:, values(contact_caches[1]).follow_surf.side_nodes[:, edge]] -# @show Cthonios._compute_normal(X[:, values(contact_caches[1]).follow_surf.side_nodes[:, edge]]) +# # Cthonios._compute_closest_distance_to_each_side( +# # contact_caches, X, U +# # ) +# function constraint_func(contact_caches, U, p) +# Cthonios._compute_closest_distance_to_each_side(contact_caches, p.h1_coords, U) # end + +# c = constraint_func(contact_caches, U, p) +# κ = 4. * ones(length(c)) +# λ = 1e-4 * abs.(κ .* c) + +# # constraint_func(nothing, U, p) +# mortar = Cthonios.MortarContact(0.05, 0.1) +# penalty = Cthonios.PenaltyContact() +# al_cache = Cthonios.ConstrainedObjectiveCache( +# objective_cache, constraint_func, +# contact_caches +# ) +# # Cthonios.gradient(al_objective_cache, U, p, λ0, κ0) +# solver = Cthonios.AugmentedLagrangeSolver(al_cache, p) +# # Cthonios.value(al_cache, U, p, λ, κ) +# # Cthonios.solve!(solver, U, p, λ, κ) + +# cenergy = zeros(1) +# cvector = create_field(al_cache.objective_cache) +# Cthonios.assemble_contact_scalar!( +# cenergy, mortar, penalty, contact_caches, X, U +# ) +# Cthonios.assemble_contact_vector!( +# cvector, mortar, penalty, contact_caches, X, U +# ) diff --git a/examples/wip/clamped/script_explicit_dynamics.jl b/examples/wip/clamped/script_explicit_dynamics.jl new file mode 100644 index 0000000..865210f --- /dev/null +++ b/examples/wip/clamped/script_explicit_dynamics.jl @@ -0,0 +1,67 @@ +using ConstitutiveModels +using Cthonios +using FiniteElementContainers + +# Mesh +mesh_file = Base.source_dir() * "/clamped.g" +output_file = splitext(mesh_file)[1] * "-output.exo" + +# Times +times = TimeStepper(0., 1.e-5, 101) + +# Physics +physics = (; + clamped = SolidMechanics( + ThreeDimensional(), LinearElastic() + ) +) +props = (; + clamped = Dict{String, Any}( + "density" => 1000.0, + "Young's modulus" => 1.0e+9, + "Poisson's ratio" => 0.0 + ) +) + +# Initial conditions +a = 0.01 +s = 0.02 +func_1(x) = a * exp(-x[3] * x[3] / s / s / 2) + +displ_ics = InitialCondition[ + InitialCondition("displ_z", "clamped", func_1) +] + +# Boundary Conditions +func_2(x, t) = 0.0 + +dirichlet_bcs = DirichletBC[ + DirichletBC("displ_x", func_2; nodeset_name = "nsx-") + DirichletBC("displ_x", func_2; nodeset_name = "nsx+") + DirichletBC("displ_y", func_2; nodeset_name = "nsy-") + DirichletBC("displ_y", func_2; nodeset_name = "nsy+") + DirichletBC("displ_z", func_2; nodeset_name = "nsz-") + DirichletBC("displ_z", func_2; nodeset_name = "nsz+") + +] + +# Simulation setup +sim = SingleDomainSimulation( + mesh_file, output_file, + times, physics, props; + dirichlet_bcs=dirichlet_bcs +) +objective = Cthonios.ExplicitDynamicsObjective() +objective_cache, U, p = Cthonios.setup_caches(objective, sim, 0.2) + +# small hack we should remove +mesh = UnstructuredMesh(mesh_file) +displ_ics = InitialConditions(mesh, objective_cache.assembler.dof, displ_ics) +# small hack above we should removes + +Cthonios.initialize!(objective_cache, U, p; displ_ics = displ_ics) +solver = Cthonios.ExplicitSolver(objective_cache, p) +Cthonios.run!( + solver, objective_cache, U, p, sim; + output_exodus_every = 1 +) diff --git a/src/Cthonios.jl b/src/Cthonios.jl index 97fc036..44ea616 100644 --- a/src/Cthonios.jl +++ b/src/Cthonios.jl @@ -28,6 +28,7 @@ RuntimeGeneratedFunctions.init(@__MODULE__) # Re-exports export DirichletBC +export InitialCondition export PlaneStrain export ThreeDimensional export TimerOutput diff --git a/src/PostProcessor.jl b/src/PostProcessor.jl index de00ead..3c88add 100644 --- a/src/PostProcessor.jl +++ b/src/PostProcessor.jl @@ -11,7 +11,8 @@ end function PostProcessor(objective_cache, U, p, sim) mesh = UnstructuredMesh(sim.mesh_file) V_q = FunctionSpace(mesh, L2QuadratureField, Lagrange) - disp_var = objective_cache.assembler.dof.var + # disp_var = objective_cache.assembler.dof.var + disp_var = assembler(objective_cache).dof.var mat_outputs, mat_vars = create_material_output(objective_cache, V_q, StandardMaterialOutput{Float64}) update_material_output!(mat_outputs, objective_cache, U, p) diff --git a/src/cli/Parser.jl b/src/cli/Parser.jl index d6adcf7..8991cad 100644 --- a/src/cli/Parser.jl +++ b/src/cli/Parser.jl @@ -20,7 +20,7 @@ function _parse_dirichclet_bcs(settings) func = @RuntimeGeneratedFunction(Meta.parse(bc[:function])) for field in bc[:fields] for sset in bc[:sidesets] - push!(bcs, DirichletBC(field, sset, func)) + push!(bcs, DirichletBC(field, func; sideset_name = sset)) end end end diff --git a/src/contact/AD.jl b/src/contact/AD.jl new file mode 100644 index 0000000..a2831e5 --- /dev/null +++ b/src/contact/AD.jl @@ -0,0 +1,74 @@ +function contact_residual_side_a( + type::AbstractContactFormulation, + X_a, X_b, U_a, U_b, + max_overlap_dist, + rel_smoothing_size +) + return ForwardDiff.gradient( + z -> contact_energy( + type, + X_a, X_b, z, U_b, + max_overlap_dist, rel_smoothing_size + ), U_a + ) +end + +function contact_residual_side_b( + type::AbstractContactFormulation, + X_a, X_b, U_a, U_b, + max_overlap_dist, + rel_smoothing_size +) + return ForwardDiff.gradient( + z -> contact_energy( + type, + X_a, X_b, U_a, z, + max_overlap_dist, rel_smoothing_size + ), U_b + ) +end + +function contact_stiffness_side_a_side_a( + type::AbstractContactFormulation, + X_a, X_b, U_a, U_b, + max_overlap_dist, + rel_smoothing_size +) + return ForwardDiff.hessian( + z -> contact_energy( + type, + X_a, X_b, z, U_b, + max_overlap_dist, rel_smoothing_size + ), U_a + ) +end + +function contact_stiffness_side_b_side_b( + type::AbstractContactFormulation, + X_a, X_b, U_a, U_b, + max_overlap_dist, + rel_smoothing_size +) + return ForwardDiff.hessian( + z -> contact_energy( + type, + X_a, X_b, U_a, z, + max_overlap_dist, rel_smoothing_size + ), U_b + ) +end + +function contact_stiffness_side_a_side_b( + type::AbstractContactFormulation, + X_a, X_b, U_a, U_b, + max_overlap_dist, + rel_smoothing_size +) + return ForwardDiff.jacobian( + z -> contact_residual_side_a( + type, + X_a, X_b, U_a, z, + max_overlap_dist, rel_smoothing_size + ), U_b + ) +end diff --git a/src/contact/ClosestPointProjection.jl b/src/contact/ClosestPointProjection.jl index 3ac10f6..16cfc87 100644 --- a/src/contact/ClosestPointProjection.jl +++ b/src/contact/ClosestPointProjection.jl @@ -45,4 +45,32 @@ function _cpp_distance(edge, p) end return d -end \ No newline at end of file +end + +function _compute_closest_distance_to_each_side( + contact_caches, X, U +) + # loop over contact pairs + cpps = Float64[] + for cpair in values(contact_caches) + # loop over interactions in current contact pair + for n in axes(cpair.interactions, 2) + side = cpair.side_b.sides[n] + X_I = facet_coordinates(cpair.side_b, X, U, n) + # loop over quadrature points + cpp = 1e64 # some really big number + for q in 1:num_quadrature_points(cpair.side_b.ref_fe) + X_Q = facet_quadrature_point_coordinates(cpair.side_b, X_I, q, side) + # loop over neighbors in interaction list + for m in axes(cpair.interactions, 1) + X_M = facet_coordinates(cpair.side_a, X, U, cpair.interactions[m, n]) + cpp = min(cpp, _cpp_distance(X_M, X_Q)) + end + # display(X_Q) + end + # @show cpp + push!(cpps, cpp) + end + end + cpps +end diff --git a/src/contact/Contact.jl b/src/contact/Contact.jl index 06b3cdf..dcd4c64 100644 --- a/src/contact/Contact.jl +++ b/src/contact/Contact.jl @@ -5,18 +5,6 @@ abstract type AbstractContactEnforcement end abstract type AbstractContactFormulation end -# only for 2D -# assumes.... -# closest point projection distance -include("ContactSurface.jl") - -include("ContactPair.jl") - -include("PenaltyContact.jl") - -include("Mortar.jl") -include("Search.jl") - function assemble_contact_matrix!( enforcement::AbstractContactEnforcement, formulation::AbstractContactFormulation, @@ -29,107 +17,54 @@ function assemble_contact_matrix!( end function assemble_contact_scalar!( + energy, enforcement::AbstractContactEnforcement, formulation::AbstractContactFormulation, caches::NamedTuple, + # caches, X, U ) for cache in values(caches) - assemble_contact_scalar!(enforcement, formulation, cache, X, U) + assemble_contact_scalar!(energy, enforcement, formulation, cache, X, U) end end function assemble_contact_vector!( + vector, enforcement::AbstractContactEnforcement, formulation::AbstractContactFormulation, caches::NamedTuple, X, U ) for cache in values(caches) - assemble_contact_vector!(enforcement, formulation, cache, X, U) + assemble_contact_vector!(vector, enforcement, formulation, cache, X, U) end end -function contact_residual_side_a( - ::AbstractContactFormulation, - X_a, X_b, U_a, U_b, - max_overlap_dist, - rel_smoothing_size -) - return ForwardDiff.gradient( - z -> contact_energy( - X_a, X_b, z, U_b, - max_overlap_dist, rel_smoothing_size - ), U_a - ) -end - -function contact_residual_side_b( - ::AbstractContactFormulation, - X_a, X_b, U_a, U_b, - max_overlap_dist, - rel_smoothing_size -) - return ForwardDiff.gradient( - z -> contact_energy( - X_a, X_b, U_a, z, - max_overlap_dist, rel_smoothing_size - ), U_b - ) -end - -function contact_stiffness_side_a_side_a( - ::AbstractContactFormulation, - X_a, X_b, U_a, U_b, - max_overlap_dist, - rel_smoothing_size -) - return ForwardDiff.hessian( - z -> contact_energy( - X_a, X_b, z, U_b, - max_overlap_dist, rel_smoothing_size - ), U_a - ) -end +# only for 2D +# assumes.... +# closest point projection distance +include("AD.jl") +include("ContactSurface.jl") +include("ContactPair.jl") -function contact_stiffness_side_b_side_b( - ::AbstractContactFormulation, - X_a, X_b, U_a, U_b, - max_overlap_dist, - rel_smoothing_size -) - return ForwardDiff.hessian( - z -> contact_energy( - X_a, X_b, U_a, z, - max_overlap_dist, rel_smoothing_size - ), U_b - ) -end +include("ClosestPointProjection.jl") +include("Search.jl") -function contact_stiffness_side_a_side_b( - ::AbstractContactFormulation, - X_a, X_b, U_a, U_b, - max_overlap_dist, - rel_smoothing_size -) - return ForwardDiff.jacobian( - z -> ForwardDiff.gradient( - contact_residual_side_a( - X_a, X_b, U_a, z, - max_overlap_dist, rel_smoothing_size - ), U_b - ) - ) -end +include("PenaltyContact.jl") +include("Mortar.jl") # TODO contact_stiffness_side_a_side_b # need to be careful when facet element types differ # used to construct these guys -function create_contact_pair_caches(mesh, contact_pairs) +function create_contact_pair_caches( + mesh, dof, contact_pairs; + kwargs... +) @assert num_dimensions(mesh) == 2 "Contact is only supported in 2D currently" syms = map(x -> Symbol("contact_pair_$(x.side_a)_$(x.side_b)"), contact_pairs) - caches = map(x -> ContactPairCache(mesh, x), contact_pairs) + caches = map(x -> ContactPairCache(mesh, dof, :displ, x; kwargs...), contact_pairs) return NamedTuple{tuple(syms...)}(tuple(caches...)) end diff --git a/src/contact/ContactSurface.jl b/src/contact/ContactSurface.jl index cee03a4..5d89a19 100644 --- a/src/contact/ContactSurface.jl +++ b/src/contact/ContactSurface.jl @@ -37,7 +37,7 @@ end @inline function facet_coordinates(surf::ContactSurface, X, U, id::Int) RT = eltype(X) ND = size(X, 1) - NNPF = 2 # hardcoded for now + NNPF = 2 # hardcoded for now TODO coords = SMatrix{ND, NNPF, RT, ND * NNPF}( @views X[:, surf.side_nodes[:, id]] @@ -52,7 +52,7 @@ end @inline function facet_coordinates(surf::ContactSurface, X, U, node::Int, facet::Int) RT = eltype(X) ND = size(X, 1) - NNPF = 2 # hardcoded for now + NNPF = 2 # hardcoded for now TODO coords = SVector{ND, RT}( @views X[:, surf.side_nodes[node, facet]] diff --git a/src/contact/Mortar.jl b/src/contact/Mortar.jl index 57e6d17..49d5996 100644 --- a/src/contact/Mortar.jl +++ b/src/contact/Mortar.jl @@ -43,16 +43,16 @@ function assemble_contact_matrix!( end function assemble_contact_scalar!( + energy, m::MortarContact, type::PenaltyContact, cache::ContactPairCache, X, U ) + block_energy = zero(eltype(X)) side_a = cache.side_a side_b = cache.side_b interactions = cache.interactions - - energy = zero(eltype(X)) for b in axes(side_b.side_nodes, 2) for a in axes(interactions, 1) index_a = interactions[a, b] @@ -68,15 +68,18 @@ function assemble_contact_scalar!( X_b = facet_field(side_b, X, b) U_a = facet_field(side_a, U, index_a) U_b = facet_field(side_b, U, b) - energy += contact_energy(type, X_a, X_b, U_a, U_b, m.max_overlap_dist, m.rel_smoothing_size) + block_energy += contact_energy(type, X_a, X_b, U_a, U_b, m.max_overlap_dist, m.rel_smoothing_size) end # @assert false end end - return energy + # fill!(contact_energy, energy) + energy .+= block_energy + return nothing end function assemble_contact_vector!( + vector, m::MortarContact, type::PenaltyContact, cache::ContactPairCache, @@ -85,7 +88,7 @@ function assemble_contact_vector!( side_a = cache.side_a side_b = cache.side_b interactions = cache.interactions - + n_dofs = size(U, 1) energy = zero(eltype(X)) for b in axes(side_b.side_nodes, 2) for a in axes(interactions, 1) @@ -97,19 +100,31 @@ function assemble_contact_vector!( # do nothing for this case @warn "Have an invalid index" else - X_a = facet_field(side_a, X, a) + X_a = facet_field(side_a, X, index_a) X_b = facet_field(side_b, X, b) - U_a = facet_field(side_a, U, a) + U_a = facet_field(side_a, U, index_a) U_b = facet_field(side_b, U, b) - # energy += contact_energy(X_a, X_b, U_a, U_b, m.max_overlap_dist, m.rel_smoothing_size) R_d = contact_residual_side_a(type, X_a, X_b, U_a, U_b, m.max_overlap_dist, m.rel_smoothing_size) R_f = contact_residual_side_b(type, X_a, X_b, U_a, U_b, m.max_overlap_dist, m.rel_smoothing_size) - # display(R_d) - # display(R_f) - # R_e = R_d + R_f - display(R_e) + + # now assemble the two contributions + local_dof = 1 + for node in axes(side_a.side_nodes, 1) + for dof in 1:n_dofs + vector[dof, side_a.side_nodes[node, index_a]] = R_d.data[local_dof] + local_dof = local_dof + 1 + end + end + + local_dof = 1 + for node in axes(side_b.side_nodes, 1) + for dof in 1:n_dofs + vector[dof, side_b.side_nodes[node, b]] = R_f.data[local_dof] + local_dof = local_dof + 1 + end + end end end end return energy -end \ No newline at end of file +end diff --git a/src/contact/PenaltyContact.jl b/src/contact/PenaltyContact.jl index 3cace82..47a4819 100644 --- a/src/contact/PenaltyContact.jl +++ b/src/contact/PenaltyContact.jl @@ -24,102 +24,81 @@ function contact_energy( # TODO move below to method l_a = norm(X_a[SVector{2, Int}(1, 2)] - X_a[SVector{2, Int}(3, 4)]) - integral = _integrate_gap(xi_a, g, max_overlap_dist) + integral = integrate_gap(xi_a, g, max_overlap_dist) return l_a * integral * _smooth_heaviside_at_zero(scaling, 1.) end -function _integrate_gap(xi, g, delta) - return _integrate_normalized_gap(xi, g / delta) -end - -# need these otherwise we allocate for some weird reason -# maybe we should use @setindex tools -_set_first(g::SVector{2, RT}, v::RT) where RT = SVector{2, RT}(v, g[2]) -_set_second(g::SVector{2, RT}, v::RT) where RT = SVector{2, RT}(g[1], v) - -function _integrate_normalized_gap(xi::SVector{2, RT}, g::SVector{2, RT}) where RT <: Number - g - dxi = xi[2] - xi[1] - g1_large = g[1] >= one(RT) - g2_large = g[2] >= one(RT) - dxi_invalid = (dxi <= 1e-14) || (g1_large && g2_large) - - # a = (g[2] - g[1]) / jnp.where(dxi_invalid, 1.0, dxi) - a = (g[2] - g[1]) - if !dxi_invalid - a = a / dxi +function _smooth_heaviside_at_zero(x, eps) + r = x / eps + # return jnp.where(x < eps, -2*r*r*r+3*r*r, 1.0) + if x < eps + return -2. * r * r * r + 3. * r * r + else + return 1. end +end - # g = jnp.where(g0_large, g.at[0].set(1.0), g) - if g1_large - # g = SVector{2, RT}(1., g[2]) - g = _set_first(g, one(eltype(g))) - end - # g = jnp.where(g1_large, g.at[1].set(1.0), g) - if g2_large - # g = SVector{2, RT}(g[1], 1.) - g = _set_second(g, one(eltype(g))) - end +@inline function integrate_normalized_gap( + ξ::AbstractVector{T}, + g::AbstractVector{T} +) where {T} + dξ = ξ[2] - ξ[1] + g0_large = g[1] >= one(eltype(g)) + g1_large = g[2] >= one(eltype(g)) + dξ_invalid = (dξ <= 1e-14) || (g0_large && g1_large) + + a = (g[2] - g[1]) / ifelse(dξ_invalid, one(eltype(g)), dξ) + g = ifelse(g0_large, typeof(g)(one(eltype(g)), g[2]), g) + g = ifelse(g1_large, typeof(g)(g[1], one(eltype(g))), g) a_is_zero = abs(g[1] - g[2]) < 1e-14 - # a = jnp.where(a_is_zero, 1.0, a) - if a_is_zero - a = 1. - end - - g1_old = g[1] - g2_old = g[2] - - # xi = jnp.where(g0_large & ~a_is_zero, - # xi.at[0].set((1.-g1_old) / a + xi[1]), - # xi) - - if g1_large && !a_is_zero - xi = typeof(xi)((1. - g2_old) / a + xi[2], xi[2]) - end - - # xi = jnp.where(g1_large & ~a_is_zero, - # xi.at[1].set((1.0-g0_old) / a + xi[0]), - # xi) - if g2_large && !a_is_zero - xi = typeof(xi)(xi[1], (1. - g1_old) / a + xi[1]) + a = ifelse(a_is_zero, one(eltype(a)), a) + g0_old = g[1] + g1_old = g[2] + + ξ = ifelse( + g0_large && !a_is_zero, + typeof(ξ)((1. - g1_old) / a + ξ[2], ξ[2]), + ξ + ) + ξ = ifelse( + g1_large && !a_is_zero, + typeof(ξ)(ξ[1], (1. - g0_old) / a + ξ[1]), + ξ + ) + dξ = ξ[2] - ξ[1] + + ξ0 = ξ[1] + ξ1 = ξ[2] + + safe_log = if g[2] / g[1] < 0. + 0. + else + log(g[2] / g[1]) end - dxi = xi[2] - xi[1] - - xi1 = xi[1] - xi2 = xi[2] - - intgl = log(g[2] / g[1]) / a + (0.5 * a * (xi2 * xi2 - xi1 * xi1) + g[1] * xi2 - g[2] * xi1) - 2.0 * dxi - + intgl = safe_log / a + + (0.5 * a * (ξ1 * ξ1 - ξ0 * ξ0)) + + g[1] * ξ1 - g[2] * ξ0 - + 2. * dξ + gbar = 0.5 * (g[1] + g[2]) dg = g[2] - g[1] - # intgl = jnp.where(jnp.abs(g[0]-g[1]) < 2.0e-4 * gbar, (1.0/gbar + gbar - 2.0 + dg * dg / (12.0 * gbar*gbar*gbar)) * dxi, intgl) - - if abs(g[1] - g[2]) < 2.e-4 * gbar - intgl = (1.0 / gbar + gbar - 2.0 + dg * dg / (12.0 * gbar * gbar * gbar)) * dxi - end - - # any_neg = any(g <= 0) - any_neg = any(x -> x <= zero(eltype(g)), g) - # intgl = jnp.where(any_neg, jnp.inf, intgl) - if any_neg - intgl = Inf - end - # intgl = jnp.where(dxi_invalid, 0.0 * intgl, intgl) - if dxi_invalid - intgl = 0. * intgl - end - + intgl = ifelse( + abs(g[1] - g[2]) < 2.e-4 * gbar, + (1.0 / gbar + gbar - 2.0 + dg * dg / (12.0 * gbar * gbar * gbar)) * dξ, + intgl + ) + any_neg = any(x -> x <= 0, g) + intgl = ifelse(any_neg, Inf, intgl) + intgl = ifelse(dξ_invalid, 0. * intgl, intgl) return intgl end -function _smooth_heaviside_at_zero(x, eps) - r = x / eps - # return jnp.where(x < eps, -2*r*r*r+3*r*r, 1.0) - if x < eps - return -2. * r * r * r + 3. * r * r - else - return 1. - end +@inline function integrate_gap( + xi::AbstractVector{T}, + g::AbstractVector{T}, + delta +) where {T} + return integrate_normalized_gap(xi, g ./ delta) end diff --git a/src/objectives/ConstrainedObjective.jl b/src/objectives/ConstrainedObjective.jl new file mode 100644 index 0000000..8b000e7 --- /dev/null +++ b/src/objectives/ConstrainedObjective.jl @@ -0,0 +1,43 @@ +struct ConstrainedObjective{ + C, + O <: AbstractObjective +} <: AbstractObjective{O} + constraint_func::C + objective_func::O +end + +struct ConstrainedObjectiveCache{ + A, O, RT, RV, + ObjCache <: AbstractObjectiveCache{A, O, RT, RV}, + F <: Function, CP +} <: AbstractObjectiveCache{A, O, RT, RV} + objective_cache::ObjCache + constraint_func::F + contact_pairs::CP +end + +function _penalty(c, λ, κ) + return ifelse( + λ >= κ * c, + -c * λ + 0.5 * κ * c * c, + -0.5 * λ * λ / κ + ) +end + +function _constraint_value!( + val, contact_caches, + constraint_func, u, p, λ, κ +) + c = constraint_func(contact_caches, u, p) + # c .= _penalty() +end + +function value(objective_cache, u, p, λ, κ) + c = objective_cache.constraint_func(objective_cache.contact_pairs, u, p) + c .= _penalty.(c, λ, κ) + return value(objective_cache.objective_cache, u, p) + sum(c) +end + +# function gradient(objective_cache, u, p, λ, κ) +# return gradient(objective_cache.objective_cache, u, p) +# end diff --git a/src/objectives/ContactObjective.jl b/src/objectives/ContactObjective.jl new file mode 100644 index 0000000..154faaa --- /dev/null +++ b/src/objectives/ContactObjective.jl @@ -0,0 +1,58 @@ +struct PenaltyContactObjectiveCache{ + F, C, + A, O, RT, RV, + P <: AbstractObjectiveCache{A, O, RT, RV} +} <: AbstractObjectiveCache{A, O, RT, RV} + contact_formulation::F + contact_pairs::C + physics_cache::P +end + +FiniteElementContainers.create_unknowns(o::PenaltyContactObjectiveCache) = +create_unknowns(o.physics_cache) + +assembler(o::PenaltyContactObjectiveCache) = assembler(o.physics_cache) + +function gradient(o::PenaltyContactObjectiveCache, U, p) + g = gradient(o.physics_cache, U, p) + c = create_field(o.physics_cache) + assemble_contact_vector!(c, o.contact_formulation, PenaltyContact(), o.contact_pairs, p.h1_coords, p.h1_field) + return g .+ 1e3 * c.data +end + +function hessian(o::PenaltyContactObjectiveCache, U, p; symmetric = true) + # TODO eventually tack on contact contribution + return hessian(o.physics_cache, U, p; symmetric = symmetric) +end + +function hvp(o::PenaltyContactObjectiveCache, U, V, p) + return hvp(o.physics_cache, U, V, p) +end + +function value(o::PenaltyContactObjectiveCache, U, p) + e = value(o.physics_cache, U, p) + c = zeros(1) + assemble_contact_scalar!(c, o.contact_formulation, PenaltyContact(), o.contact_pairs, p.h1_coords, p.h1_field) + return e + 1e3 * c[1] +end + +function initialize!(o::PenaltyContactObjectiveCache, U, p) + initialize!(o.physics_cache, U, p) + return nothing +end + +function step!(solver, o::PenaltyContactObjectiveCache, U, p; verbose = true) + FiniteElementContainers.update_time!(p) + FiniteElementContainers.update_bc_values!(p) + update_nearest_neighbors!(o.contact_pairs, p.h1_coords, U) + _step_begin_banner(o.physics_cache, p; verbose = verbose) + # do solve TODO + solve!(solver, U.data, p) + update_field_dirichlet_bcs!(U, p.dirichlet_bcs) + + # update values at end of step + gradient(o, U, p) + value(o, U, p) + # update gradient/vale at end of step TODO + +end diff --git a/src/objectives/DesignObjective.jl b/src/objectives/DesignObjective.jl index 8f94d03..4910fec 100644 --- a/src/objectives/DesignObjective.jl +++ b/src/objectives/DesignObjective.jl @@ -1,8 +1,9 @@ struct DesignObjective{ F <: Function, - Q, + Q <: AbstractArray, P, - S + S <: AbstractArray{<:AbstractField, 1}, + QS } func::F qois::Q @@ -11,20 +12,65 @@ struct DesignObjective{ dstored_parameters::P stored_solutions::S dstored_solutions::S + qoi_storages::QS + dqoi_storages::QS end + function DesignObjective(func, qois, U, p) dqois = Enzyme.make_zero(qois) stored_parameters = typeof(p)[] dstored_parameters = typeof(p)[] stored_solutions = typeof(U)[] dstored_solutions = typeof(U)[] + + # if this works, we'll likely need to reconcile + # how to deal with multiple qois + qoi_storages = typeof(qois[1].storage)[] + dqoi_storages = typeof(qois[1].storage)[] + + @show typeof(qoi_storages) return DesignObjective( func, qois, dqois, stored_parameters, dstored_parameters, - stored_solutions, dstored_solutions + stored_solutions, dstored_solutions, + qoi_storages, dqoi_storages ) end +function forward_problem!(obj::DesignObjective, solver, objective_cache, U, p) + empty!(obj.stored_parameters) + empty!(obj.dstored_parameters) + empty!(obj.stored_solutions) + empty!(obj.dstored_solutions) + empty!(obj.qoi_storages) + empty!(obj.dqoi_storages) + initialize!(objective_cache, U, p) + + time_end = sum(p.times.time_end) + + while FiniteElementContainers.current_time(p.times) < time_end - 1e3 * eps(time_end) + step!(solver, objective_cache, U, p; verbose=true) + # might need to reset BCs here + FiniteElementContainers.update_field_dirichlet_bcs!(U, p.dirichlet_bcs) + + temp = similar(U) + # copyto!(temp.data, U.data) + temp.data .= U.data + # push!(obj.stored_solutions, deepcopy(U)) + push!(obj.stored_solutions, temp) + push!(obj.stored_parameters, deepcopy(p)) + end + + # now resize stored qoi storages + # for val in obj.qoi_storages + # TODO currently will only work with one QOI + for _ in 1:length(obj.stored_solutions) + push!(obj.qoi_storages, Enzyme.make_zero(obj.qois[1].storage)) + push!(obj.dqoi_storages, Enzyme.make_zero(obj.qois[1].storage)) + end + return nothing +end + function _gradient_and_value_init!(obj::DesignObjective, ::Enzyme.ReverseMode) Enzyme.make_zero!(obj.dqois) @@ -35,13 +81,18 @@ function _gradient_and_value_init!(obj::DesignObjective, ::Enzyme.ReverseMode) for val in obj.stored_solutions push!(obj.dstored_solutions, Enzyme.make_zero(val)) end + + for val in obj.dqoi_storages + Enzyme.make_zero!(val) + end return nothing end function _gradient_and_value!( design_obj, f, df, - qois, dqois, + qoi_storages, dqoi_storages, + qois, Us, dUs, ps, dps, params, dparams @@ -50,7 +101,8 @@ function _gradient_and_value!( Reverse, design_obj, Duplicated(f, df), - Duplicated(qois, dqois), + Duplicated(qoi_storages, dqoi_storages), + Const(qois), Duplicated(Us, dUs), Duplicated(ps, dps), Duplicated(params, dparams) @@ -64,7 +116,8 @@ function gradient_and_value(obj::DesignObjective, params) _gradient_and_value_init!(obj, Reverse) _gradient_and_value!( obj.func, f, df, - obj.qois, obj.dqois, + obj.qoi_storages, obj.dqoi_storages, + obj.qois, obj.stored_solutions, obj.dstored_solutions, obj.stored_parameters, obj.dstored_parameters, params, dparams @@ -73,7 +126,7 @@ function gradient_and_value(obj::DesignObjective, params) end function value!(f, obj::DesignObjective, params) - obj.func(f, obj.qois, obj.stored_solutions, obj.stored_parameters, params) + obj.func(f, obj.qoi_storages, obj.qois, obj.stored_solutions, obj.stored_parameters, params) return nothing end diff --git a/src/objectives/ExplicitDynamicsObjective.jl b/src/objectives/ExplicitDynamicsObjective.jl index ee40f44..cd6c02d 100644 --- a/src/objectives/ExplicitDynamicsObjective.jl +++ b/src/objectives/ExplicitDynamicsObjective.jl @@ -1,14 +1,23 @@ +struct ExplicitDynamicsObjective{ + F1 <: Function, + F2 <: Function, + F3 <: Function +} <: AbstractObjective{F1} + value::F1 + gradient_u::F2 + hessian_u::F3 +end + function ExplicitDynamicsObjective() - return QuadratureLevelObjective(energy, residual, stiffness) + return ExplicitDynamicsObjective(energy, residual, stiffness) end -struct ExplicitDynamicsObjectiveCache{ - A, O, P, +mutable struct ExplicitDynamicsObjectiveCache{ + A, O, RT, RV <: AbstractArray{RT, 1}, NF -} <: AbstractObjectiveCache{A, O, P, RT, RV} +} <: AbstractObjectiveCache{A, O, RT, RV} assembler::A objective::O - parameters::P # β::RT γ::RT @@ -20,9 +29,6 @@ struct ExplicitDynamicsObjectiveCache{ external_force::H1Field{RT, RV, NF} inertial_force::H1Field{RT, RV, NF} internal_force::H1Field{RT, RV, NF} - predictor_solution::H1Field{RT, RV, NF} - predictor_solution_rate::H1Field{RT, RV, NF} - solution::H1Field{RT, RV, NF} solution_rate::H1Field{RT, RV, NF} solution_rate_rate::H1Field{RT, RV, NF} solution_old::H1Field{RT, RV, NF} @@ -31,20 +37,29 @@ struct ExplicitDynamicsObjectiveCache{ value::RV gradient::H1Field{RT, RV, NF} lumped_hessian::RV + lumped_mass::RV + # + step_number::Int # timer::TimerOutput end -function ExplicitDynamicsObjectiveCache(sim, CFL, β=0.25, γ=0.5) - objective = ImplicitDynamicsObjectiveNew() - assembler, parameters = _setup_simulation_common( - sim, nothing; - return_post_processor=false, - use_condensed=true - ) +function ExplicitDynamicsObjectiveCache( + assembler, + objective::ExplicitDynamicsObjective, + CFL, β=0.25, γ=0.5 +) + # objective = ImplicitDynamicsObjectiveNew() + # assembler, parameters = _setup_simulation_common( + # sim, nothing; + # return_post_processor=false, + # use_condensed=true + # ) - RT = eltype(parameters.h1_coords) - backend = KA.get_backend(parameters.h1_coords) + # RT = eltype(parameters.h1_coords) + # backend = KA.get_backend(parameters.h1_coords) + RT = eltype(assembler.constraint_storage) + backend = KA.get_backend(assembler) external_energy = KA.zeros(backend, RT, 1) internal_energy = KA.zeros(backend, RT, 1) @@ -54,10 +69,6 @@ function ExplicitDynamicsObjectiveCache(sim, CFL, β=0.25, γ=0.5) inertial_force = create_field(assembler) internal_force = create_field(assembler) - predictor_solution = create_field(assembler) - predictor_solution_rate = create_field(assembler) - - solution = create_field(assembler) solution_rate = create_field(assembler) solution_rate_rate = create_field(assembler) @@ -67,19 +78,157 @@ function ExplicitDynamicsObjectiveCache(sim, CFL, β=0.25, γ=0.5) value = KA.zeros(backend, RT, 1) gradient = create_field(assembler) lumped_hessian = KA.zeros(backend, RT, length(gradient)) + lumped_mass = KA.zeros(backend, RT, length(gradient)) timer = TimerOutput() return ExplicitDynamicsObjectiveCache( - assembler, objective, parameters, + assembler, objective, β, γ, CFL, external_energy, internal_energy, kinetic_energy, external_force, inertial_force, internal_force, - predictor_solution, predictor_solution_rate, - solution, solution_rate, solution_rate_rate, + solution_rate, solution_rate_rate, solution_old, solution_rate_old, - value, gradient, lumped_hessian, - timer + value, gradient, + lumped_hessian, lumped_mass, + 0, timer ) end +function setup_cache(assembler, objective::ExplicitDynamicsObjective, args...) + return ExplicitDynamicsObjectiveCache(assembler, objective, args...) +end + +function initialize!( + o::ExplicitDynamicsObjectiveCache, U, p; + lumped_mass_style = :row_sum, + displ_ics = nothing, + vel_ics = nothing, + acc_ics = nothing +) + fill!(p.times.time_current, zero(eltype(p.times.time_current))) + + M = mass_matrix(o, U, p) + + # need to lump the mass + if lumped_mass_style == :row_sum + n = size(M, 1) + d = vec(sum(M, dims=2)) + o.lumped_mass .= d + elseif lumped_mass_style == :diagonal_extraction + # return spdiagm(0 => diag(M)) + o.lumped_mass .= diag(M) + elseif lumped_mass_style == :scaled_diagonal + d = diag(M) + total_mass = sum(M) + scale = total_mass / sum(d) + # o.lumped_mass.= spdiagm(0 => scale .* d) + o.lumped_mass .= scale .* d + elseif lumped_mass_style == :block_lumping + @assert ndim == 2 || ndim == 3 + ndofs = size(M, 1) + @assert ndofs % ndim == 0 + + nnodes = div(ndofs, ndim) + Ml = zeros(Float64, ndofs) + + # Row-sum per dof + rowsums = vec(sum(M, dims=2)) + + # Accumulate nodal mass + for a in 1:nnodes + idx = (a-1)*ndim + 1 : a*ndim + m_node = sum(rowsums[idx]) + Ml[idx] .= m_node / ndim + end + + o.lumped_mass .= spdiagm(0 => Ml) + else + @assert false + end + + # apply initial conditions + # TODO make this better in FEC + if displ_ics !== nothing + FiniteElementContainers.update_ic_values!(displ_ics, p.h1_coords) + FiniteElementContainers.update_field_ics!(U, displ_ics) + FiniteElementContainers.update_field_ics!(o.solution_old, displ_ics) + end + + if vel_ics !== nothing + FiniteElementContainers.update_ic_values!(vel_ics, p.h1_coords) + FiniteElementContainers.update_field_ics!(o.solution_rate, vel_ics) + FiniteElementContainers.update_field_ics!(o.solution_rate_old, vel_ics) + end + + if acc_ics !== nothing + FiniteElementContainers.update_ic_values!(acc_ics, p.h1_coords) + FiniteElementContainers.update_field_ics!(o.solution_rate_rate, acc_ics) + end + + # calculate mechanics stuff at beginning + internal_energy = _internal_energy(o, U, p) + kinetic_energy = 0.5 * dot(o.lumped_mass, o.solution_rate.data .* o.solution_rate.data) + fill!(o.internal_energy, internal_energy) + fill!(o.kinetic_energy, kinetic_energy) + + internal_force = _internal_force(o, U, p) + o.internal_force .= internal_force + + # TODO external force + external_force = create_field(o.assembler) + o.external_force .= external_force + o.inertial_force .= internal_force .- external_force + o.solution_rate_rate.data .= o.inertial_force.data ./ o.lumped_mass +end + +function _internal_energy(o::ExplicitDynamicsObjectiveCache, U, p) + assemble_scalar!(assembler(o), o.objective.value, U, p) + val = mapreduce(sum, +, values(assembler(o).scalar_quadrature_storage)) + fill!(o.internal_energy, val) + # assemble_scalar!(assembler(o), #neumann energy, U, p) + + # TODO need an actual neumann energy method in FEContainers + fill!(o.external_energy, dot(o.external_force, U)) + o.value .= o.external_energy .+ o.internal_energy + return sum(o.value) +end + +function _internal_force(o::ExplicitDynamicsObjectiveCache, U, p) + assemble_vector!(assembler(o), o.objective.gradient_u, U, p) + # return assembler(o).residual_storage + o.internal_force .= assembler(o).residual_storage +end + +function mass_matrix(o::ExplicitDynamicsObjectiveCache, U, p) + # assemble_mass!(assembler(o), o.objective.hessian_u, U, p) + # M = FiniteElementContainers.mass(assembler(o)) |> Symmetric + # really passing displacement below instead of velocity to save + # on having to store the solution_rate for quasi-statics + # this method is mainly so you can construct a mass matrix + # if you want one. Should add abstract types in the future + assemble_mass!(assembler(o), mass, U, p) + M = FiniteElementContainers.mass(assembler(o)) + return M #|> Symmetric +end + +function step!(solver, o::ExplicitDynamicsObjectiveCache, U, p; verbose = true) + FiniteElementContainers.update_time!(p) + FiniteElementContainers.update_bc_values!(p) + solve!(solver, U, p) + _step_log(o, p) +end + +function _step_log(o::ExplicitDynamicsObjectiveCache, p) + if o.step_number % 1000 == 0 + @info "Step Time Internal External Kinetic" + @info "Number Increment Energy Energy Energy" + end + + if o.step_number % 100 == 0 + str = @sprintf "%8d %.4e %.4e %.4e %.4e" o.step_number p.times.Δt[1] o.internal_energy[1] o.external_energy[1] o.kinetic_energy[1] + @info str + end + + o.step_number = o.step_number + 1 +end diff --git a/src/objectives/Objectives.jl b/src/objectives/Objectives.jl index 6f9ce2e..814cddc 100644 --- a/src/objectives/Objectives.jl +++ b/src/objectives/Objectives.jl @@ -1,9 +1,7 @@ """ User facing """ -abstract type AbstractObjective{ - F1 <: Function -} end +abstract type AbstractObjective{F1} end # new cache implementation abstract type AbstractObjectiveCache{ @@ -29,7 +27,9 @@ end # function parameters(o::AbstractObjectiveCache) # return o.parameters # end - +include("ContactObjective.jl") +include("ConstrainedObjective.jl") include("DesignObjective.jl") +include("ExplicitDynamicsObjective.jl") # include("ImplicitDynamicsObjective.jl") include("QuasiStaticObjective.jl") diff --git a/src/objectives/QuasiStaticObjective.jl b/src/objectives/QuasiStaticObjective.jl index 3c820de..7c85d33 100644 --- a/src/objectives/QuasiStaticObjective.jl +++ b/src/objectives/QuasiStaticObjective.jl @@ -13,21 +13,17 @@ function QuasiStaticObjective() end struct QuasiStaticObjectiveCache{ - # A, O, P, A, O, RT, RV <: AbstractArray{RT, 1}, NF } <: AbstractObjectiveCache{A, O, RT, RV} -# } <: AbstractObjectiveCache{A, O, P, RT, RV} assembler::A objective::O - # parameters::P # external_energy::RV internal_energy::RV external_force::H1Field{RT, RV, NF} internal_force::H1Field{RT, RV, NF} - # solution::H1Field{RT, RV, NF} solution_old::H1Field{RT, RV, NF} # solver helpers value::RV @@ -36,36 +32,17 @@ struct QuasiStaticObjectiveCache{ timer::TimerOutput end -# function QuasiStaticObjectiveCache( -# objective::QuasiStaticObjective, -# sim; -# kwargs... -# ) -# # objective = QuasiStaticObjective() -# # assembler, parameters = _setup_simulation_common( -# # sim, nothing; -# # # return_post_processor=false, -# # use_condensed=true, -# # kwargs... -# # ) -# assembler = _setup_assembler( -# sim; -# use_condensed=true, -# kwargs... -# ) function QuasiStaticObjectiveCache( assembler, objective::QuasiStaticObjective ) RT = eltype(assembler.constraint_storage) - # backend = KA.get_backend(parameters.h1_coords) backend = KA.get_backend(assembler) external_energy = KA.zeros(backend, RT, 1) internal_energy = KA.zeros(backend, RT, 1) external_force = create_field(assembler) internal_force = create_field(assembler) - # solution = create_field(assembler) solution_old = create_field(assembler) value = KA.zeros(backend, RT, 1) @@ -74,10 +51,9 @@ function QuasiStaticObjectiveCache( timer = TimerOutput() return QuasiStaticObjectiveCache( - assembler, objective, #parameters, + assembler, objective, external_energy, internal_energy, external_force, internal_force, - # solution, solution_old, solution_old, value, gradient, timer @@ -85,14 +61,10 @@ function QuasiStaticObjectiveCache( end # cache hook -function setup_cache(assembler, objective) +function setup_cache(assembler, objective::QuasiStaticObjective) return QuasiStaticObjectiveCache(assembler, objective) end -# function setup_cache(objective::O, sim; kwargs...) where O <: QuasiStaticObjective -# return QuasiStaticObjectiveCache(objective, sim; kwargs...) -# end - # objective hooks function gradient(o::QuasiStaticObjectiveCache, U, p) RT = eltype(U) @@ -128,9 +100,13 @@ function gradient(o::QuasiStaticObjectiveCache, U, p) end end -function hessian(o::QuasiStaticObjectiveCache, U, p) +function hessian(o::QuasiStaticObjectiveCache, U, p; symmetric = true) assemble_stiffness!(assembler(o), o.objective.hessian_u, U, p) - H = stiffness(assembler(o)) |> Symmetric + if symmetric + H = stiffness(assembler(o)) |> Symmetric + else + H = stiffness(assembler(o)) + end return H end @@ -153,7 +129,7 @@ end function value(o::QuasiStaticObjectiveCache, U, p) assemble_scalar!(assembler(o), o.objective.value, U, p) - val = mapreduce(sum, sum, values(assembler(o).scalar_quadrature_storage)) + val = mapreduce(sum, +, values(assembler(o).scalar_quadrature_storage)) fill!(o.internal_energy, val) # assemble_scalar!(assembler(o), #neumann energy, U, p) @@ -164,35 +140,24 @@ function value(o::QuasiStaticObjectiveCache, U, p) end # integrator hooks -# function initialize!(o::QuasiStaticObjectiveCache, sim) -# # p = o.parameters -# p = create_parameters( -# mesh, -# assembler, sim.physics, sim.properties; -# dirichlet_bcs=sim.dirichlet_bcs, -# neumann_bcs=sim.neumann_bcs, -# times=sim.times -# ) -# fill!(p.times.time_current, zero(eltype(p.times.time_current))) -# return nothing -# end function initialize!(o::QuasiStaticObjectiveCache, U, p) fill!(p.times.time_current, zero(eltype(p.times.time_current))) return nothing end -# function step!(o::QuasiStaticObjectiveCache, solver; verbose=true) function step!(solver, o::QuasiStaticObjectiveCache, U, p; verbose = true) FiniteElementContainers.update_time!(p) FiniteElementContainers.update_bc_values!(p) _step_begin_banner(o, p; verbose = verbose) solve!(solver, U.data, p) + update_field_dirichlet_bcs!(U, p.dirichlet_bcs) # update values at end of step gradient(o, U, p) value(o, U, p) # update old solution + update_field_dirichlet_bcs!(U, p.dirichlet_bcs) copyto!(o.solution_old, U) _step_end_banner(o, p; verbose = verbose) @@ -201,7 +166,6 @@ end # logging hooks function _step_begin_banner(::QuasiStaticObjectiveCache, p; verbose::Bool = true) - # p = o.parameters if verbose time_curr = FiniteElementContainers.current_time(p.times) time_start = sum(p.times.time_start) diff --git a/src/physics/MaterialOutput.jl b/src/physics/MaterialOutput.jl index d1173e1..fc4a239 100644 --- a/src/physics/MaterialOutput.jl +++ b/src/physics/MaterialOutput.jl @@ -6,6 +6,7 @@ struct StandardMaterialOutput{ algorithmic_energy::RT cauchy_stress::SymmetricTensor{2, 3, RT, 6} displacement_gradient::Tensor{2, 3, RT, 9} + # state_variables::Svector{NS, RT} end @inline function standard_material_output( @@ -36,9 +37,10 @@ function update_material_output!( objective_cache, U, p ) FiniteElementContainers.assemble_quadrature_quantity!( - mat_output, assembler(objective_cache).dof, + mat_output, nothing, assembler(objective_cache).dof, standard_material_output, - U, p + U, p, + FiniteElementContainers.AssembledStruct{StandardMaterialOutput}() ) return nothing end diff --git a/src/physics/SolidMechanics.jl b/src/physics/SolidMechanics.jl index a278b29..54d4470 100644 --- a/src/physics/SolidMechanics.jl +++ b/src/physics/SolidMechanics.jl @@ -13,7 +13,10 @@ function SolidMechanics(formulation, model) end function FiniteElementContainers.create_properties(physics::SolidMechanics, inputs) - return ConstitutiveModels.initialize_props(physics.constitutive_model, inputs) + density = inputs["density"] + mat_model_props = ConstitutiveModels.initialize_props(physics.constitutive_model, inputs) + mat_model_props = Array(mat_model_props) + return pushfirst!(mat_model_props, density) end function FiniteElementContainers.create_initial_state(physics::SolidMechanics) @@ -35,10 +38,12 @@ end # kinematics ∇u_q = modify_field_gradients(physics.formulation, ∇u_q) + mat_props = @views props_el[2:end] + # constitutive θ = 0.0 # TODO ψ_q = ConstitutiveModels.helmholtz_free_energy( - physics.constitutive_model, props_el, dt, ∇u_q, θ, state_old_q, state_new_q + physics.constitutive_model, mat_props, dt, ∇u_q, θ, state_old_q, state_new_q ) return JxW * ψ_q end @@ -55,10 +60,12 @@ end # kinematics ∇u_q = modify_field_gradients(physics.formulation, ∇u_q) + mat_props = @views props_el[2:end] + # constitutive θ = 0.0 # TODO ψ_q = func( - physics.constitutive_model, props_el, dt, ∇u_q, θ, state_old_q, state_new_q + physics.constitutive_model, mat_props, dt, ∇u_q, θ, state_old_q, state_new_q ) return ψ_q end @@ -74,10 +81,12 @@ end # kinematics ∇u_q = modify_field_gradients(physics.formulation, ∇u_q) + mat_props = @views props_el[2:end] + # constitutive θ = 0.0 # TODO ψ_q = func( - physics.constitutive_model, props_el, dt, ∇u_q, θ, state_old_q, state_new_q + physics.constitutive_model, mat_props, dt, ∇u_q, θ, state_old_q, state_new_q ) return JxW * ψ_q end @@ -92,10 +101,12 @@ end # kinematics ∇u_q = modify_field_gradients(physics.formulation, ∇u_q) + mat_props = @views props_el[2:end] + # constitutive θ = 0.0 # TODO P_q = ConstitutiveModels.pk1_stress( - physics.constitutive_model, props_el, dt, ∇u_q, θ, state_old_q, state_new_q + physics.constitutive_model, mat_props, dt, ∇u_q, θ, state_old_q, state_new_q ) # turn into voigt notation P_q = extract_stress(physics.formulation, P_q) @@ -114,10 +125,12 @@ end # kinematics ∇u_q = modify_field_gradients(physics.formulation, ∇u_q) + mat_props = @views props_el[2:end] + # constitutive θ = 0. # TODO A_q = ConstitutiveModels.material_tangent( - physics.constitutive_model, props_el, dt, ∇u_q, θ, state_old_q, state_new_q + physics.constitutive_model, mat_props, dt, ∇u_q, θ, state_old_q, state_new_q ) # turn into voigt notation K_q = extract_stiffness(physics.formulation, A_q) @@ -134,7 +147,7 @@ end v_q = interpolate_field_values(physics, interps, v_el) # TODO - rho = 1. + rho = props_el[1] return 0.5 * JxW * rho * dot(v_q, v_q) end @@ -148,6 +161,6 @@ function mass( physics::SolidMechanics, interps, x_el, t, dt, v_el, v_el_old, state_old_q, state_new_q, props_el ) return ForwardDiff.hessian(z -> _kinetic_energy( - physics, interps, x_el, t, dt, z, u_el_old, state_old_q, state_new_q, props_el + physics, interps, x_el, t, dt, z, v_el_old, state_old_q, state_new_q, props_el ), v_el) end diff --git a/src/qoi_extractors/QOIExtractors.jl b/src/qoi_extractors/QOIExtractors.jl index f98373a..156e243 100644 --- a/src/qoi_extractors/QOIExtractors.jl +++ b/src/qoi_extractors/QOIExtractors.jl @@ -1,15 +1,8 @@ -# TODO -# plans for QOI extractors -# 1. have a base cache type to handle whether -# or not we need AD which is the main -# source of allocations for these structs -# 2. have input helper methods for what type of -# calc is needed, how to reduce, and what types -# to pre-allocate for storage - abstract type AbstractQOIExtractor end struct QOIExtractor{ + FieldQOI, + MatQOI, A <: FiniteElementContainers.AbstractAssembler, EvalFunc <: Function, Reduction1 <: Function, @@ -30,8 +23,12 @@ function QOIExtractor( cache_arrtype, cache_eltype; component_extractor = nothing, + is_field_qoi = false, + is_material_qoi = false, reduction_2 = identity ) + @assert is_field_qoi || is_material_qoi + if cache_arrtype isa Type{<:H1Field} @info "Requesting QOI output on reduction of H1Field" @assert cache_eltype isa Type{<:Number} @@ -60,55 +57,140 @@ function QOIExtractor( end # TODO it won't always be a mat func - if component_extractor !== nothing - function mat_func(x1, x2, x3, x4, x5, x6, x7, x8, x9, x10) - val = general_integrated_material_qoi(func, x1, x2, x3, x4, x5, x6, x7, x8, x9, x10) - return val[component_extractor...] + new_func = nothing + + if is_field_qoi + if component_extractor !== nothing + @assert false "Finish me" + else + # field_func = (x1, x2, x3, x4, x5, x6, x7, x8, x9, x10) + new_func = func end - else - mat_func = (x1, x2, x3, x4, x5, x6, x7, x8, x9, x10) -> - general_integrated_material_qoi(func, x1, x2, x3, x4, x5, x6, x7, x8, x9, x10) end - return QOIExtractor( + if is_material_qoi + if component_extractor !== nothing + function mat_func(x1, x2, x3, x4, x5, x6, x7, x8, x9, x10) + val = general_integrated_material_qoi(func, x1, x2, x3, x4, x5, x6, x7, x8, x9, x10) + return val[component_extractor...] + end + else + mat_func = (x1, x2, x3, x4, x5, x6, x7, x8, x9, x10) -> + general_integrated_material_qoi(func, x1, x2, x3, x4, x5, x6, x7, x8, x9, x10) + end + + new_func = mat_func + end + + return QOIExtractor{ + is_field_qoi, is_material_qoi, + typeof(objective_cache.assembler), + typeof(new_func), typeof(reduction_1), typeof(reduction_2), + typeof(storage) + }( objective_cache.assembler, - mat_func, reduction_1, reduction_2, + new_func, reduction_1, reduction_2, storage ) end -# TODO will need to specialize for GPU -function _value!( - f, storage, asm, func, U, p, reduction_1, reduction_2 -) - # TODO - # ideally below line is not necessary - # but we need to handle the differences between - # condensed or not in update_for_assembly! - # FiniteElementContainers.update_field_dirichlet_bcs!(U, p.dirichlet_bcs) - FiniteElementContainers.assemble_quadrature_quantity!( - storage, asm.dof, func, U, p - ) - f_temp = mapreduce(x -> reduce(reduction_1, x), reduction_2, values(storage)) +# TODO need to add component extractors for _perform_reductions! methods +function _perform_reductions!( + f, + qoi::QOIExtractor{false, true, A, F, R1, R2, S}, + storage +) where {A, F, R1, R2, S} + f_temp = mapreduce(x -> reduce(qoi.reduction_1, x), qoi.reduction_2, values(storage)) fill!(f, f_temp) return nothing end -function _value!( - f, qoi::QOIExtractor, U, p +function _perform_reductions!( + f, + qoi::QOIExtractor{true, false, A, F, R1, R2, S}, + storage +) where {A, F, R1, R2, S} + # f_temp = mapreduce(x -> reduce(qoi.reduction_1, x), qoi.reduction_2, values(qoi.storage)) + # fill!(f, f_temp) + # TODO do something useful here + # currently not actually do any reductions + copyto!(f.data, storage.data) + return nothing +end + +function _update_storage!( + asm_method, + storage, dof, func, U, p ) - FiniteElementContainers.assemble_quadrature_quantity!( + asm_method(storage, dof, func, U, p) + return nothing +end + +# deprecate +function _update_storage!( + qoi::QOIExtractor{false, true, A, F, R1, R2, S}, + U, p +) where {A, F, R1, R2, S} + _update_storage!( + FiniteElementContainers.assemble_quadrature_quantity!, qoi.storage, qoi.assembler.dof, qoi.func, U, p ) - f_temp = mapreduce(x -> reduce(qoi.reduction_1, x), qoi.reduction_2, values(qoi.storage)) - fill!(f, f_temp) return nothing end +# deprecate +function _update_storage!( + qoi::QOIExtractor{true, false, A, F, R1, R2, S}, + U, p +) where {A, F, R1, R2, S} + _update_storage!( + FiniteElementContainers.assemble_vector!, + qoi.storage, qoi.assembler.dof, qoi.func, U, p + ) + return nothing +end + +function _update_storage!( + storage, + qoi::QOIExtractor{false, true, A, F, R1, R2, S}, + U, p +) where {A, F, R1, R2, S} + _update_storage!( + FiniteElementContainers.assemble_quadrature_quantity!, + storage, qoi.assembler.dof, qoi.func, U, p + ) + return nothing +end + +function _update_storage!( + storage, + qoi::QOIExtractor{true, false, A, F, R1, R2, S}, + U, p +) where {A, F, R1, R2, S} + _update_storage!( + FiniteElementContainers.assemble_vector!, + storage, qoi.assembler.dof, qoi.func, U, p + ) + return nothing +end + +function value!(f, qoi::QOIExtractor, U, p) + _update_storage!(qoi, U, p) + _perform_reductions!(f, qoi, qoi.storage) + return nothing +end + +# for use in design objective +function value!(f, storage, qoi::QOIExtractor, U, p) + _update_storage!(storage, qoi, U, p) + _perform_reductions!(f, qoi, storage) + return nothing +end + +# TODO this out of place method will currently only +# work for scalar reductions function value(qoi::QOIExtractor, U, p) f = zeros(1) - # asm = assembler(qoi.objective_cache) - # _value!(f, qoi.storage, qoi.assembler, qoi.func, U, p, qoi.reduction_1, qoi.reduction_2) - _value!(f, qoi, U, p) + value!(f, qoi, U, p) return f[1] end diff --git a/src/simulations/Simulations.jl b/src/simulations/Simulations.jl index 3674019..f6c6823 100644 --- a/src/simulations/Simulations.jl +++ b/src/simulations/Simulations.jl @@ -17,7 +17,10 @@ abstract type AbstractSimulation end # abstract type AbstractSimulationCache end # function run!(sim, objective_type, solver_type) -function run!(solver, objective_cache, U, p, sim) +function run!( + solver, objective_cache, U, p, sim; + output_exodus_every = 1 +) pp = PostProcessor(objective_cache, U, p, sim) post_process(pp, objective_cache, U, p, 1) @@ -29,7 +32,10 @@ function run!(solver, objective_cache, U, p, sim) try while FiniteElementContainers.current_time(p.times) < time_end - 1e3 * eps(time_end) step!(solver, objective_cache, U, p; verbose=solver.settings.verbose) - post_process(pp, objective_cache, U, p, n) + + if n % output_exodus_every == 0 + post_process(pp, objective_cache, U, p, n) + end n = n + 1 end finally @@ -63,6 +69,7 @@ function _setup_assembler_and_parameters( mesh, assembler, sim.physics, sim.properties; dirichlet_bcs=sim.dirichlet_bcs, + ics=sim.solution_ics, neumann_bcs=sim.neumann_bcs, times=sim.times ) @@ -70,19 +77,13 @@ function _setup_assembler_and_parameters( return assembler, parameters end -# function setup_parameters( -# mesh -# ) - -# end - -function setup_caches(obj, sim; kwargs...) +function setup_caches(obj, sim, args...; kwargs...) asm, p = _setup_assembler_and_parameters( sim; use_condensed = true, kwargs... ) - objective_cache = setup_cache(asm, obj) + objective_cache = setup_cache(asm, obj, args...) U = create_field(asm) return objective_cache, U, p end diff --git a/src/simulations/SingleDomainSimulation.jl b/src/simulations/SingleDomainSimulation.jl index ab23405..7eefa69 100644 --- a/src/simulations/SingleDomainSimulation.jl +++ b/src/simulations/SingleDomainSimulation.jl @@ -3,6 +3,9 @@ struct SingleDomainSimulation{ RV <: AbstractArray{RT, 1}, P1 <: NamedTuple, P2 <: NamedTuple, + I1, + I2, + I3, D, N, C @@ -12,6 +15,9 @@ struct SingleDomainSimulation{ times::TimeStepper{RV} physics::P1 properties::P2 + solution_ics::I1 + solution_rate_ics::I2 + solution_rate_rate_ics::I3 dirichlet_bcs::D neumann_bcs::N contact_pairs::C @@ -23,6 +29,9 @@ function SingleDomainSimulation( times::TimeStepper, physics::NamedTuple, properties::NamedTuple; + solution_ics::Vector{<:InitialCondition} = InitialCondition[], + solution_rate_ics::Vector{<:InitialCondition} = InitialCondition[], + solution_rate_rate_ics::Vector{<:InitialCondition} = InitialCondition[], dirichlet_bcs::Vector{<:DirichletBC} = DirichletBC[], neumann_bcs::Vector{<:NeumannBC} = NeumannBC[], contact_pairs::Vector{<:ContactPair} = ContactPair[] @@ -32,6 +41,7 @@ function SingleDomainSimulation( return SingleDomainSimulation( mesh_file, output_file, times, physics, properties, + solution_ics, solution_rate_ics, solution_rate_rate_ics, dirichlet_bcs, neumann_bcs, contact_pairs ) end diff --git a/src/solvers/AugmentedLagrange.jl b/src/solvers/AugmentedLagrange.jl new file mode 100644 index 0000000..7da12ad --- /dev/null +++ b/src/solvers/AugmentedLagrange.jl @@ -0,0 +1,95 @@ +Base.@kwdef struct AugmentedLagrangeSolverSettings + penalty_scaling::Float64 = 4. + target_constraint_decrease_factor::Float64 = 0.75 + relative_gmres_tol::Float64 = 2.e-2 + max_gmres_iters::Int = 100 + use_second_order_update::Bool = true + use_newton_only::Bool = false + num_initial_low_order_iterations::Int = 3 + inverse_ncp_hessian_bound::Float64 = 1.e-2 + max_al_iters::Int = 100 + tol::Float64 = 1.e-8 +end + +struct AugmentedLagrangeSolver + settings::AugmentedLagrangeSolverSettings +end + +function AugmentedLagrangeSolver(objective_cache, p; kwargs...) + settings = AugmentedLagrangeSolverSettings(; kwargs...) + return AugmentedLagrangeSolver(settings) +end + +function linear_update() + +end + +function solve!(solver::AugmentedLagrangeSolver, Uu, p, λ, κ) + # TODO add warm start + + # TODO add update preconditioner + + huge_val = 1.e64 + ncp_error = huge_val * ones(length(λ)) + error_norm = huge_val + + initial_tol_scaling = 100.0 + tol_ramp_down_iters = 3 + + for it in 1:solver.settings.max_al_iters + # TODO finish this up... we'll need a way + # to dynamically change equation solver settings + # if it < tol_ramp_down_iters + # tol_scaling = initial_tol_scaling^((tol_ramp_down_iters - it) / tol_ramp_down_iters) + + # end + + update_precond = false + if (solver.settings.use_second_order_update && it >= solver.settings.num_initial_low_order_iterations) || + solver.settings.use_newton_only + + # do linear update + # TODO finish this + dx, dl, solve_success = linear_update(objective_cache, Uu, p, λ, κ) + + if !solve_success + update_precond = true + end + + λ_save = copy(λ) + + # TODO move to settings + max_line_search_iters = 10 + for linesearch in 1:max_line_search_iters + if linesearch == max_line_search_iters + update_precond = true + end + + y = Uu + dx + λ .+= dl + + # TODO need total residual call here + trial_error_norm = norm() + + if norm(trial_error_norm < error_norm) + @info "Total error after 2nd order update = $trial_error_norm" + error_norm = trial_error_norm + copyto!(Uu, y) + else + @info "Total error after 2nd order update = $trial_error_norm, no improvement" + copyto!(λ, λ_save) + dx .*= 0.2 + dl .*= 0.2 + end + end + end + + if update_precond + # TODO make call to update precond + end + + if !solver.settings.use_newton_only + + end + end +end diff --git a/src/solvers/ExplicitSolver.jl b/src/solvers/ExplicitSolver.jl new file mode 100644 index 0000000..5d23eea --- /dev/null +++ b/src/solvers/ExplicitSolver.jl @@ -0,0 +1,91 @@ +Base.@kwdef struct ExplicitSolverSettings + verbose::Bool = true +end + +struct ExplicitSolver{O} + objective_cache::O + settings::ExplicitSolverSettings +end + +function ExplicitSolver(objective_cache, p) + return ExplicitSolver(objective_cache, ExplicitSolverSettings()) +end + +function solve!(solver::ExplicitSolver, u, p) + Δt = p.times.Δt[1] + o = solver.objective_cache + v_old = o.solution_rate_old + v = o.solution_rate + a = o.solution_rate_rate + M = o.lumped_mass + + internal_energy = _internal_energy(o, u, p) + kinetic_energy = 0.5 * mapreduce(i -> M[i] * v[i]^2, +, eachindex(v)) + # TODO do external force (u dot f_ext is the energy) + fill!(o.internal_energy, internal_energy) + fill!(o.kinetic_energy, kinetic_energy) + + f_i = _internal_force(o, u, p) + # TODO external forces + a.data .= -f_i.data ./ M + @. v = v_old + Δt * a + @. u += Δt * v + + copyto!(v_old, v) + return nothing +end + + +# function predict!(solver::ExplicitSolver, u, p) +# # unpack a bunch of stuff +# Δt = p.times.Δt[1] +# γ = solver.objective_cache.γ +# v = solver.objective_cache.solution_rate +# a = solver.objective_cache.solution_rate_rate + +# u .= u .+= Δt * v .+ 0.5 * Δt * Δt * a +# v .= v .+= (1. - γ) * Δt * a + +# # ensure bcs are correct +# FiniteElementContainers.update_field_dirichlet_bcs!(u_pre, v_pre, p.dirichlet_bcs) +# return nothing +# end + +# function _update_acceleration(solver::ExplicitSolver, u, p) +# # solver.objective_cache.inertial_force .= +# # unpack stuff +# o = solver.objective_cache +# v = o.solution_rate +# a = o.solution_rate_rate +# M = o.lumped_mass + +# internal_energy = _internal_energy(o, u, p) +# kinetic_energy = 0.5 * dot(M, v .* v) +# energy = internal_energy + kinetic_energy +# # TODO do external force (u dot f_ext is the energy) +# fill!(o.internal_energy, internal_energy) +# fill!(o.kinetic_energy, kinetic_energy) + +# f_i = _internal_force(o, u, p) +# f_k = M .* a +# # TODO need external force +# o.gradient .= f_i .+ f_k + + +# end + +# function correct!(solver::ExplicitSolver, u, p) +# # unpack a bunch of stuff +# Δt = p.times.Δt[1] +# γ = solver.objective_cache.γ +# v = solver.objective_cache.solution_rate +# a = solver.objective_cache.solution_rate_rate + +# # need to "solve" for acceleration + +# # a .= (u - u_pre) / β / Δt / Δt +# v .+= γ * Δt * a + +# # ensure bcs are correct +# FiniteElementContainers.update_field_dirichlet_bcs!(u, v, a, p.dirichlet_bcs) +# end diff --git a/src/solvers/NewtonSolver.jl b/src/solvers/NewtonSolver.jl index 02f7829..7991e5d 100644 --- a/src/solvers/NewtonSolver.jl +++ b/src/solvers/NewtonSolver.jl @@ -1,7 +1,7 @@ Base.@kwdef struct NewtonSolverSettings max_iters::Int = 50 - rel_tol::Float64 = 1e-5 - abs_tol::Float64 = 1e-6 + rel_tol::Float64 = 1e-10 + abs_tol::Float64 = 1e-10 verbose::Bool = true end @@ -16,13 +16,7 @@ end function solve!(solver::NewtonSolver, Uu, p) objective_cache = solver.objective_cache - # R = gradient(objective_cache, Uu, p) - # res_norm0 = norm(R) - - # if res_norm0 == 0.0 - # res_norm0 = 1. - # end - res_norm0 = Inf + res_norm0 = one(eltype(Uu)) n = 1 if solver.settings.verbose @@ -35,7 +29,10 @@ function solve!(solver::NewtonSolver, Uu, p) R = gradient(objective_cache, Uu, p) if n == 1 - res_norm0 = norm(R) + norm_R = norm(R) + if norm_R != 0. + res_norm0 = norm(R) + end end K = hessian(objective_cache, Uu, p) diff --git a/src/solvers/Preconditioners.jl b/src/solvers/Preconditioners.jl index 95fce06..9f99a8a 100644 --- a/src/solvers/Preconditioners.jl +++ b/src/solvers/Preconditioners.jl @@ -56,9 +56,9 @@ function CholeskyPreconditioner(obj::AbstractObjectiveCache, p, timer) # TODO need stuff here # inefficiency here by creating these copies Uu = create_unknowns(asm) - # H = hessian(obj, Uu, p) + H = hessian(obj, Uu, p) # P = cholesky(H) - assemble_stiffness!(asm, obj.objective.hessian_u, Uu, p) + # assemble_stiffness!(asm, obj.objective.hessian_u, Uu, p) # NOTE: # note assembling since the stiffness is assembled @@ -99,8 +99,9 @@ function update_preconditioner!(P::CholeskyPreconditioner, obj, Uu, p; verbose=f # H = hessian!(P.assembler, obj, Uu, p) # H = hessian(obj, Uu, p) asm = assembler(obj) - assemble_stiffness!(asm, obj.objective.hessian_u, Uu, p) - H = stiffness(asm) + # assemble_stiffness!(asm, obj.objective.hessian_u, Uu, p) + # H = stiffness(asm) + H = hessian(obj, Uu, p; symmetric = false) attempt = 1 while attempt < 10 if verbose diff --git a/src/solvers/Solvers.jl b/src/solvers/Solvers.jl index 40535cc..c8db541 100644 --- a/src/solvers/Solvers.jl +++ b/src/solvers/Solvers.jl @@ -2,5 +2,8 @@ include("Preconditioners.jl") include("WarmStart.jl") include("EigenSolver.jl") +include("ExplicitSolver.jl") include("NewtonSolver.jl") include("TrustRegionSolver.jl") + +include("AugmentedLagrange.jl") diff --git a/src/solvers/TrustRegionSolver.jl b/src/solvers/TrustRegionSolver.jl index 106dfb4..d5b7885 100644 --- a/src/solvers/TrustRegionSolver.jl +++ b/src/solvers/TrustRegionSolver.jl @@ -426,6 +426,8 @@ function solve!(solver::TrustRegionSolver, Uu, p) cumulative_cg_iters = 0 print_banner(solver) + tried_new_precond = false + # begin loop over trust region iterations for i in 1:settings.max_trust_iters if i > 1 && i % 25 == 0 @@ -443,8 +445,11 @@ function solve!(solver::TrustRegionSolver, Uu, p) hess_vec_func = v -> hvp(solver.objective_cache, Uu, v, p) # TODO need to fix below - # mult_by_approx_hessian = v -> hvp(solver.objective_cache, Uu, p, v) + # mult_by_approx_hessian = v -> hvp(solver.objective_cache, Uu, v, p) mult_by_approx_hessian = v -> v + # mult_by_approx_hessian = v -> P.preconditioner.L * (P.preconditioner.L' * v) + # mult_by_approx_hessian = v -> P.preconditioner \ v + # mult_by_approx_hessian = v -> hessian(solver.objective_cache, Uu, p) * v # calculate cauchy point cp_norm_sq = calculate!( @@ -456,13 +461,13 @@ function solve!(solver::TrustRegionSolver, Uu, p) # minimize model problem if cp_norm_sq >= tr_size * tr_size - @warn "unpreconditioned gradient cauchy point outside trust region at dist = $(sqrt.(cauchyPointNormSquared))" + @warn "unpreconditioned gradient cauchy point outside trust region at dist = $(sqrt.(cp_norm_sq))" # cauchyPoint *= (tr_size / sqrt(cp_norm_sq)) copyto!(solver.cauchy_point.cauchy_point, (tr_size / sqrt(cp_norm_sq)) * solver.cauchy_point.cauchy_point) cp_norm_sq = tr_size * tr_size q_newton_point = solver.cauchy_point.cauchy_point - stepType = boundaryString - cgIters = 1 + step_type = boundaryString + cg_iters = 1 else # @assert false "Implement model problem" q_newton_point, _, step_type, cg_iters = @@ -498,7 +503,7 @@ function solve!(solver::TrustRegionSolver, Uu, p) rho = real_improve / model_improve - if model_objective > 0 + if model_objective > 0 && model_objective > eps(typeof(model_objective)) @warn "Found a positive model objective increase. Debug if you see this." @warn "modelObjective = $model_objective" @warn "realObjective = $real_objective" @@ -516,7 +521,7 @@ function solve!(solver::TrustRegionSolver, Uu, p) real_res_norm = norm(gy) will_accept = rho >= settings.η_1 || - (rho >= -0 && realResNorm <= g_norm) + (rho >= -0 && real_res_norm <= g_norm) print_min_banner( solver, real_objective, model_objective, @@ -559,8 +564,9 @@ function solve!(solver::TrustRegionSolver, Uu, p) if !tried_new_precond @info "The trust region is too small, updating precond and trying again." - update_preconditioner!(solver.preconditioner, solver.objective, Uu, p) - P = solver.preconditioner.preconditioner + update_preconditioner!(solver.preconditioner, solver.objective_cache, Uu, p) + # P = solver.preconditioner.preconditioner + P = solver.preconditioner cumulative_cg_iters = 0 tried_new_precond = true happy_about_tr_size = true diff --git a/src/solvers/WarmStart.jl b/src/solvers/WarmStart.jl index 3dabc6e..bc3af68 100644 --- a/src/solvers/WarmStart.jl +++ b/src/solvers/WarmStart.jl @@ -85,6 +85,7 @@ function solve!( Forward, assemble_vector!, Duplicated(R, dR), + Const(asm.vector_pattern), Const(asm.dof), Const(residual), Duplicated(Uu, dUu), diff --git a/test/TestContact.jl b/test/TestContact.jl index 995d2e2..9c7f66b 100644 --- a/test/TestContact.jl +++ b/test/TestContact.jl @@ -148,6 +148,147 @@ function test_neareast_neighbor_search_two_block_tri3_mesh(mn) end end +function test_integral_in_overlap() + penalty_length = 0.1 + edge_smoothing = 0.2 + ξ = SVector{2, Float64}(0.2, 0.9) + δ = 0.5 + @test Cthonios.integrate_gap(ξ, SVector{2, Float64}(-0.1, 0.1), δ) ≈ Inf + @test Cthonios.integrate_gap(ξ, SVector{2, Float64}(0.1, -0.1), δ) ≈ Inf + @test Cthonios.integrate_gap(ξ, SVector{2, Float64}(-0.1, -0.1), δ) ≈ Inf + @test Cthonios.integrate_gap(ξ, SVector{2, Float64}(-0.1, -0.2), δ) ≈ Inf + @test Cthonios.integrate_gap(ξ, SVector{2, Float64}(-0.1, -0.1 + 1e-12), δ) ≈ Inf + @test Cthonios.integrate_gap(ξ, SVector{2, Float64}(-0.2, 0.6), δ) ≈ Inf + @test Cthonios.integrate_gap(ξ, SVector{2, Float64}(0.7, -0.3), δ) ≈ Inf + @test Cthonios.integrate_gap(ξ, SVector{2, Float64}(0.7, 0.0), δ) ≈ Inf + @test Cthonios.integrate_gap(ξ, SVector{2, Float64}(0.0, 0.7), δ) ≈ Inf + @test Cthonios.integrate_gap(ξ, SVector{2, Float64}(0.2, 0.0), δ) ≈ Inf + @test Cthonios.integrate_gap(ξ, SVector{2, Float64}(0.0, 0.3), δ) ≈ Inf + @test Cthonios.integrate_gap(ξ, SVector{2, Float64}(0.0, 0.0), δ) ≈ Inf +end + +# function _integrate_gap_numeric(xi, g, delta) +# N = 10_000 +# xig = range(0.5/N, stop = 1.0 - 0.5/N, length = N) + +# dxi = xi[2] - xi[1] +# w = dxi / N + +# gap(x) = g[1] + x * (g[2] - g[1]) + +# p(x) = begin +# v = gap(x) +# v < delta ? (v / delta + delta / v - 2) : 0.0 +# end + +# # return sum(p, xig) * w +# return sum(map(p, xig)) * w +# end + +function _integrate_gap_numeric(xi, g, delta) + N = 10_000 + + ξ0, ξ1 = xi + g0, g1 = g + + # enforce positive orientation (important!) + if ξ1 < ξ0 + ξ0, ξ1 = ξ1, ξ0 + g0, g1 = g1, g0 + end + + dξ = ξ1 - ξ0 + w = dξ / N + + # midpoint rule on [0,1] + xg = range(0.5/N, stop = 1.0 - 0.5/N, length = N) + + function gap_at_x(x) + ξ = ξ0 + x * dξ + return g0 + (ξ - ξ0) / dξ * (g1 - g0) + end + + function p(x) + v = gap_at_x(x) + if v <= 0 + return Inf + elseif v < delta + return v / delta + delta / v - 2 + else + return 0.0 + end + end + + return sum(p(x) for x in xg) * w +end + +function test_integral_both_in_contact() + ξ = SVector{2, Float64}(0.2, 0.9) + g1 = SVector{2, Float64}(0.1, 0.3) + g2 = SVector{2, Float64}(0.3, 0.1) + δ = 0.5 + integral_1 = Cthonios.integrate_gap(ξ, g1, δ) + integral_2 = Cthonios.integrate_gap(ξ, g2, δ) + @test integral_1 ≈ integral_2 + + integral_3 = _integrate_gap_numeric(ξ, g1, δ) + integral_4 = _integrate_gap_numeric(ξ, g2, δ) + + @test isapprox(integral_1, integral_3, rtol = 1e-7) + @test isapprox(integral_2, integral_4, rtol = 1e-7) +end + +function test_integral_both_equal_in_contact() + ξ = SVector{2, Float64}(0.2, 0.9) + g = SVector{2, Float64}(0.3, 0.3) + δ = 0.5 + integral_1 = _integrate_gap_numeric(ξ, g, δ) + integral_2 = Cthonios.integrate_gap(ξ, g, δ) + @test integral_1 ≈ integral_2 + # g[2] = g[2] + 1e-8 + g = SVector{2, Float64}(0.3, 0.3 + 1e-8) + integral_3 = Cthonios.integrate_gap(ξ, g, δ) + @test isapprox(integral_2, integral_3, atol=1e-7) +end + +function test_integral_both_out_of_contact() + ξ = SVector{2, Float64}(0.2, 0.9) + g1 = SVector{2, Float64}(0.6, 0.8) + g2 = SVector{2, Float64}(0.8, 0.6) + δ = 0.5 + integral_1 = Cthonios.integrate_gap(ξ, g1, δ) + integral_2 = Cthonios.integrate_gap(ξ, g2, δ) + @test integral_1 ≈ 0.0 + @test integral_2 ≈ 0.0 + integral_3 = _integrate_gap_numeric(ξ, g1, δ) + @test isapprox(integral_1, integral_3, atol=1e-7) +end + +function test_integral_both_equal_out_of_contact() + ξ = SVector{2, Float64}(0.2, 0.9) + g1 = SVector{2, Float64}(0.8, 0.8) + g2 = SVector{2, Float64}(0.8, 0.8) + δ = 0.5 + integral_1 = Cthonios.integrate_gap(ξ, g1, δ) + integral_2 = Cthonios.integrate_gap(ξ, g2, δ) + @test integral_1 ≈ 0.0 + @test integral_2 ≈ 0.0 + integral_3 = _integrate_gap_numeric(ξ, g1, δ) + @test isapprox(integral_1, integral_3, atol = 1e-7) +end + +function test_integral_in_out_of_contact() + ξ = SVector{2, Float64}(0.2, 0.9) + g1 = SVector{2, Float64}(0.1, 0.7) + g2 = SVector{2, Float64}(0.7, 0.1) + δ = 0.5 + integral_1 = Cthonios.integrate_gap(ξ, g1, δ) + integral_2 = Cthonios.integrate_gap(ξ, g2, δ) + @test integral_1 ≈ integral_2 + integral_3 = Cthonios.integrate_gap(ξ, g1, δ) + @test isapprox(integral_1, integral_3, atol = 1e-7) +end + @testset "ContactPair" begin test_contact_pair_constructor() end @@ -163,3 +304,12 @@ end test_neareast_neighbor_search_two_block_tri3_mesh(mn) end end + +@testset "Integral penalty contact tests" begin + test_integral_in_overlap() + test_integral_both_in_contact() + test_integral_both_equal_in_contact() + test_integral_both_out_of_contact() + test_integral_both_equal_out_of_contact() + test_integral_in_out_of_contact() +end diff --git a/test/TestPatchTests.jl b/test/TestPatchTests.jl index 1ebc630..6c5663d 100644 --- a/test/TestPatchTests.jl +++ b/test/TestPatchTests.jl @@ -8,6 +8,7 @@ function linear_patch_test_dirichlet(mesh_file, q_degree) ) props = (; var"" = Dict{String, Any}( + "density" => 1.0, "bulk modulus" => 1.0, "shear modulus" => 0.25 ) @@ -20,14 +21,14 @@ function linear_patch_test_dirichlet(mesh_file, q_degree) func_2(x, t) = 0.5 * t dirichlet_bcs = [ - DirichletBC("displ_x", "sset_1", func_1) - DirichletBC("displ_x", "sset_2", func_1) - DirichletBC("displ_x", "sset_3", func_1) - DirichletBC("displ_x", "sset_4", func_1) - DirichletBC("displ_y", "sset_1", func_1) - DirichletBC("displ_y", "sset_2", func_1) - DirichletBC("displ_y", "sset_3", func_1) - DirichletBC("displ_y", "sset_4", func_1) + DirichletBC("displ_x", func_1; sideset_name = "sset_1") + DirichletBC("displ_x", func_1; sideset_name = "sset_2") + DirichletBC("displ_x", func_1; sideset_name = "sset_3") + DirichletBC("displ_x", func_1; sideset_name = "sset_4") + DirichletBC("displ_y", func_1; sideset_name = "sset_1") + DirichletBC("displ_y", func_1; sideset_name = "sset_2") + DirichletBC("displ_y", func_1; sideset_name = "sset_3") + DirichletBC("displ_y", func_1; sideset_name = "sset_4") ] sim = SingleDomainSimulation( diff --git a/test/TestPhysics.jl b/test/TestPhysics.jl index c627378..ce95336 100644 --- a/test/TestPhysics.jl +++ b/test/TestPhysics.jl @@ -1,7 +1,8 @@ function neohookean_props() props = Dict( - "bulk modulus" => 1000., - "shear modulus" => 1. + "density" => 1.0, + "bulk modulus" => 1000., + "shear modulus" => 1.0 ) return props end @@ -13,8 +14,9 @@ function test_sm_material_model_input_neohookean() physics = SolidMechanics(formulation, model) props = FiniteElementContainers.create_properties(physics, props) - @test props[1] ≈ 1000. - @test props[2] ≈ 1. + @test props[1] ≈ 1.0 + @test props[2] ≈ 1000. + @test props[3] ≈ 1. state = FiniteElementContainers.create_initial_state(physics) @test state == SVector{0, Float64}() diff --git a/test/TestSimulations.jl b/test/TestSimulations.jl index 36146b7..83cdb9b 100644 --- a/test/TestSimulations.jl +++ b/test/TestSimulations.jl @@ -11,8 +11,9 @@ function test_single_domain_simulation() ) props = (; block_1 = Dict{String, Any}( - "bulk modulus" => 10.0, - "shear modulus" => 1.0 + "density" => 1.0, + "bulk modulus" => 10.0, + "shear modulus" => 1.0 ) ) @@ -20,12 +21,12 @@ function test_single_domain_simulation() func_2(x, t) = -0.5 * t dirichlet_bcs = [ - DirichletBC("displ_x", "sset_y_negative", func_1) - DirichletBC("displ_y", "sset_y_negative", func_1) - DirichletBC("displ_z", "sset_y_negative", func_1) - DirichletBC("displ_x", "sset_y_positive", func_1) - DirichletBC("displ_z", "sset_y_positive", func_1) - DirichletBC("displ_y", "sset_y_positive", func_2) + DirichletBC("displ_x", func_1; sideset_name = "sset_y_negative") + DirichletBC("displ_y", func_1; sideset_name = "sset_y_negative") + DirichletBC("displ_z", func_1; sideset_name = "sset_y_negative") + DirichletBC("displ_x", func_1; sideset_name = "sset_y_positive") + DirichletBC("displ_z", func_1; sideset_name = "sset_y_positive") + DirichletBC("displ_y", func_2; sideset_name = "sset_y_positive") ] sim = SingleDomainSimulation( @@ -38,7 +39,7 @@ function test_single_domain_simulation() @test sim.output_file == output_file @test sim.times == times @test sim.physics == physics - @test sim.properties == (; block_1 = [10., 1.]) + @test sim.properties == (; block_1 = [1., 10., 1.]) # TODO add more testing on ics, bcs, and contact pairs end diff --git a/test/objectives/TestDesignObjective.jl b/test/objectives/TestDesignObjective.jl index dceaf16..b7035bb 100644 --- a/test/objectives/TestDesignObjective.jl +++ b/test/objectives/TestDesignObjective.jl @@ -1,74 +1,154 @@ +using Cthonios +using ConstitutiveModels +using Distributions +using FiniteElementContainers +using LinearAlgebra +using Random +using Test + Random.seed!(123) -function build_direction_vector(n_design_vars) +# function build_direction_vector(n_design_vars) +function build_direction_vector(params) + if typeof(params) <: H1Field + n_design_vars = length(params.data) + elseif typeof(params) <: NamedTuple + n_design_vars = mapreduce(length, +, values(params)) + else + @assert false + end dir_vec = rand(Uniform(-1., 1.), n_design_vars) return dir_vec / norm(dir_vec) end +function generic_dot(a, b) + if typeof(b) <: H1Field + return dot(a, b.data) + elseif typeof(b) <: NamedTuple + # NOTE this will likely fail for multi-block problems + # b_vec = reduce(vcat, values(b)) + # return dot(a, b_vec) + return dot(a, values(b)[1]) + else + @assert false + end +end + +function perturb_params(params, step_size, direction_vec) + if typeof(params) <: H1Field + return H1Field(params + step_size * reshape(direction_vec, size(params))) + elseif typeof(params) <: NamedTuple + # return map( + # (x, y) -> x .+ step_size * y, + # values(params), direction_vec + # ) + # NOTE this will likely fail for multi-block problems + # new_params = deepcopy(params) + # for val in values(new_params) + # val .=+ step_size * direction_vec + # end + new_params = NamedTuple{keys(params)}(( + values(params)[1] + step_size * direction_vec, + )) + return new_params + else + @assert false + end +end + # TODO this is hardcoded for x Sensitivity function finite_difference_error( - forward, - objective, sim, qoi, - step_size + design_obj, physics_obj, sim, + step_size, param_sym ) - objective_cache, U, p = setup_caches(objective, sim; q_degree=1) - # sens = forward(objective_cache, U, p, qoi, 2) - des_obj = forward(objective_cache, U, p, qoi, 2) - X = p.h1_coords + objective_cache, U, p = setup_caches(physics_obj, sim; q_degree=1) + solver = TrustRegionSolver(objective_cache, p; verbose=false) + des_obj = design_obj(objective_cache, U, p) + Cthonios.forward_problem!(des_obj, solver, objective_cache, U, p) + + params = getfield(p, param_sym) - # og_obj, og_grad = Cthonios.gradient_x_and_value(sens, U, p) - # og_obj, og_grad = Cthonios.gradient_and_value(des_obj, U, p) - og_obj, og_grad = Cthonios.gradient_and_value(des_obj, X) + og_obj, og_grad = Cthonios.gradient_and_value(des_obj, params) - direction_vec = build_direction_vector(length(X.data)) - directional_deriv = dot(direction_vec, og_grad.data) + direction_vec = build_direction_vector(params) + directional_deriv = generic_dot(direction_vec, og_grad) + params_perturbed = perturb_params(params, step_size, direction_vec) - X_perturbed = H1Field(X + step_size * reshape(direction_vec, size(X))) + objective_cache, U, p = setup_caches(physics_obj, sim; q_degree=1) + solver = TrustRegionSolver(objective_cache, p; verbose=false) + des_obj = design_obj(objective_cache, U, p) + Cthonios.forward_problem!(des_obj, solver, objective_cache, U, p) - objective_cache, U, p = setup_caches(objective, sim; q_degree=1) - copyto!(p.h1_coords, X_perturbed) + obj_perturbed = Cthonios.value(des_obj, params_perturbed) - # sens = forward(objective_cache, U, p, qoi, 2) - des_obj = forward(objective_cache, U, p, qoi, 2) - # obj_perturbed = Cthonios.value(sens, U, p) - # obj_perturbed = Cthonios.value(des_obj, U, p) - obj_perturbed = Cthonios.value(des_obj, X_perturbed) + display(og_obj) + display(og_grad) + display(obj_perturbed) fd_value = (obj_perturbed - og_obj) / step_size error = abs(directional_deriv - fd_value) return error end -function strain_energy_objective_with_coords!(f, qois, Us, ps, params) +function strain_energy_objective_with_coords!(f, qoi_storages, qois, Us, ps, params) U, p = Us[end], ps[end] copyto!(p.h1_coords.data, params.data) - Cthonios._value!(f, qois[1], U, p) + Cthonios.value!(f, qoi_storages[end], qois[1], U, p) return nothing end -function neohookean_forward_problem!(objective_cache, U, p, qoi, n_steps) - # sens = Cthonios.Sensitivity(qoi(objective_cache), U, p) - solver = TrustRegionSolver(objective_cache, p) - Cthonios.initialize!(objective_cache, U, p) +function strain_energy_objective_with_props!(f, qoi_storages, qois, Us, ps, params) + U, p = Us[end], ps[end] + for (n, val) in enumerate(values(params)) + copyto!(values(p.properties)[n], val) + end + Cthonios.value!(f, qoi_storages[end], qois[1], U, p) + return nothing +end - design_obj = Cthonios.DesignObjective( - strain_energy_objective_with_coords!, - [qoi(objective_cache)], U, p - ) +function total_work_objective_with_coords!(f, qoi_storages, qois, Us, ps, params) + # getting the y bc dofs on sset_3 + dofs = ps[end].dirichlet_bcs.bc_caches[4].dofs - for _ in 1:n_steps - Cthonios.step!(solver, objective_cache, U, p) - push!(design_obj.stored_solutions, deepcopy(U)) - push!(design_obj.stored_parameters, deepcopy(p)) + # update qoi storages + for n in 1:length(Us) + copyto!(ps[n].h1_coords.data, params.data) + Cthonios._update_storage!(qoi_storages[n], qois[1], Us[n], ps[n]) end - return design_obj + val = 0. + for n in 2:length(Us) + U, U_prev = Us[n], Us[n - 1] + force, force_prev = qoi_storages[n], qoi_storages[n - 1] + + # u, u_prev = sum(U.data[dofs]) / length(dofs), sum(U_prev.data[dofs]) / length(dofs) + # f_r, f_r_prev = sum(force.data[dofs]), sum(force_prev.data[dofs]) + + # val += 0.5 * (f_r + f_r_prev) * (u - u_prev) + + # @show dot( + # (force.data[dofs] + force_prev.data[dofs]), + # (U.data[dofs] - U_prev.data[dofs]) + # ) + val += @views 0.5 * dot( + (force.data[dofs] + force_prev.data[dofs]), + (U.data[dofs] - U_prev.data[dofs]) + ) + # val += @views 0.5 * dot( + # (force + force_prev), + # (U - U_prev) + # ) + # val += 0.5 * dot(temp1, temp2) + # val += 0.5 * dot(temp1, temp2) + end + fill!(f, val) + return nothing end function neohookean_sim_helper() mesh_file = dirname(Base.source_dir()) * "/mesh/patch_test_mesh_quad4.g" output_file = splitext(mesh_file)[1] * "-output.exo" - times = TimeStepper(0., 1., 2) + times = TimeStepper(0., 1., 10) physics = (; var"" = SolidMechanics( # Block1 = SolidMechanics( @@ -87,10 +167,14 @@ function neohookean_sim_helper() func_2(x, t) = -0.05 * t dirichlet_bcs = [ - DirichletBC("displ_x", "sset_1", func_1) - DirichletBC("displ_y", "sset_1", func_1) - DirichletBC("displ_x", "sset_3", func_1) - DirichletBC("displ_y", "sset_3", func_2) + # DirichletBC("displ_x", "sset_1", func_1) + # DirichletBC("displ_y", "sset_1", func_1) + # DirichletBC("displ_x", "sset_3", func_1) + # DirichletBC("displ_y", "sset_3", func_2) + DirichletBC("displ_x", func_1; sideset_name = "sset_3") + DirichletBC("displ_y", func_1; sideset_name = "sset_3") + DirichletBC("displ_x", func_1; sideset_name = "sset_1") + DirichletBC("displ_y", func_2; sideset_name = "sset_1") ] objective = Cthonios.QuasiStaticObjective() sim = SingleDomainSimulation( @@ -102,32 +186,107 @@ function neohookean_sim_helper() end function neohookean_self_adjoint() - objective, sim = neohookean_sim_helper() - qoi = x -> Cthonios.QOIExtractor( + physics_obj, sim = neohookean_sim_helper() + qoi = x -> QOIExtractor( x, helmholtz_free_energy, +, FiniteElementContainers.L2QuadratureField, Float64; + is_material_qoi = true, reduction_2 = + ) + design_obj_1 = (x, u, p) -> Cthonios.DesignObjective( + strain_energy_objective_with_coords!, + [qoi(x)], u, p + ) + design_obj_2 = (x, u, p) -> Cthonios.DesignObjective( + strain_energy_objective_with_props!, + [qoi(x)], u, p + ) - - fd_err = finite_difference_error( - neohookean_forward_problem!, - objective, sim, qoi, 1e-7 + for (design_obj, params_sym) in zip( + [design_obj_1, design_obj_2], + [:h1_coords, :properties] ) - @test fd_err <= 1e-7 - - # # TODO setup test for checking v shape - # step_sizes = [1e-4, 1e-5, 1e-6, 1e-7, 1e-8, 1e-9, 1e-10] - # fd_errs = Float64[] - # for step_size in step_sizes - # fd_err = finite_difference_error( - # neohookean_forward_problem!, - # objective, sim, qoi, step_size - # ) - # push!(fd_errs, fd_err) - # end - - # display(fd_errs) + fd_err = finite_difference_error( + design_obj, physics_obj, sim, 1e-7, + params_sym + ) + @test fd_err <= 1e-7 + + # TODO setup test for checking v shape + step_sizes = [1e-4, 1e-5, 1e-6, 1e-7, 1e-8, 1e-9, 1e-10] + fd_errs = Float64[] + for step_size in step_sizes + fd_err = finite_difference_error( + design_obj, physics_obj, sim, step_size, + params_sym + ) + push!(fd_errs, fd_err) + end + display(fd_errs) + end +end + +# TODO this test should be failing +# we might have something wrong in our derivative shadowing +# since we're recycling the storage in the qois +function neohookean_non_self_adjoint() + physics_obj, sim = neohookean_sim_helper() + qoi = x -> QOIExtractor( + x, residual, identity, + H1Field, Float64; + is_field_qoi = true, + reduction_2 = identity + ) + qoi_energy = x -> QOIExtractor( + x, helmholtz_free_energy, +, + FiniteElementContainers.L2QuadratureField, Float64; + is_material_qoi = true, + reduction_2 = + + ) + design_obj_1 = (x, u, p) -> Cthonios.DesignObjective( + total_work_objective_with_coords!, + [qoi(x)], u, p + ) + + for (design_obj, params_sym) in zip( + # [design_obj_1, design_obj_2], + # [:h1_coords, :properties] + [design_obj_1], + [:h1_coords] + ) + fd_err = finite_difference_error( + design_obj, physics_obj, sim, 1e-7, + params_sym + ) + @show fd_err + @test_throws Exception fd_err <= 1e-7 + + # TODO setup test for checking v shape + # step_sizes = [1e-4, 1e-5, 1e-6, 1e-7, 1e-8, 1e-9, 1e-10] + # fd_errs = Float64[] + # for step_size in step_sizes + # fd_err = finite_difference_error( + # design_obj, physics_obj, sim, step_size, + # params_sym + # ) + # push!(fd_errs, fd_err) + # end + # display(fd_errs) + end + + # objective_cache, U, p = setup_caches(physics_obj, sim; q_degree=1) + # solver = TrustRegionSolver(objective_cache, p; verbose=false) + # des_obj = design_obj_1(objective_cache, U, p) + # Cthonios.forward_problem!(des_obj, solver, objective_cache, U, p) + + # f = zeros(1) + # params = p.h1_coords + # Cthonios.value!(f, des_obj, params) + # display(f) + + # qoi_temp = qoi_energy(objective_cache) + # Cthonios.value(qoi_temp, U, p) end neohookean_self_adjoint() +neohookean_non_self_adjoint() diff --git a/test/runtests.jl b/test/runtests.jl index 73ea22d..525a106 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -22,6 +22,7 @@ function sim_helper() ) props = (; block_1 = Dict{String, Any}( + "density" => 1.0, "bulk modulus" => 10.0, "shear modulus" => 1.0 ) @@ -31,12 +32,12 @@ function sim_helper() func_2(x, t) = -0.5 * t dirichlet_bcs = [ - DirichletBC("displ_x", "sset_y_negative", func_1) - DirichletBC("displ_y", "sset_y_negative", func_1) - DirichletBC("displ_z", "sset_y_negative", func_1) - DirichletBC("displ_x", "sset_y_positive", func_1) - DirichletBC("displ_z", "sset_y_positive", func_1) - DirichletBC("displ_y", "sset_y_positive", func_2) + DirichletBC("displ_x", func_1; sideset_name = "sset_y_negative") + DirichletBC("displ_y", func_1; sideset_name = "sset_y_negative") + DirichletBC("displ_z", func_1; sideset_name = "sset_y_negative") + DirichletBC("displ_x", func_1; sideset_name = "sset_y_positive") + DirichletBC("displ_z", func_1; sideset_name = "sset_y_positive") + DirichletBC("displ_y", func_2; sideset_name = "sset_y_positive") ] sim = SingleDomainSimulation( @@ -56,7 +57,7 @@ end end @testset "Objectives" begin - include("objectives/TestDesignObjective.jl") + # include("objectives/TestDesignObjective.jl") include("objectives/TestObjectives.jl") end