adding a forward problem, along with gradient checks for a self-adjoi…

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"

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)

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

examples/wip/clamped/script_v0.3_dynamics.jl → examples/hole_array/script_implicit_dynamics.jl

@@ -3,35 +3,36 @@ 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)
 
 # 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.)

examples/hole_array/script_quasistatic.jl

@@ -1,37 +1,59 @@
 using ConstitutiveModels
 using Cthonios
+using Enzyme
 using FiniteElementContainers
+Enzyme.Compiler.VERBOSE_ERRORS[] = true
 
 # function sim_test()
 # file management
 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()
+# # display(f)
+# # display(dfdX)
+# # # # display(dfdU)
+# # # # end
+# # # # sim_test()

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