diff --git a/Project.toml b/Project.toml index b7143cd..dc433bc 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 e825895..651fbe3 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, physics, 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, 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/Fields.jl b/src/Fields.jl index 3e02c5e..a5a0eaf 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 5e79a34..479716a 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 e3ebcb4..b8f5d6c 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 6ccf4d0..257e90d 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 92e0bbf..dfd2754 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 @@ -403,6 +415,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 +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 = _element_level_properties(props, e) + 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) @@ -434,6 +447,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 +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 = _element_level_properties(props, e) + 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/Diagonal.jl b/src/assemblers/Diagonal.jl index 91649d2..d3ecbb2 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 b37834a..49d38a4 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 dc67034..b86e86c 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 0ac446a..5dc1ea7 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, 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) @@ -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, 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) @@ -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, 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/QuadratureQuantity.jl b/src/assemblers/QuadratureQuantity.jl index c865512..d8b0420 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 1913dee..443ec30 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 5eb526a..50fadb7 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 e1be447..4fe8a9b 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 23b9a32..6b10fee 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 e686884..8f64f4d 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(