From dd47257baaa8b2ca7ccf7402c8a116da30cc4212 Mon Sep 17 00:00:00 2001 From: "Craig M. Hamel" Date: Thu, 6 Aug 2026 23:34:22 -0400 Subject: [PATCH 1/2] Adding PropertyFieldView type so we don't need to resort to StaticArrays. The only things StaticArrays offers us isbits-ness. But this PropertyFieldView type is also isbits. If we only need properties for e.g. elastoplastic models, then StaticArrays is just fine since we'll usually be under ~10 parameters. However, for viscoelastic models, we could easily get into the dozens or even over a hundred individual parameters where StaticArrays become less competetive and can lead to subtsantial compile times. This PR slightly reverts some things but not by much. Tested against the whole Carina.jl test suite. --- src/Enzyme.jl | 4 ++-- src/Fields.jl | 16 +++++++++++----- src/assemblers/Assemblers.jl | 10 +++++++--- src/assemblers/MatrixAction.jl | 2 +- 4 files changed, 21 insertions(+), 11 deletions(-) diff --git a/src/Enzyme.jl b/src/Enzyme.jl index 651fbe3..e6cc8a8 100644 --- a/src/Enzyme.jl +++ b/src/Enzyme.jl @@ -56,7 +56,7 @@ function _assemble_scalar_block_enzyme_safe!( for e in axes(state_old, 3) conn = connectivity(ref_fe, conns, e, coffset) x_el, u_el, u_el_old = element_level_fields(ref_fe, conn, X, U, U_old) - props_el = properties(props, physics, e, b) + props_el = properties(props, e, b) # val_el = _element_scratch(return_type, ref_fe, U) for q in 1:num_cell_quadrature_points(ref_fe) @@ -168,7 +168,7 @@ function _assemble_vector_block_enzyme_safe!( conn = connectivity(ref_fe, conns, e, coffset) x_el, u_el, u_el_old = element_level_fields(ref_fe, conn, X, U, U_old) - props_el = properties(props, physics, e, b) + props_el = properties(props, e, b) # # val_el = _element_scratch(return_type, ref_fe, U) for q in 1:num_cell_quadrature_points(ref_fe) diff --git a/src/Fields.jl b/src/Fields.jl index 284038e..f7a2001 100644 --- a/src/Fields.jl +++ b/src/Fields.jl @@ -552,16 +552,22 @@ function properties(field::PropertyField, e::Int, b::Int) nfields = num_fields(field, b) if field.isblockconstant[b] == PROPS_CONST start = offset - finish = offset + nfields - 1 elseif field.isblockconstant[b] == PROPS_ELEMS start = offset + nfields * (e - 1) - finish = offset + nfields * e - 1 - else - @assert false "Should never happen" end - return view(field.data, start:finish) + return PropertyFieldView(field.data, start, nfields) end +struct PropertyFieldView{D <: AbstractVector} <: AbstractVector{eltype(D)} + data::D + start::Int + len::Int +end + +Base.size(v::PropertyFieldView) = (v.len,) +Base.length(v::PropertyFieldView) = v.len +Base.@propagate_inbounds Base.getindex(v::PropertyFieldView, i::Int) = v.data[v.start + i - 1] + ###################################################################################################### # StateVariableField ###################################################################################################### diff --git a/src/assemblers/Assemblers.jl b/src/assemblers/Assemblers.jl index dfd2754..e6580b1 100644 --- a/src/assemblers/Assemblers.jl +++ b/src/assemblers/Assemblers.jl @@ -199,7 +199,11 @@ $(TYPEDSIGNATURES) offset = field.offsets[b] ec = ifelse(field.isblockconstant[b] == PROPS_CONST, 1, e) base = offset + NP * (ec - 1) - return SVector{NP, eltype(field)}(ntuple(i -> field.data[base + i - 1], NP)) + # return SVector{NP, eltype(field)}(ntuple(i -> field.data[base + i - 1], NP)) + nprops = num_fields(field, b) + # return view(field.data, base:base + NP - 1) + # return PropertyFieldView(field.data, base, NP) + return PropertyFieldView(field.data, base, nprops) end @inline function _element_level_properties(props::AbstractArray, ::Int) @@ -431,7 +435,7 @@ function _assemble_block!( fec_foraxes(state_old, 3) do e conn = connectivity(ref_fe, conns, e, coffset) x_el, u_el, u_el_old = element_level_fields(ref_fe, conn, X, U, U_old) - props_el = properties(props, physics, e, b) + props_el = properties(props, e, b) val_el = _element_scratch(return_type, ref_fe, U) for q in 1:num_cell_quadrature_points(ref_fe) interps = _cell_interpolants(ref_fe, q) @@ -463,7 +467,7 @@ function _assemble_block!( fec_foraxes(state_old, 3) do e conn = connectivity(ref_fe, conns, e, coffset) x_el, u_el, u_el_old = element_level_fields(ref_fe, conn, X, U, U_old) - props_el = properties(props, physics, e, b) + props_el = properties(props, e, b) for q in 1:num_cell_quadrature_points(ref_fe) interps = _cell_interpolants(ref_fe, q) state_old_q = _quadrature_level_state(state_old, q, e) diff --git a/src/assemblers/MatrixAction.jl b/src/assemblers/MatrixAction.jl index 5dc1ea7..e45d3cd 100644 --- a/src/assemblers/MatrixAction.jl +++ b/src/assemblers/MatrixAction.jl @@ -66,7 +66,7 @@ function _assemble_block_matrix_free_action!( fec_foraxes(state_old, 3) do e conn = connectivity(ref_fe, conns, e, coffset) x_el, u_el, u_el_old, v_el = element_level_fields(ref_fe, conn, X, U, U_old, V) - props_el = properties(props, physics, e, b) + props_el = properties(props, e, b) Kv_el = _element_scratch(AssembledVector(), ref_fe, U) for q in 1:num_cell_quadrature_points(ref_fe) interps = _cell_interpolants(ref_fe, q) From c0e4cfea66835a8341f45ab4caa570e066b48403 Mon Sep 17 00:00:00 2001 From: Alejandro Mota Date: Fri, 7 Aug 2026 11:43:58 -0700 Subject: [PATCH 2/2] Give PropertyFieldView a concrete eltype and a bounds check `PropertyFieldView{D} <: AbstractVector{eltype(D)}` does not do what it reads like: `eltype` of an unbound TypeVar is `Any`, so the supertype was `AbstractVector{Any}` and `eltype(properties(...))` came back `Any`. That is not cosmetic downstream. ConstitutiveModels' `module_props` builds `SVector{NP, eltype(props)}`, and every `Hyperelastic` and `LinearElastic` evaluation routes through it, so an `Any` eltype turned the whole constitutive call into a boxed, dynamically dispatched one. Measured through Carina's assembly loop on torsion-qs (160k elements, 527,877 DOF, 12 threads), against this branch's parent: stiffness_action 25.19 ms -> 551.38 ms (21.9x slower) residual 23.77 ms -> 414.50 ms (17.4x slower) GPU was unaffected -- 9.83 ms vs 9.70 ms, checksums bit-identical -- because GPUCompiler inlines aggressively enough that SROA recovers the concrete type before anything can box. Carrying the element type as a parameter restores CPU to 24.84 ms and 21.63 ms, with the same checksums. Separately, `getindex` carried `@propagate_inbounds` but had no `@boundscheck`, so no bounds check existed at any optimization level. Since every block's properties share one flat vector, an out-of-range read silently returned the *next* block's properties rather than failing: on a three-property block, `p[4]` handed back block 2's first property. Blocks may now have different property counts, so this is reachable from an ordinary off-by-one in a material model. Tests cover both, and `IndexStyle` is declared linear so the generic AbstractArray fallbacks stop routing through CartesianIndices. Signed-off-by: Alejandro Mota --- src/Fields.jl | 8 ++++++-- test/TestFields.jl | 35 +++++++++++++++++++++++++++++++++++ 2 files changed, 41 insertions(+), 2 deletions(-) diff --git a/src/Fields.jl b/src/Fields.jl index f7a2001..3d2841f 100644 --- a/src/Fields.jl +++ b/src/Fields.jl @@ -558,7 +558,7 @@ function properties(field::PropertyField, e::Int, b::Int) return PropertyFieldView(field.data, start, nfields) end -struct PropertyFieldView{D <: AbstractVector} <: AbstractVector{eltype(D)} +struct PropertyFieldView{T, D <: AbstractVector{T}} <: AbstractVector{T} data::D start::Int len::Int @@ -566,7 +566,11 @@ end Base.size(v::PropertyFieldView) = (v.len,) Base.length(v::PropertyFieldView) = v.len -Base.@propagate_inbounds Base.getindex(v::PropertyFieldView, i::Int) = v.data[v.start + i - 1] +Base.IndexStyle(::Type{<:PropertyFieldView}) = IndexLinear() +Base.@propagate_inbounds function Base.getindex(v::PropertyFieldView, i::Int) + @boundscheck checkbounds(v, i) + return @inbounds v.data[v.start + i - 1] +end ###################################################################################################### # StateVariableField diff --git a/test/TestFields.jl b/test/TestFields.jl index 96aa5a0..847d453 100644 --- a/test/TestFields.jl +++ b/test/TestFields.jl @@ -281,6 +281,41 @@ end @test all(FiniteElementContainers.properties(props, 1, 1) .≈ original) end +@testitem "Fields - test_property_field_view_eltype_is_concrete" begin + using StaticArrays + # `PropertyFieldView` must carry the element type as a parameter. Declaring + # it as `AbstractVector{eltype(D)}` over a bare `D` silently yields + # `AbstractVector{Any}`, because `eltype` of an unbound TypeVar is `Any`. + # Downstream that is not cosmetic: ConstitutiveModels' `module_props` builds + # `SVector{NP, eltype(props)}`, so an `Any` eltype turns every constitutive + # evaluation into a boxed, dynamically dispatched call. + props = FiniteElementContainers.PropertyField([[1.0, 2.0, 3.0]]) + view = FiniteElementContainers.properties(props, 1, 1) + @test eltype(view) === Float64 + @test view isa AbstractVector{Float64} + @test eltype(collect(view)) === Float64 + @test isconcretetype(eltype(SVector{2, eltype(view)}(view[1], view[2]))) + @test Base.IndexStyle(typeof(view)) === IndexLinear() +end + +@testitem "Fields - test_property_field_view_is_bounds_checked" begin + # Every block's properties live in one flat vector, so an unchecked + # out-of-range read returns the *next* block's properties instead of + # failing. Two blocks with different property counts make that concrete. + props = FiniteElementContainers.PropertyField([[1.0, 2.0, 3.0], [10.0, 20.0, 30.0, 40.0]]) + block1 = FiniteElementContainers.properties(props, 1, 1) + @test length(block1) == 3 + @test block1[3] ≈ 3.0 + # Without a bounds check this returns 10.0 -- block 2's first property. + @test_throws BoundsError block1[4] + @test_throws BoundsError block1[0] + + block2 = FiniteElementContainers.properties(props, 1, 2) + @test length(block2) == 4 + @test block2[4] ≈ 40.0 + @test_throws BoundsError block2[5] +end + @testitem "Fields - test_state_variable_field" begin a1 = rand(2, 3, 40) a2 = rand(3, 4, 10)