From 0d3ab53bdb93d5a9e953301cadccebc39bf39c67 Mon Sep 17 00:00:00 2001 From: "Craig M. Hamel" Date: Thu, 6 Aug 2026 00:05:46 -0400 Subject: [PATCH 1/2] Hooking up PropertyField to be the required interface moving forward. It stores all properties as a flat vector and can allow for some or all of the blocks in the mesh to have element specific properties rather than the properties being constant across the whole block. The indexing behavior that supports this also lives on the GPU as Cu/RocM/Vectors of Ints --- Project.toml | 5 +- src/Enzyme.jl | 24 +++++---- src/Fields.jl | 73 ++++++++++++++++++++------- src/Parameters.jl | 72 ++++++++++++++------------ src/Physics.jl | 4 +- src/Utils.jl | 4 +- src/assemblers/Assemblers.jl | 6 ++- src/assemblers/Diagonal.jl | 9 ++-- src/assemblers/LumpedMass.jl | 6 ++- src/assemblers/Matrix.jl | 9 ++-- src/assemblers/MatrixAction.jl | 31 ++++++++---- src/assemblers/QuadratureQuantity.jl | 15 +++--- src/assemblers/Vector.jl | 11 ++-- test/TestAssemblers.jl | 2 +- test/TestBlockOrdering.jl | 2 - test/TestFields.jl | 24 +++++++++ test/mechanics/TestMechanicsCommon.jl | 2 +- 17 files changed, 200 insertions(+), 99 deletions(-) diff --git a/Project.toml b/Project.toml index b7143cd9..dc433bc2 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "FiniteElementContainers" uuid = "d08262e4-672f-4e7f-a976-f2cea5767631" -version = "0.15.0" +version = "0.15.1" authors = ["Craig M. Hamel and contributors"] [deps] @@ -73,6 +73,7 @@ TimerOutputs = "0.5" julia = "1" [extras] +AMDGPU = "21141c5a-9bdb-4563-92ae-f87d6854732e" Gmsh = "705231aa-382f-11e9-3f0c-b7cb4346fdeb" JuliaC = "acedd4c2-ced6-4a15-accc-2607eb759ba2" PartitionedArrays = "5a9dfac6-5c52-46f7-8278-5e2210713be9" @@ -81,4 +82,4 @@ TestItemRunner = "f8b46487-2199-4994-9208-9a1283c18c0a" TestItems = "1c621080-faea-4a02-84b6-bbd5e436b8fe" [targets] -test = ["Gmsh", "JuliaC", "PartitionedArrays", "Test", "TestItemRunner", "TestItems"] +test = ["AMDGPU", "Gmsh", "JuliaC", "PartitionedArrays", "Test", "TestItemRunner", "TestItems"] diff --git a/src/Enzyme.jl b/src/Enzyme.jl index e825895a..e6cc8a8e 100644 --- a/src/Enzyme.jl +++ b/src/Enzyme.jl @@ -13,20 +13,21 @@ function assemble_scalar_enzyme_safe!( conns = fspace.elem_conns # foreach_block(fspace, p) do physics, props, ref_fe, b for (b, ( - block_physics, ref_fe, props + block_physics, ref_fe )) in enumerate(zip( - values(p.physics), values(fspace.ref_fes), - values(p.properties) + values(p.physics), values(fspace.ref_fes) )) _assemble_scalar_block_enzyme_safe!( KA.CPU(), block_view(storage, b), conns.data, conns.offsets[b], func, + b, block_physics, ref_fe, X, t, Δt, U, U_old, - block_view(p.state_old, b), block_view(p.state_new, b), props, + block_view(p.state_old, b), block_view(p.state_new, b), + p.properties, return_type ) end @@ -38,6 +39,7 @@ function _assemble_scalar_block_enzyme_safe!( field, conns::Conn, coffset::Int, func::Function, + b::Int, physics::AbstractPhysics, ref_fe::ReferenceFE, X::AbstractField, t::T, dt::T, U::Solution, U_old::Solution, @@ -54,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 = _element_level_properties(props, e) + 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) @@ -111,10 +113,9 @@ function assemble_vector_enzyme_safe!( # return_type = AssembledVector() conns = fspace.elem_conns for (b, ( - block_physics, ref_fe, props + block_physics, ref_fe )) in enumerate(zip( - values(p.physics), values(fspace.ref_fes), - values(p.properties) + values(p.physics), values(fspace.ref_fes) )) _assemble_vector_block_enzyme_safe!( # KA.get_backend(storage), @@ -122,10 +123,12 @@ function assemble_vector_enzyme_safe!( storage, conns.data, conns.offsets[b], func, + b, block_physics, ref_fe, X, t, Δt, U, U_old, - block_view(p.state_old, b), block_view(p.state_new, b), props, + block_view(p.state_old, b), block_view(p.state_new, b), + p.properties # return_type ) end @@ -147,6 +150,7 @@ function _assemble_vector_block_enzyme_safe!( field, conns::Conn, coffset::Int, func::Function, + b::Int, physics::AbstractPhysics, ref_fe::ReferenceFE, X::AbstractField, t::T, dt::T, U::Solution, U_old::Solution, @@ -164,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 = _element_level_properties(props, e) + 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 3e02c5e9..a5a0eaf2 100644 --- a/src/Fields.jl +++ b/src/Fields.jl @@ -463,42 +463,78 @@ const PROPS_ELEMS = -2 struct PropertyField{ T <: Number, - D <: AbstractVector{T} + D <: AbstractVector{T}, + I <: AbstractVector{Int} } <: AbstractDiscontinuousField{T, D} data::D - isblockconstant::Vector{Int} + isblockconstant::I nblocks::Int - nfields::Vector{Int} - nepes::Vector{Int} - nelems::Vector{Int} - offsets::Vector{Int} + nepes::I + nelems::I + offsets::I - function PropertyField{T, D}(data, isblockconstant, nblocks, nfields, nepes, nelems, offsets) where {T, D} - new{T, D}(data, isblockconstant, nblocks, nfields, nepes, nelems, offsets) + function PropertyField{T, D, I}(data, isblockconstant, nblocks, nepes, nelems, offsets) where {T, D, I} + new{T, D, I}(data, isblockconstant, nblocks, nepes, nelems, offsets) end + # case where all blocks have properties the same throughout function PropertyField(arrs::Vector{<:Vector{T}}) where T <: Number data = reduce(vcat, arrs) isblockconstant = PROPS_CONST * ones(Int, length(arrs)) nblocks = length(arrs) - nfields = map(length, arrs) - nepes = -1 * ones(Int, nblocks) + nepes = map(length, arrs) nelems = -1 * ones(Int, nblocks) offsets = Vector{Int}(undef, 0) offset = 1 for n in axes(arrs, 1) push!(offsets, offset) - offset += nfields[n] + offset += nepes[n] end - return PropertyField{T, typeof(data)}(data, isblockconstant, nblocks, nfields, nepes, nelems, offsets) + return PropertyField{T, typeof(data), typeof(isblockconstant)}(data, isblockconstant, nblocks, nepes, nelems, offsets) + end + + # case where we have mixed constant or none constant + function PropertyField(arrs::Vector{<:Array{T}}) where T <: Number + data = mapreduce(vec, vcat, arrs) + isblockconstant = map(x -> begin + if isa(x, Matrix) + return PROPS_ELEMS + elseif isa(x, Vector) + return PROPS_CONST + else + @assert false "Property arrays should be Vector or Matrix" + end + end, arrs) + nblocks = length(arrs) + out = map(x -> begin + if isa(x, Matrix) + return size(x) + elseif isa(x, Vector) + return (length(x), -1) + end + end, arrs) + nepes = map(x -> x[1], out) + nelems = map(x -> x[2], out) + offsets = Vector{Int}(undef, 0) + offset = 1 + for n in axes(arrs, 1) + push!(offsets, offset) + if isa(arrs[n], Matrix) + offset += nepes[n] * nelems[n] + elseif isa(arrs[n], Vector) + offset += nepes[n] + end + end + return PropertyField{T, typeof(data), typeof(isblockconstant)}(data, isblockconstant, nblocks, nepes, nelems, offsets) end end function Adapt.adapt_structure(to, field::PropertyField) data = adapt(to, field.data) - return PropertyField{eltype(field), typeof(data)}( + isblockconstant = adapt(to, field.isblockconstant) + return PropertyField{eltype(field), typeof(data), typeof(isblockconstant)}( data, - adapt(to, field.isblockconstant), + isblockconstant, field.nblocks, adapt(to, field.nepes), adapt(to, field.nelems), @@ -507,19 +543,22 @@ function Adapt.adapt_structure(to, field::PropertyField) end function num_fields(field::PropertyField, b::Int) - return field.nfields[b] + return field.nepes[b] end function properties(field::PropertyField, e::Int, b::Int) offset = field.offsets[b] nfields = num_fields(field, b) if field.isblockconstant[b] == PROPS_CONST - return view(field.data, offset:offset + nfields - 1) + start = offset + finish = offset + nfields - 1 elseif field.isblockconstant[b] == PROPS_ELEMS - @assert false finish me + start = offset + nfields * (e - 1) + finish = offset + nfields * e - 1 else @assert false "Should never happen" end + return view(field.data, start:finish) end ###################################################################################################### diff --git a/src/Parameters.jl b/src/Parameters.jl index 5e79a347..479716aa 100644 --- a/src/Parameters.jl +++ b/src/Parameters.jl @@ -51,6 +51,26 @@ function _align_blocks(fspace, x, what) return NamedTuple{names}(ntuple(_ -> x, length(names))) end +# a single properties object shared by every block +# needs to be constant props, can't be element level +# unless we have one block, but let's not specialize that muc +function _setup_properties(fspace, props::Vector) + return PropertyField(map(_ -> props, block_names(fspace))) +end + +# namedtuple case that should become deprecated soon +function _setup_properties(fspace, props::NamedTuple) + names = tuple(Symbol.(block_names(fspace))...) + _check_block_keys(keys(props), names, "properties") + return PropertyField([map(x -> getfield(props, x), names)...]) +end + +function _setup_properties(fspace, props::Dict{String}) + names = block_names(fspace) + _check_block_keys(keys(props), names, "properties") + return PropertyField([map(x -> props[x], names)...]) +end + function _setup_state_variables(fspace, physics) state_old = Array{Float64, 3}[] state_new = Array{Float64, 3}[] @@ -88,6 +108,7 @@ $(TYPEDSIGNATURES) $(TYPEDFIELDS) """ struct Parameters{ + D, # dimension IT <: Integer, RT <: Number, IV <: AbstractVector{IT}, @@ -103,8 +124,6 @@ struct Parameters{ PBCFuncs <: AbstractVector, RBCFuncs <: AbstractVector, Phys, - Props, - Coords <: AbstractField, Field <: AbstractField } <: AbstractParameters ics::InitialConditions{ICFuncs, IV, RV} @@ -115,10 +134,10 @@ struct Parameters{ sources::Sources{SRCFuncs, RM4} times::TimeStepper{RT} physics::Phys - properties::Props + properties::PropertyField{RT, RV, IV} state_old::StateVariableField{RT, RV} state_new::StateVariableField{RT, RV} - coords::Coords + coords::H1Field{RT, RV, D} field::Field field_old::Field # scratch fields @@ -150,7 +169,7 @@ function Parameters( # for mixed spaces we'll need to do this more carefully physics = _align_blocks(fspace, physics, "physics") - properties = _align_blocks(fspace, properties, "properties") + properties = _setup_properties(fspace, properties) # setup state variables state_old, state_new = _setup_state_variables(fspace, physics) @@ -172,19 +191,6 @@ function Parameters( end function Adapt.adapt_structure(to, p::Parameters) - - # need to handle props specially - props = [] - for p in values(p.properties) - if isa(p, SArray) - push!(props, p) - else - push!(props, adapt(to, p)) - end - end - - props = NamedTuple{keys(p.properties)}(props) - return Parameters( adapt(to, p.ics), adapt(to, p.dirichlet_bcs), @@ -194,7 +200,7 @@ function Adapt.adapt_structure(to, p::Parameters) adapt(to, p.sources), adapt(to, p.times), adapt(to, p.physics), - props, + adapt(to, p.properties), adapt(to, p.state_old), adapt(to, p.state_new), adapt(to, p.coords), @@ -233,7 +239,7 @@ function KA.get_backend(p::Parameters) end struct TypeStableParameters{ - # Funcs <: AbstractVector, + D, # dimension SFuncT, VFuncT, IT <: Integer, @@ -242,8 +248,6 @@ struct TypeStableParameters{ RV <: AbstractVector{RT}, RM <: AbstractMatrix{<:SVector}, Phys, - Props, - Coords <: AbstractField, Field <: AbstractField } <: AbstractParameters ics::InitialConditions{Vector{InitialConditionFunction{SFuncT}}, IV, RV} @@ -254,16 +258,16 @@ struct TypeStableParameters{ sources::Sources{Vector{SourceFunction{VFuncT}}, RM} times::TimeStepper{RT} physics::Phys - properties::Props + properties::PropertyField{RT, RV, IV} state_old::StateVariableField{RT, RV} state_new::StateVariableField{RT, RV} - coords::Coords + coords::H1Field{RT, RV, D} field::Field field_old::Field # scratch fields hvp_scratch_field::Field - function TypeStableParameters{SF, VF}(mesh, assembler, physics, props, ics, dbcs, nbcs, pbcs, srcs, times) where {SF, VF} + function TypeStableParameters{D, SF, VF}(mesh, assembler, physics, props, ics, dbcs, nbcs, pbcs, srcs, times) where {D, SF, VF} dof = assembler.dof ND = size(dof, 1) fspace = function_space(dof) @@ -274,7 +278,8 @@ struct TypeStableParameters{ srcs = Sources{VF}(mesh, dof, srcs) physics = _align_blocks(fspace, physics, "physics") - props = _align_blocks(fspace, props, "properties") + # props = _align_blocks(fspace, props, "properties") + props = _setup_properties(fspace, props) state_old, state_new = _setup_state_variables(fspace, physics) @@ -287,8 +292,8 @@ struct TypeStableParameters{ update_dofs!(assembler, dbcs, pbcs) new{ - SF, VF, Int, Float64, Vector{Int}, Vector{Float64}, Matrix{SVector{ND, Float64}}, - typeof(physics), typeof(props), typeof(mesh.nodal_coords), typeof(field) + D, SF, VF, Int, Float64, Vector{Int}, Vector{Float64}, Matrix{SVector{ND, Float64}}, + typeof(physics), typeof(field) }( ics, dbcs, nbcs, pbcs, srcs, times, @@ -296,10 +301,10 @@ struct TypeStableParameters{ ) end - function TypeStableParameters{SF, VF}( + function TypeStableParameters{D, SF, VF}( mesh, assembler, physics, props, state_old, state_new, ics, dbcs, nbcs, pbcs, srcs, times - ) where {SF, VF} + ) where {D, SF, VF} dof = assembler.dof ND = size(dof, 1) fspace = function_space(dof) @@ -310,7 +315,8 @@ struct TypeStableParameters{ srcs = Sources{VF}(mesh, dof, srcs) physics = _align_blocks(fspace, physics, "physics") - props = _align_blocks(fspace, props, "properties") + # props = _align_blocks(fspace, props, "properties") + props = _setup_properties(fspace, props) coords = mesh.nodal_coords field = create_field(assembler) @@ -321,8 +327,8 @@ struct TypeStableParameters{ update_dofs!(assembler, dbcs, pbcs) new{ - SF, VF, Int, Float64, Vector{Int}, Vector{Float64}, Matrix{SVector{ND, Float64}}, - typeof(physics), typeof(props), typeof(mesh.nodal_coords), typeof(field) + D, SF, VF, Int, Float64, Vector{Int}, Vector{Float64}, Matrix{SVector{ND, Float64}}, + typeof(physics), typeof(field) }( ics, dbcs, nbcs, pbcs, srcs, times, diff --git a/src/Physics.jl b/src/Physics.jl index e3ebcb4e..b8f5d6c9 100644 --- a/src/Physics.jl +++ b/src/Physics.jl @@ -30,7 +30,7 @@ end # default function create_initial_state(::AbstractPhysics{NF, NP, 0}) where {NF, NP} - return SVector{0, Float64}() + return zeros(0) end """ @@ -43,7 +43,7 @@ end # default function create_properties(::AbstractPhysics{NF, 0, NS}) where {NF, NS} - return SVector{0, Float64}() + return zeros(0) end # default diff --git a/src/Utils.jl b/src/Utils.jl index 6ccf4d0b..257e90d2 100644 --- a/src/Utils.jl +++ b/src/Utils.jl @@ -180,7 +180,7 @@ end for k in 1:N push!(stmts.args, quote f( - values(p.physics)[$k], values(p.properties)[$k], + values(p.physics)[$k], block_reference_element(fspace, $k), $k ) end) @@ -201,7 +201,7 @@ end ref_dispatches = map(1:n_refs) do j quote if fspace.block_to_ref_fe_id[$i] == $j - f(p.physics[$i], p.properties[$i], fspace.ref_fes[$j], $i) + f(p.physics[$i], fspace.ref_fes[$j], $i) end end end diff --git a/src/assemblers/Assemblers.jl b/src/assemblers/Assemblers.jl index 92e0bbfc..9a65faf8 100644 --- a/src/assemblers/Assemblers.jl +++ b/src/assemblers/Assemblers.jl @@ -403,6 +403,7 @@ function _assemble_block!( field, conns::Conn, coffset::Int, func::Function, + b::Int, # block index physics::AbstractPhysics, ref_fe::ReferenceFE, X::AbstractField, t::T, dt::T, U::Solution, U_old::Solution, @@ -418,7 +419,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 = _element_level_properties(props, e) + 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) @@ -434,6 +435,7 @@ end function _assemble_block!( field, func!::Function, + b::Int, # block index physics::AbstractPhysics, t::T, Δt::T, props::P, state_old::S, state_new::S, @@ -449,7 +451,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 = _element_level_properties(props, e) + 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/Diagonal.jl b/src/assemblers/Diagonal.jl index 91649d2a..d3ecbb2c 100644 --- a/src/assemblers/Diagonal.jl +++ b/src/assemblers/Diagonal.jl @@ -43,14 +43,15 @@ function assemble_diagonal!( _update_for_assembly!(p, dof, Uu) return_type = AssembledDiagonal() conns = fspace.elem_conns - foreach_block(fspace, p) do physics, props, ref_fe, b + foreach_block(fspace, p) do physics, ref_fe, b if use_inplace_methods _assemble_block!( storage, func, + b, physics, t, Δt, - props, + p.properties, block_view(p.state_old, b), block_view(p.state_new, b), conns.data, conns.offsets[b], ref_fe, X, U, U_old ) @@ -59,10 +60,12 @@ function assemble_diagonal!( storage, conns.data, conns.offsets[b], func, + b, physics, ref_fe, X, t, Δt, U, U_old, - block_view(p.state_old, b), block_view(p.state_new, b), props, + block_view(p.state_old, b), block_view(p.state_new, b), + p.properties, return_type ) end diff --git a/src/assemblers/LumpedMass.jl b/src/assemblers/LumpedMass.jl index b37834a5..49d38a4b 100644 --- a/src/assemblers/LumpedMass.jl +++ b/src/assemblers/LumpedMass.jl @@ -41,15 +41,17 @@ function assemble_lumped_mass!( _update_for_assembly!(p, dof, Uu) return_type = AssembledVector() conns = fspace.elem_conns - foreach_block(fspace, p) do physics, props, ref_fe, b + foreach_block(fspace, p) do physics, ref_fe, b _assemble_block!( storage, conns.data, conns.offsets[b], func, + b, physics, ref_fe, X, t, Δt, U, U_old, - block_view(p.state_old, b), block_view(p.state_new, b), props, + block_view(p.state_old, b), block_view(p.state_new, b), + p.properties, return_type ) end diff --git a/src/assemblers/Matrix.jl b/src/assemblers/Matrix.jl index dc670344..b86e86c9 100644 --- a/src/assemblers/Matrix.jl +++ b/src/assemblers/Matrix.jl @@ -46,14 +46,15 @@ function assemble_matrix!( _update_for_assembly!(p, dof, Uu) return_type = AssembledMatrix() conns = fspace.elem_conns - foreach_block(fspace, p) do physics, props, ref_fe, b + foreach_block(fspace, p) do physics, ref_fe, b if use_inplace_methods _assemble_block!( block_view(storage, pattern, b), func, + b, physics, t, dt, - props, + p.properties, block_view(p.state_old, b), block_view(p.state_new, b), conns.data, conns.offsets[b], ref_fe, X, U, U_old ) @@ -62,10 +63,12 @@ function assemble_matrix!( block_view(storage, pattern, b), conns.data, conns.offsets[b], func, + b, physics, ref_fe, X, t, dt, U, U_old, - block_view(p.state_old, b), block_view(p.state_new, b), props, + block_view(p.state_old, b), block_view(p.state_new, b), + p.properties, return_type ) end diff --git a/src/assemblers/MatrixAction.jl b/src/assemblers/MatrixAction.jl index 0ac446a5..b7109a22 100644 --- a/src/assemblers/MatrixAction.jl +++ b/src/assemblers/MatrixAction.jl @@ -33,15 +33,17 @@ function assemble_matrix_free_action!( V = p.hvp_scratch_field _update_for_assembly!(p, dof, Uu, Vu) conns = fspace.elem_conns - foreach_block(fspace, p) do physics, props, ref_fe, b + foreach_block(fspace, p) do physics, ref_fe, b _assemble_block_matrix_free_action!( storage, conns.data, conns.offsets[b], func_action, + b, physics, ref_fe, X, t, Δt, U, U_old, V, - block_view(p.state_old, b), block_view(p.state_new, b), props + block_view(p.state_old, b), block_view(p.state_new, b), + p.properties ) end end @@ -50,6 +52,7 @@ function _assemble_block_matrix_free_action!( field::AbstractField, conns::Conn, coffset, func_action::Function, + b::Int, physics::AbstractPhysics, ref_fe::ReferenceFE, X::AbstractField, t::T, Δt::T, U::Solution, U_old::Solution, V::Solution, @@ -63,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 = _element_level_properties(props, e) + 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) @@ -123,15 +126,17 @@ function assemble_matrix_free_action_full!( V = p.hvp_scratch_field _update_for_assembly_full!(p, U_full, v_full) conns = fspace.elem_conns - foreach_block(fspace, p) do physics, props, ref_fe, b + foreach_block(fspace, p) do physics, ref_fe, b _assemble_block_matrix_free_action!( storage, conns.data, conns.offsets[b], func_action, + b, physics, ref_fe, X, t, Δt, U, U_old, V, - block_view(p.state_old, b), block_view(p.state_new, b), props + block_view(p.state_old, b), block_view(p.state_new, b), + p.properties ) end # The free-DOF entry point above relies on `p.hvp_scratch_field`'s BC @@ -180,14 +185,15 @@ function assemble_matrix_action!( V = p.hvp_scratch_field _update_for_assembly!(p, dof, Uu, Vu) conns = fspace.elem_conns - foreach_block(fspace, p) do physics, props, ref_fe, b + foreach_block(fspace, p) do physics, ref_fe, b if use_inplace_methods _assemble_block_matrix_action!( storage, func, + b, physics, t, Δt, - props, + p.properties, block_view(p.state_old, b), block_view(p.state_new, b), conns.data, conns.offsets[b], ref_fe, X, U, U_old, V ) @@ -196,10 +202,12 @@ function assemble_matrix_action!( storage, conns.data, conns.offsets[b], func, + b, physics, ref_fe, X, t, Δt, U, U_old, V, - block_view(p.state_old, b), block_view(p.state_new, b), props + block_view(p.state_old, b), block_view(p.state_new, b), + p.properties ) end end @@ -209,6 +217,7 @@ function _assemble_block_matrix_action!( field::AbstractField, conns::Conn, coffset, func::Function, + b::Int, physics::AbstractPhysics, ref_fe::ReferenceFE, X::AbstractField, t::T, Δt::T, U::Solution, U_old::Solution, V::Solution, @@ -222,7 +231,8 @@ function _assemble_block_matrix_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 = _element_level_properties(props, e) + # props_el = _element_level_properties(props, e) + props_el = properties(props, e, b) K_el = _element_scratch(AssembledMatrix(), ref_fe, U) for q in 1:num_cell_quadrature_points(ref_fe) interps = _cell_interpolants(ref_fe, q) @@ -240,6 +250,7 @@ end function _assemble_block_matrix_action!( field, func!::Function, + b::Int, physics::AbstractPhysics, t::T, Δt::T, props::P, state_old::S, state_new::S, @@ -255,7 +266,7 @@ function _assemble_block_matrix_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 = _element_level_properties(props, e) + 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/QuadratureQuantity.jl b/src/assemblers/QuadratureQuantity.jl index c865512c..d8b04209 100644 --- a/src/assemblers/QuadratureQuantity.jl +++ b/src/assemblers/QuadratureQuantity.jl @@ -29,15 +29,17 @@ function assemble_quadrature_quantity!( U_old = p.field_old _update_for_assembly!(p, dof, Uu) conns = fspace.elem_conns - foreach_block(fspace, p) do physics, props, ref_fe, b + foreach_block(fspace, p) do physics, ref_fe, b _assemble_block!( block_view(storage, b), conns.data, conns.offsets[b], func, + b, physics, ref_fe, X, t, Δt, U, U_old, - block_view(p.state_old, b), block_view(p.state_new, b), props, + block_view(p.state_old, b), block_view(p.state_new, b), + p.properties, return_type ) end @@ -61,21 +63,22 @@ function assemble_quadrature_quantity!( conns = fspace.elem_conns for (b, ( block_storage, - block_physics, ref_fe, props + block_physics, ref_fe )) in enumerate(zip( values(storage), - values(p.physics), values(fspace.ref_fes), - values(p.properties) + values(p.physics), values(fspace.ref_fes) )) _assemble_block!( # backend, block_storage, conns.data, conns.offsets[b], func, + b, block_physics, ref_fe, X, t, Δt, U, U_old, - block_view(p.state_old, b), block_view(p.state_new, b), props, + block_view(p.state_old, b), block_view(p.state_new, b), + p.properties, return_type ) end diff --git a/src/assemblers/Vector.jl b/src/assemblers/Vector.jl index 1913dee4..443ec306 100644 --- a/src/assemblers/Vector.jl +++ b/src/assemblers/Vector.jl @@ -38,7 +38,8 @@ function assemble_vector!( return_type = AssembledVector() conns = fspace.elem_conns # foreach_block(conns, p.physics, p.properties, fspace.ref_fes) do physics, props, ref_fe, b - foreach_block(fspace, p) do physics, props, ref_fe, b + # foreach_block(fspace, p) do physics, props, ref_fe, b + foreach_block(fspace, p) do physics, ref_fe, b # if use_sparse_vector # field = block_view(storage, pattern, b) # else @@ -50,9 +51,11 @@ function assemble_vector!( _assemble_block!( field, func, + b, physics, t, Δt, - props, + # props, + p.properties, block_view(p.state_old, b), block_view(p.state_new, b), conns.data, conns.offsets[b], ref_fe, X, U, U_old ) @@ -61,10 +64,12 @@ function assemble_vector!( field, conns.data, conns.offsets[b], func, + b, physics, ref_fe, X, t, Δt, U, U_old, - block_view(p.state_old, b), block_view(p.state_new, b), props, + block_view(p.state_old, b), block_view(p.state_new, b), + p.properties, return_type ) end diff --git a/test/TestAssemblers.jl b/test/TestAssemblers.jl index 5eb526a3..50fadb73 100644 --- a/test/TestAssemblers.jl +++ b/test/TestAssemblers.jl @@ -50,7 +50,7 @@ end f(X, _) = 2. * π^2 * sin(π * X[1]) * sin(π * X[2]) bc_func(_, _) = 0. physics = Poisson(f) - props = SVector{0, Float64}() + props = zeros(0) u = ScalarFunction(V, "u") dbcs = DirichletBC[ DirichletBC("u", bc_func; sideset_name = "boundary") diff --git a/test/TestBlockOrdering.jl b/test/TestBlockOrdering.jl index e1be4477..4fe8a9b6 100644 --- a/test/TestBlockOrdering.jl +++ b/test/TestBlockOrdering.jl @@ -157,7 +157,6 @@ end # Single material for the whole mesh: replicated, and keyed by block name. p = create_parameters(mesh, asm, one_physics, one_props) @test collect(keys(p.physics)) == Symbol.(BLOCK_NAMES) - @test collect(keys(p.properties)) == Symbol.(BLOCK_NAMES) # State variables are allocated per block by walking `values(physics)` # against `block_quadrature_size(fspace, b)`, so their element counts are a @@ -169,7 +168,6 @@ end scrambled_props = (b3 = one_props, b1 = one_props, b2 = one_props) p = create_parameters(mesh, asm, scrambled_physics, scrambled_props) @test collect(keys(p.physics)) == Symbol.(BLOCK_NAMES) - @test collect(keys(p.properties)) == Symbol.(BLOCK_NAMES) @test [block_size(p.state_old, b)[3] for b in 1:length(BLOCK_NAMES)] == BLOCK_SIZES # And a mismatch stops the run instead of producing a plausible wrong answer. diff --git a/test/TestFields.jl b/test/TestFields.jl index 23b9a32a..6b10fee5 100644 --- a/test/TestFields.jl +++ b/test/TestFields.jl @@ -213,6 +213,30 @@ end @test all(FiniteElementContainers.properties(props, 100, 2) .≈ props_2) end +@testitem "Fields - test_property_field_,mixed_constant_and_element_level" begin + props_1 = rand(3) + props_2 = rand(4, 20) + props = FiniteElementContainers.PropertyField([props_1, props_2]) + @test all(FiniteElementContainers.properties(props, 1, 1) .≈ props_1) + @test all(FiniteElementContainers.properties(props, 100, 1) .≈ props_1) + + for e in axes(props_2, 2) + @test all(FiniteElementContainers.properties(props, e, 2) .≈ props_2[:, e]) + end +end + +@testitem "Fields - test_property_field_all_element_level" begin + props_1 = rand(3, 10) + props_2 = rand(4, 20) + props = FiniteElementContainers.PropertyField([props_1, props_2]) + for e in axes(props_1, 2) + @test all(FiniteElementContainers.properties(props, e, 1) .≈ props_1[:, e]) + end + for e in axes(props_2, 2) + @test all(FiniteElementContainers.properties(props, e, 2) .≈ props_2[:, e]) + end +end + @testitem "Fields - test_state_variable_field" begin a1 = rand(2, 3, 40) a2 = rand(3, 4, 10) diff --git a/test/mechanics/TestMechanicsCommon.jl b/test/mechanics/TestMechanicsCommon.jl index e6868843..8f64f4d7 100644 --- a/test/mechanics/TestMechanicsCommon.jl +++ b/test/mechanics/TestMechanicsCommon.jl @@ -8,7 +8,7 @@ function FiniteElementContainers.create_properties(::Mechanics) ρ = 1e3 K = 10.e9 G = 1.e9 - return SVector{3, Float64}(ρ, K, G) + return [ρ, K, G] end @inline function strain_energy( From 42a21a49f2890e63be9de04a0606ca4270edd1a7 Mon Sep 17 00:00:00 2001 From: Alejandro Mota Date: Thu, 6 Aug 2026 11:39:17 -0700 Subject: [PATCH 2/2] Make the property accessor GPU safe `properties(field, e, b)` returns `view(field.data, start:finish)`. Building that SubArray inside a device kernel lowers to a dynamic call, so every assembler that reads properties fails to compile for the GPU. GPUCompiler does not surface it as an InvalidIRError either -- it segfaults in check_ir! while decoding the jl_invoke operand -- so the failure reads as a crash rather than a diagnostic. Take the property count from the physics type and return a statically sized SVector, mirroring the GPU-safe `connectivity(ref_fe, conn_data, e, boffset)` that takes its node count from the reference element. The host-side three-argument accessor is unchanged and still serves callers that slice a block outside a kernel. Verified on gfx1102 against Carina's torsion benchmarks: the explicit and quasi-static GPU paths compile and run again, physics bit-identical, and the cost relative to FEC main is within 1.5% on CPU and 0.5% on GPU. Signed-off-by: Alejandro Mota --- src/Enzyme.jl | 4 ++-- src/assemblers/Assemblers.jl | 16 ++++++++++++++-- src/assemblers/MatrixAction.jl | 6 +++--- 3 files changed, 19 insertions(+), 7 deletions(-) diff --git a/src/Enzyme.jl b/src/Enzyme.jl index e6cc8a8e..651fbe3b 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, e, b) + props_el = properties(props, physics, 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, e, b) + props_el = properties(props, physics, 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/assemblers/Assemblers.jl b/src/assemblers/Assemblers.jl index 9a65faf8..dfd27545 100644 --- a/src/assemblers/Assemblers.jl +++ b/src/assemblers/Assemblers.jl @@ -190,6 +190,18 @@ end """ $(TYPEDSIGNATURES) """ +# GPU safe: the property count comes from the physics type, so the slice is a +# statically sized SVector. `view(field.data, range)` builds a SubArray whose +# construction lowers to a dynamic call inside a device kernel. +@inline function properties( + field::PropertyField, ::AbstractPhysics{NF, NP, NS}, e::Int, b::Int +) where {NF, NP, NS} + 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)) +end + @inline function _element_level_properties(props::AbstractArray, ::Int) return props end @@ -419,7 +431,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, e, b) + props_el = properties(props, physics, 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) @@ -451,7 +463,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, e, b) + props_el = properties(props, physics, 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 b7109a22..5dc1ea7f 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, e, b) + props_el = properties(props, physics, 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) @@ -232,7 +232,7 @@ function _assemble_block_matrix_action!( 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 = _element_level_properties(props, e) - props_el = properties(props, e, b) + props_el = properties(props, physics, e, b) K_el = _element_scratch(AssembledMatrix(), ref_fe, U) for q in 1:num_cell_quadrature_points(ref_fe) interps = _cell_interpolants(ref_fe, q) @@ -266,7 +266,7 @@ function _assemble_block_matrix_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, e, b) + props_el = properties(props, physics, 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)