From 5ac789481058296cb4cec047a782813cb541b461 Mon Sep 17 00:00:00 2001 From: "Craig M. Hamel" Date: Sun, 9 Aug 2026 18:03:52 -0400 Subject: [PATCH 1/2] stashing some stuff. --- Project.toml | 4 +- .../hole_array/script_quasistatic.jl | 18 +++++-- src/Physics.jl | 51 +++++++++++-------- test/TestExplicitDynamics.jl | 4 +- 4 files changed, 50 insertions(+), 27 deletions(-) diff --git a/Project.toml b/Project.toml index 1bd511f..ef01a4c 100644 --- a/Project.toml +++ b/Project.toml @@ -26,10 +26,10 @@ TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" [compat] Aqua = "0.8" Arpack = "0.5" -ConstitutiveModels = "0.3" +ConstitutiveModels = "0.4" Distributions = "0.25" DocStringExtensions = "0.9" -FiniteElementContainers = "0.14" +FiniteElementContainers = "0.15" ForwardDiff = "1" KernelAbstractions = "0.9" Krylov = "0.10" diff --git a/examples/mechanics/hole_array/script_quasistatic.jl b/examples/mechanics/hole_array/script_quasistatic.jl index ba18300..aac9427 100644 --- a/examples/mechanics/hole_array/script_quasistatic.jl +++ b/examples/mechanics/hole_array/script_quasistatic.jl @@ -8,11 +8,14 @@ mesh_file = Base.source_dir() * "/mesh/hole_array.exo" output_file = splitext(mesh_file)[1] * "-output.exo" # Times -times = TimeStepper(0., 1., 20) +t_end = 100.0 +times = TimeStepper(0., 2. * t_end, 40) # Physics physics = (; - Block1 = SolidMechanics(PlaneStrain(), NeoHookean()) + # Block1 = SolidMechanics(PlaneStrain(), Hyperelastic(NeoHookean())) + Block1 = SolidMechanics(PlaneStrain(), FeFv()) + # Block1 = SolidMechanics(PlaneStrain(), Hyperelastic(Hencky())) ) props = (; Block1 = Dict{String, Any}( @@ -20,11 +23,20 @@ props = (; "Young's modulus" => 1., "Poisson's ratio" => 0.48, # "Jm" => 3. + "G neq" => 50.0, + "relaxation time" => 50.0 ) ) # Boundary Conditions -func_1(x, t) = -7.5 * t +# func_1(x, t) = -9. * (t / t_end) +function func_1(x, t) + if t >= t_end + return -9. + else + return -9. * (t / t_end) + end +end func_2(x, t) = 0.0 dirichlet_bcs = [ diff --git a/src/Physics.jl b/src/Physics.jl index e2fd66e..1909b48 100644 --- a/src/Physics.jl +++ b/src/Physics.jl @@ -25,11 +25,13 @@ 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, props_el, dt, ∇u_q, θ, state_old_q, state_new_q + physics.constitutive_model, props_el, state_old_q, state_new_q, dt, ∇u_q, θ ) σ_q = ConstitutiveModels.cauchy_stress( - physics.constitutive_model, props_el, dt, ∇u_q, θ, state_old_q, state_new_q + # physics.constitutive_model, props_el, dt, ∇u_q, θ, state_old_q, state_new_q + physics.constitutive_model, props_el, state_old_q, state_new_q, dt, ∇u_q, θ ) σ_q = σ_q |> symmetric return StandardMaterialOutput(ψ_q, σ_q, ∇u_q) @@ -193,12 +195,13 @@ end # kinematics ∇u_q = modify_field_gradients(physics.formulation, ∇u_q) - mat_props = @views props_el[2:end] + # mat_props = @views props_el[2:end] # constitutive θ = 0.0 # TODO ψ_q = ConstitutiveModels.helmholtz_free_energy( - physics.constitutive_model, mat_props, dt, ∇u_q, θ, state_old_q, state_new_q + # physics.constitutive_model, mat_props, dt, ∇u_q, θ, state_old_q, state_new_q + physics.constitutive_model, props_el, state_old_q, state_new_q, dt, ∇u_q, θ ) return JxW * ψ_q end @@ -215,12 +218,13 @@ end # kinematics ∇u_q = modify_field_gradients(physics.formulation, ∇u_q) - mat_props = @views props_el[2:end] + # mat_props = @views props_el[2:end] # constitutive θ = 0.0 # TODO ψ_q = func( - physics.constitutive_model, mat_props, dt, ∇u_q, θ, state_old_q, state_new_q + # physics.constitutive_model, mat_props, dt, ∇u_q, θ, state_old_q, state_new_q + physics.constitutive_model, props_el, state_old_q, state_new_q, dt, ∇u_q, θ ) return ψ_q end @@ -236,12 +240,13 @@ end # kinematics ∇u_q = modify_field_gradients(physics.formulation, ∇u_q) - mat_props = @views props_el[2:end] + # mat_props = @views props_el[2:end] # constitutive θ = 0.0 # TODO ψ_q = func( - physics.constitutive_model, mat_props, dt, ∇u_q, θ, state_old_q, state_new_q + # physics.constitutive_model, mat_props, dt, ∇u_q, θ, state_old_q, state_new_q + physics.constitutive_model, props_el, state_old_q, state_new_q, dt, ∇u_q, θ ) return JxW * ψ_q end @@ -346,12 +351,13 @@ end # kinematics ∇u_q = modify_field_gradients(physics.formulation, ∇u_q) - mat_props = @views props_el[2:end] + # mat_props = @views props_el[2:end] # constitutive θ = 0.0 # TODO P_q = ConstitutiveModels.pk1_stress( - physics.constitutive_model, mat_props, dt, ∇u_q, θ, state_old_q, state_new_q + # physics.constitutive_model, mat_props, dt, ∇u_q, θ, state_old_q, state_new_q + physics.constitutive_model, props_el, state_old_q, state_new_q, dt, ∇u_q, θ ) # turn into voigt notation P_q = extract_stress(physics.formulation, P_q) @@ -374,12 +380,13 @@ end # kinematics ∇u_q = modify_field_gradients(physics.formulation, ∇u_q) - mat_props = @views props_el[2:end] + # mat_props = @views props_el[2:end] # constitutive θ = 0.0 # TODO P_q = ConstitutiveModels.pk1_stress( - physics.constitutive_model, mat_props, dt, ∇u_q, θ, state_old_q, state_new_q + # physics.constitutive_model, mat_props, dt, ∇u_q, θ, state_old_q, state_new_q + physics.constitutive_model, props_el, state_old_q, state_new_q, dt, ∇u_q, θ ) scatter_with_gradients!(storage, physics.formulation, e, conn, ∇N_X, JxW * P_q) return nothing @@ -395,12 +402,13 @@ end # kinematics ∇u_q = modify_field_gradients(physics.formulation, ∇u_q) - mat_props = @views props_el[2:end] + # mat_props = @views props_el[2:end] # constitutive θ = 0. # TODO A_q = ConstitutiveModels.material_tangent( - physics.constitutive_model, mat_props, dt, ∇u_q, θ, state_old_q, state_new_q + # physics.constitutive_model, mat_props, dt, ∇u_q, θ, state_old_q, state_new_q + physics.constitutive_model, props_el, state_old_q, state_new_q, dt, ∇u_q, θ ) # turn into voigt notation K_q = extract_stiffness(physics.formulation, A_q) @@ -421,12 +429,13 @@ end # kinematics ∇u_q = modify_field_gradients(physics.formulation, ∇u_q) - mat_props = @views props_el[2:end] + # mat_props = @views props_el[2:end] # constitutive θ = 0. # TODO A_q = ConstitutiveModels.material_tangent( - physics.constitutive_model, mat_props, dt, ∇u_q, θ, state_old_q, state_new_q + # physics.constitutive_model, mat_props, dt, ∇u_q, θ, state_old_q, state_new_q + physics.constitutive_model, props_el, state_old_q, state_new_q, dt, ∇u_q, θ ) scatter_with_gradients_and_gradients!( storage, physics.formulation, e, conn, ∇N_X, JxW * A_q @@ -442,12 +451,13 @@ end ∇u_q = interpolate_field_gradients(physics, interps, u_el) ∇u_q = modify_field_gradients(physics.formulation, ∇u_q) - mat_props = @views props_el[2:end] + # mat_props = @views props_el[2:end] # constitutive θ = 0. # TODO A_q = ConstitutiveModels.material_tangent( - physics.constitutive_model, mat_props, dt, ∇u_q, θ, state_old_q, state_new_q + # physics.constitutive_model, mat_props, dt, ∇u_q, θ, state_old_q, state_new_q + physics.constitutive_model, props_el, state_old_q, state_new_q, dt, ∇u_q, θ ) # turn into voigt notation K_q = extract_stiffness(physics.formulation, A_q) @@ -466,12 +476,13 @@ end ∇u_q = interpolate_field_gradients(physics, interps, u_el) ∇u_q = modify_field_gradients(physics.formulation, ∇u_q) - mat_props = @views props_el[2:end] + # mat_props = @views props_el[2:end] # constitutive θ = 0. # TODO A_q = ConstitutiveModels.material_tangent( - physics.constitutive_model, mat_props, dt, ∇u_q, θ, state_old_q, state_new_q + # physics.constitutive_model, mat_props, dt, ∇u_q, θ, state_old_q, state_new_q + physics.constitutive_model, props_el, state_old_q, state_new_q, dt, ∇u_q, θ ) scatter_with_gradients_and_gradients!(storage, physics.formulation, e, conn, ∇N_X, JxW * A_q, v_el) return nothing diff --git a/test/TestExplicitDynamics.jl b/test/TestExplicitDynamics.jl index 6045e8f..5b2560b 100644 --- a/test/TestExplicitDynamics.jl +++ b/test/TestExplicitDynamics.jl @@ -8,12 +8,12 @@ function test_explicit_dynamics() # Physics physics = (; - cube = SolidMechanics( + sphere = SolidMechanics( ThreeDimensional(), NeoHookean() ) ) props = (; - cube = Dict{String, Any}( + sphere = Dict{String, Any}( "density" => 1.e3, "Young's modulus" => 1.e4, "Poisson's ratio" => 0.33 From be3a7259dfb4c82be91efef98c87d06a800084ce Mon Sep 17 00:00:00 2001 From: "Craig M. Hamel" Date: Fri, 14 Aug 2026 12:42:51 -0400 Subject: [PATCH 2/2] bumping fec and cm versions and makings some tweaks to stable time step calculations. --- src/Physics.jl | 94 ++++++++++----------- src/objectives/ExplicitDynamicsObjective.jl | 34 ++------ test/TestExplicitDynamics.jl | 2 +- test/TestHeatConduction.jl | 4 +- test/TestImplicitDynamics.jl | 6 +- test/TestQuasiStatic.jl | 2 +- 6 files changed, 55 insertions(+), 87 deletions(-) diff --git a/src/Physics.jl b/src/Physics.jl index 1909b48..0a73035 100644 --- a/src/Physics.jl +++ b/src/Physics.jl @@ -25,12 +25,10 @@ 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, props_el, state_old_q, state_new_q, dt, ∇u_q, θ ) σ_q = ConstitutiveModels.cauchy_stress( - # physics.constitutive_model, props_el, dt, ∇u_q, θ, state_old_q, state_new_q physics.constitutive_model, props_el, state_old_q, state_new_q, dt, ∇u_q, θ ) σ_q = σ_q |> symmetric @@ -159,28 +157,43 @@ function setup_function(fspace::FunctionSpace, ::SolidMechanics) return VectorFunction(fspace, "displ") end -function characteristic_element_length( - physics::SolidMechanics, interps, x_el, - t, dt, u_el, u_el_old, - state_old_q, state_new_q, props_el, -) +@inline function characteristic_element_length(x_el::SMatrix{2, N, T, N2}, u_el::SMatrix{2, N, T, N2}) where {N, T, N2} + x_cur = x_el + u_el + ndim = 2 + nnpe = length(x_cur) ÷ ndim + cx = cy = zero(T) + for i in 1:nnpe + cx += x_cur[(i - 1) * ndim + 1] + cy += x_cur[(i - 1) * ndim + 2] + end + cx /= nnpe; cy /= nnpe + total = zero(T) + for i in 1:nnpe + dx = x_cur[(i - 1) * ndim + 1] - cx + dy = x_cur[(i - 1) * ndim + 2] - cy + total += sqrt(dx * dx + dy * dy) + end + return 2 * total / nnpe +end + +@inline function characteristic_element_length(x_el::SMatrix{3, N, T, N3}, u_el::SMatrix{3, N, T, N3}) where {N, T, N3} x_cur = x_el + u_el ndim = 3 nnpe = length(x_cur) ÷ ndim - T = eltype(x_cur) + # T = eltype(x_cur) cx = cy = cz = zero(T) for i in 1:nnpe - cx += x_cur[(i-1)*ndim + 1] - cy += x_cur[(i-1)*ndim + 2] - cz += x_cur[(i-1)*ndim + 3] + cx += x_cur[(i - 1) * ndim + 1] + cy += x_cur[(i - 1) * ndim + 2] + cz += x_cur[(i - 1) * ndim + 3] end cx /= nnpe; cy /= nnpe; cz /= nnpe total = zero(T) for i in 1:nnpe - dx = x_cur[(i-1)*ndim + 1] - cx - dy = x_cur[(i-1)*ndim + 2] - cy - dz = x_cur[(i-1)*ndim + 3] - cz - total += sqrt(dx*dx + dy*dy + dz*dz) + dx = x_cur[(i - 1) * ndim + 1] - cx + dy = x_cur[(i - 1) * ndim + 2] - cy + dz = x_cur[(i - 1) * ndim + 3] - cz + total += sqrt(dx * dx + dy * dy + dz * dz) end return 2 * total / nnpe end @@ -200,7 +213,6 @@ end # constitutive θ = 0.0 # TODO ψ_q = ConstitutiveModels.helmholtz_free_energy( - # physics.constitutive_model, mat_props, dt, ∇u_q, θ, state_old_q, state_new_q physics.constitutive_model, props_el, state_old_q, state_new_q, dt, ∇u_q, θ ) return JxW * ψ_q @@ -218,12 +230,9 @@ 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, mat_props, dt, ∇u_q, θ, state_old_q, state_new_q physics.constitutive_model, props_el, state_old_q, state_new_q, dt, ∇u_q, θ ) return ψ_q @@ -240,12 +249,9 @@ 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, mat_props, dt, ∇u_q, θ, state_old_q, state_new_q physics.constitutive_model, props_el, state_old_q, state_new_q, dt, ∇u_q, θ ) return JxW * ψ_q @@ -294,19 +300,9 @@ end ) interps = map_interpolants(interps, x_el) (; X_q, N, ∇N_X, JxW) = interps - # ∇u_q = interpolate_field_gradients(physics, interps, u_el) - - # kinematics - # ∇u_q = modify_field_gradients(physics.formulation, ∇u_q) - # mat_props = @views props_el[2:end] ρ = props_el[1] - # constitutive - # θ = 0. # TODO - # A_q = ConstitutiveModels.material_tangent( - # physics.constitutive_model, mat_props, dt, ∇u_q, θ, state_old_q, state_new_q - # ) scatter_with_values_and_values!( storage, physics.formulation, e, conn, N, JxW * ρ ) @@ -351,17 +347,13 @@ 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, mat_props, dt, ∇u_q, θ, state_old_q, state_new_q physics.constitutive_model, props_el, state_old_q, state_new_q, dt, ∇u_q, θ ) # turn into voigt notation P_q = extract_stress(physics.formulation, P_q) - # P_q = tovoigt(SVector, P_q) G_q = discrete_gradient(physics.formulation, ∇N_X) f_q = G_q * P_q return JxW * f_q[:] @@ -380,18 +372,30 @@ 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, mat_props, dt, ∇u_q, θ, state_old_q, state_new_q physics.constitutive_model, props_el, state_old_q, state_new_q, dt, ∇u_q, θ ) scatter_with_gradients!(storage, physics.formulation, e, conn, ∇N_X, JxW * P_q) return nothing end +@inline function stable_time_step( + physics::SolidMechanics, interps, x_el, t, dt, u_el, u_el_old, state_old_q, state_new_q, props_el +) + N, D = size(interps.∇N_ξ) + x_el = SMatrix{D, N, eltype(x_el), N * D}(x_el.data) + u_el = SMatrix{D, N, eltype(u_el), N * D}(u_el.data) + l = characteristic_element_length(x_el, u_el) + ρ = props_el[1] + M = p_wave_modulus(physics.constitutive_model, props_el) + c_p = sqrt(M / ρ) + CFL = 0.1 + stable_dt = CFL * l / c_p + return stable_dt +end + @inline function FiniteElementContainers.stiffness( physics::SolidMechanics, interps, x_el, t, dt, u_el, u_el_old, state_old_q, state_new_q, props_el ) @@ -402,12 +406,9 @@ 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, mat_props, dt, ∇u_q, θ, state_old_q, state_new_q physics.constitutive_model, props_el, state_old_q, state_new_q, dt, ∇u_q, θ ) # turn into voigt notation @@ -429,12 +430,9 @@ 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, mat_props, dt, ∇u_q, θ, state_old_q, state_new_q physics.constitutive_model, props_el, state_old_q, state_new_q, dt, ∇u_q, θ ) scatter_with_gradients_and_gradients!( @@ -451,12 +449,9 @@ end ∇u_q = interpolate_field_gradients(physics, interps, u_el) ∇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, mat_props, dt, ∇u_q, θ, state_old_q, state_new_q physics.constitutive_model, props_el, state_old_q, state_new_q, dt, ∇u_q, θ ) # turn into voigt notation @@ -476,12 +471,9 @@ end ∇u_q = interpolate_field_gradients(physics, interps, u_el) ∇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, mat_props, dt, ∇u_q, θ, state_old_q, state_new_q physics.constitutive_model, props_el, state_old_q, state_new_q, dt, ∇u_q, θ ) scatter_with_gradients_and_gradients!(storage, physics.formulation, e, conn, ∇N_X, JxW * A_q, v_el) diff --git a/src/objectives/ExplicitDynamicsObjective.jl b/src/objectives/ExplicitDynamicsObjective.jl index 80d48a0..f701489 100644 --- a/src/objectives/ExplicitDynamicsObjective.jl +++ b/src/objectives/ExplicitDynamicsObjective.jl @@ -171,41 +171,17 @@ end # copied from Carina.jl function _compute_stable_dt(asm, p, CFL) fspace = FEC.function_space(asm.dof) - - # Pre-allocate per-block storage for element char lengths (nq × nelem) - # char_len_storage = Array{Float64, 3}[] - # for (b, ref_fe) in enumerate(fspace.ref_fes) - # for b in 1:FEC.num_blocks(fspace) - # ref_fe = FEC.block_reference_element(fspace, b) - # nquad = num_cell_quadrature_points(ref_fe) - # nelem = num_elements(fspace, b) - # push!(char_len_storage, zeros(Float64, 1, nquad, nelem)) - # end - # char_len_storage = NamedTuple{keys(fspace.ref_fes)}(char_len_storage) - char_len_storage = L2Field(undef, Float64, 1, FEC.block_quadrature_sizes(fspace)) - fill!(char_len_storage, 0.0) + stable_dts = L2Field(undef, Float64, 1, FEC.block_quadrature_sizes(fspace)) + fill!(stable_dts, 0.0) # Assemble per-element char lengths on device U_zeros = zeros(Float64, length(asm.dof.unknown_dofs)) FEC.assemble_quadrature_quantity!( - char_len_storage, nothing, asm.dof, - characteristic_element_length, + stable_dts, nothing, asm.dof, + stable_time_step, U_zeros, p ) - - # Min-reduction over all blocks - stable_dt = Inf - for (b, (block_physics, props)) in enumerate(zip( - values(p.physics), values(p.properties), - )) - block_storage = FEC.block_view(char_len_storage, b) - ρ = props[1] - M = p_wave_modulus(block_physics.constitutive_model, props) - c_p = sqrt(M / ρ) - h_min = minimum(block_storage) # GPU-native reduction if on device - block_dt = CFL * h_min / c_p - stable_dt = min(stable_dt, block_dt) - end + stable_dt = reduce(min, stable_dts) return stable_dt end diff --git a/test/TestExplicitDynamics.jl b/test/TestExplicitDynamics.jl index 5b2560b..ff9b3ba 100644 --- a/test/TestExplicitDynamics.jl +++ b/test/TestExplicitDynamics.jl @@ -9,7 +9,7 @@ function test_explicit_dynamics() # Physics physics = (; sphere = SolidMechanics( - ThreeDimensional(), NeoHookean() + ThreeDimensional(), Hyperelastic(NeoHookean()) ) ) props = (; diff --git a/test/TestHeatConduction.jl b/test/TestHeatConduction.jl index 75283eb..665e6ee 100644 --- a/test/TestHeatConduction.jl +++ b/test/TestHeatConduction.jl @@ -5,10 +5,10 @@ function test_heat_conduction() times = TimeStepper(0.0, 75.0, 100) physics = (; - unnamed_block_1 = Cthonios.HeatConduction() + var"" = Cthonios.HeatConduction() ) props = (; - unnamed_block_1 = Dict{String, Any}( + var"" = Dict{String, Any}( "density" => 1.0, "heat capacity" => 20.0, "thermal conductivity" => 1.0 diff --git a/test/TestImplicitDynamics.jl b/test/TestImplicitDynamics.jl index b015689..65b6e31 100644 --- a/test/TestImplicitDynamics.jl +++ b/test/TestImplicitDynamics.jl @@ -9,12 +9,12 @@ function test_implicit_dynamics() # Physics physics = (; - cube = SolidMechanics( - ThreeDimensional(), NeoHookean() + sphere = SolidMechanics( + ThreeDimensional(), Hyperelastic(NeoHookean()) ) ) props = (; - cube = Dict{String, Any}( + sphere = Dict{String, Any}( "density" => 1.e3, "Young's modulus" => 1.e4, "Poisson's ratio" => 0.33 diff --git a/test/TestQuasiStatic.jl b/test/TestQuasiStatic.jl index d6d76e6..de47889 100644 --- a/test/TestQuasiStatic.jl +++ b/test/TestQuasiStatic.jl @@ -8,7 +8,7 @@ function test_quasistatic() # Physics physics = (; - Block1 = SolidMechanics(PlaneStrain(), NeoHookean()) + Block1 = SolidMechanics(PlaneStrain(), Hyperelastic(NeoHookean())) ) props = (; Block1 = Dict{String, Any}(