diff --git a/Project.toml b/Project.toml index 0581ffc..c0afd4b 100644 --- a/Project.toml +++ b/Project.toml @@ -23,8 +23,8 @@ Tensors = "48a634ad-e948-5137-8d70-aa71f2a747f4" TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" [weakdeps] -AMDGPU = "21141c5a-9bdb-4563-92ae-f87d6854732e" AbaqusReader = "bc6b9049-e460-56d6-94b4-a597b2c0390d" +AMDGPU = "21141c5a-9bdb-4563-92ae-f87d6854732e" BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e" CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" Gmsh = "705231aa-382f-11e9-3f0c-b7cb4346fdeb" diff --git a/src/Fields.jl b/src/Fields.jl index 1cd0627..24a82e9 100644 --- a/src/Fields.jl +++ b/src/Fields.jl @@ -101,20 +101,30 @@ struct Connectivity{ nepes::Vector{T} nelems::Vector{T} offsets::Vector{T} -end +# end + + function Connectivity(mats::Vector{<:AbstractArray{<:Integer, 2}}) + nblocks = length(mats) + nepes = map(x -> size(x, 1), mats) + nelems = map(x -> size(x, 2), mats) + offsets = Vector{eltype(nepes)}(undef, 0) + offset = 1 + for (nepe, nelem) in zip(nepes, nelems) + push!(offsets, offset) + offset += nepe * nelem + end + data = mapreduce(vec, vcat, mats) + # return Connectivity(data, nblocks, nepes, nelems, offsets) + new{eltype(data), typeof(data)}(data, nblocks, nepes, nelems, offsets) + end -function Connectivity(mats::Vector{<:AbstractArray{<:Integer, 2}}) - nblocks = length(mats) - nepes = map(x -> size(x, 1), mats) - nelems = map(x -> size(x, 2), mats) - offsets = Vector{eltype(nepes)}(undef, 0) - offset = 1 - for (nepe, nelem) in zip(nepes, nelems) - push!(offsets, offset) - offset += nepe * nelem + function Connectivity(data, nblocks, nepes, nelems, offsets) + new{eltype(data), typeof(data)}(data, nblocks, nepes, nelems, offsets) + end + + function Connectivity{T, D}() where {T, D} + new{T, D}(T[], 0, T[], T[], T[]) end - data = mapreduce(vec, vcat, mats) - return Connectivity(data, nblocks, nepes, nelems, offsets) end function Adapt.adapt_structure(to, conn::Connectivity{T, D}) where {T, D} diff --git a/src/Parameters.jl b/src/Parameters.jl index 9bc53f3..189fafd 100644 --- a/src/Parameters.jl +++ b/src/Parameters.jl @@ -11,14 +11,17 @@ $(TYPEDSIGNATURES) $(TYPEDFIELDS) """ struct Parameters{ + IT <: Integer, RT <: Number, - IV <: AbstractVector{<:Integer}, + IV <: AbstractVector{IT}, RV <: AbstractVector{RT}, - RM <: AbstractMatrix, + RM1 <: AbstractMatrix, + RM2 <: AbstractMatrix, ICFuncs <: AbstractVector, DBCFuncs <: AbstractVector, SrcFuncs <: AbstractVector, - NBCs, + # NBCs, + NBCFuncs <: AbstractVector, RBCs, Phys, Props, @@ -27,9 +30,10 @@ struct Parameters{ } <: AbstractParameters ics::InitialConditions{ICFuncs, IV, RV} dirichlet_bcs::DirichletBCs{DBCFuncs, IV, RV} - neumann_bcs::NBCs + # neumann_bcs::NBCs + neumann_bcs::NeumannBCs{NBCFuncs, IT, IV, RM1} robin_bcs::RBCs - sources::Sources{SrcFuncs, RM} + sources::Sources{SrcFuncs, RM2} times::TimeStepper{RT} physics::Phys properties::Props @@ -241,7 +245,7 @@ function update_bc_values!(p::AbstractParameters, assembler) X = coordinates(p) t = current_time(p) update_bc_values!(p.dirichlet_bcs, X, t) - update_bc_values!(p.neumann_bcs, X, t) + update_bc_values!(p.neumann_bcs, assembler, X, t) update_source_values!(p.sources, assembler, X, t) # TODO how to handle Robin BCs? diff --git a/src/assemblers/WeaklyEnforcedBCs.jl b/src/assemblers/WeaklyEnforcedBCs.jl index f96ca01..a8c8bfd 100644 --- a/src/assemblers/WeaklyEnforcedBCs.jl +++ b/src/assemblers/WeaklyEnforcedBCs.jl @@ -8,7 +8,7 @@ function assemble_vector_neumann_bc!( # been conducted previously? _update_for_assembly!(p, assembler.dof, Uu) assemble_vector_weakly_enforced_bc!( - assembler.residual_storage, + assembler.residual_storage, assembler.dof, p.field, coordinates(p), p.neumann_bcs ) return nothing @@ -25,7 +25,7 @@ function assemble_vector_robin_bc!( _update_for_assembly!(p, assembler.dof, Uu) update_bc_values!(p.robin_bcs, coordinates(p), current_time(p), p.field) assemble_vector_weakly_enforced_bc!( - assembler.residual_storage, + assembler.residual_storage, assembler.dof, p.field, coordinates(p), p.neumann_bcs ) return nothing @@ -36,14 +36,16 @@ end $(TYPEDSIGNATURES) """ function assemble_vector_weakly_enforced_bc!( - storage, U, X, bcs + storage, dof, U, X, bcs ) # do not zero! - for bc in values(bcs.bc_caches) + for (n, bc) in enumerate(values(bcs.bc_caches)) + block_id = bcs.block_ids[n] + ref_fe = block_reference_element(function_space(dof), block_id) _assemble_block_vector_weakly_enforced_bc!( storage, U, X, - bc.element_conns.data, bc.ref_fe, bc.sides, bc.vals + bc.element_conns.data, ref_fe, bc.sides, bc.vals ) end end diff --git a/src/bcs/BoundaryConditions.jl b/src/bcs/BoundaryConditions.jl index d541f04..f434899 100644 --- a/src/bcs/BoundaryConditions.jl +++ b/src/bcs/BoundaryConditions.jl @@ -175,20 +175,9 @@ $(TYPEDFIELDS) abstract type AbstractWeaklyEnforcedBCContainer{ IT <: Integer, IV <: AbstractArray{IT, 1}, - RV <: AbstractArray{<:Union{<:Number, <:SVector}, 2}, - RE <: ReferenceFE + RV <: AbstractArray{<:Union{<:Number, <:SVector}, 2} } <: AbstractBCContainer{IV, RV} end -function Adapt.adapt_structure(to, bc::AbstractWeaklyEnforcedBCContainer) - el_conns = adapt(to, bc.element_conns) - elements = adapt(to, bc.elements) - sides = adapt(to, bc.sides) - ref_fe = adapt(to, bc.ref_fe) - vals = adapt(to, bc.vals) - type = eval(typeof(bc).name.name) - return type(el_conns, elements, sides, ref_fe, vals) -end - Base.length(bc::AbstractWeaklyEnforcedBCContainer) = size(bc.vals, 2) function Base.show(io::IO, bc::AbstractWeaklyEnforcedBCContainer) @@ -231,15 +220,6 @@ end """ function update_bc_values! end - -function Adapt.adapt_structure(to, bcs::AbstractBCs) - type = typeof(bcs).name.name - return eval(type)( - map(x -> adapt(to, x), bcs.bc_caches), - adapt(to, bcs.bc_funcs) - ) -end - Base.length(bcs::AbstractBCs) = length(bcs.bc_funcs) # function Base.show(io::IO, bcs::AbstractBCs) @@ -273,6 +253,10 @@ function _setup_weakly_enforced_bc_container(mesh, dof, bcs, type) fspace = function_space(dof) new_bcs = type[] new_funcs = Function[] + block_ids = Int[] + block_names = String[] + # sideset_ids = Int[] + # sideset_names = String[] for (bk, func) in zip(bks, funcs) blocks = sort(unique(bk.blocks)) @@ -283,6 +267,10 @@ function _setup_weakly_enforced_bc_container(mesh, dof, bcs, type) for block in blocks block_name = mesh.element_block_names_map[block] + block_id = findfirst(x -> x == block_name, mesh.element_block_names) + push!(block_ids, block_id) + push!(block_names, block_name) + ids = findall(x -> x == block, bk.blocks) new_blocks = bk.blocks[ids] new_elements = bk.elements[ids] @@ -311,13 +299,13 @@ function _setup_weakly_enforced_bc_container(mesh, dof, bcs, type) if type <: NeumannBCContainer new_bc = type( - conns, new_bk.elements, new_bk.sides, ref_fe, vals + conns, new_bk.elements, new_bk.sides, vals ) elseif type <: RobinBCContainer # dvalsdu = copy(vals) dvalsdu = zeros(SMatrix{ND, ND, Float64, ND * ND}, NQ, length(bk.sides)) new_bc = type( - conns, new_bk.elements, new_bk.sides, ref_fe, vals, dvalsdu + conns, new_bk.elements, new_bk.sides, vals, dvalsdu ) else _unsupported_weak_bc_error("Unsupported weak bc type $type encountered in _setup_weakly_enforced_bc_container") @@ -326,5 +314,5 @@ function _setup_weakly_enforced_bc_container(mesh, dof, bcs, type) push!(new_funcs, func) end end - return new_bcs, new_funcs + return new_bcs, new_funcs, block_ids, block_names end diff --git a/src/bcs/DirichletBCs.jl b/src/bcs/DirichletBCs.jl index 85c9601..9136eff 100644 --- a/src/bcs/DirichletBCs.jl +++ b/src/bcs/DirichletBCs.jl @@ -338,7 +338,20 @@ struct DirichletBCs{ end function Adapt.adapt_structure(to, bcs::DirichletBCs{BCF, IV, RV}) where {BCF, IV, RV} - bc_cache = adapt(to, bcs.bc_cache) + # NOTE + # below logic is needed due to improper + # adapt mapping for an empty array in julia 1.10/1.11 + # where Vector{T}(undef, 0) gets mappend to Vector{Any} + if length(bcs.bc_lengths) > 0 + bc_cache = adapt(to, bcs.bc_cache) + else + temp_int = adapt(to, zeros(Int, 0)) + temp_floats = adapt(to, zeros(Float64, 0)) + bc_cache = DirichletBCContainer{typeof(temp_int), typeof(temp_float)}( + temp_int, copy(temp_int), + temp_floats, copy(temp_floats), copy(temp_floats) + ) + end return DirichletBCs{BCF, typeof(bc_cache.dofs), typeof(bc_cache.vals)}( bc_cache, bcs.bc_funcs, diff --git a/src/bcs/NeumannBCs.jl b/src/bcs/NeumannBCs.jl index 08102d3..d3a2b78 100644 --- a/src/bcs/NeumannBCs.jl +++ b/src/bcs/NeumannBCs.jl @@ -28,14 +28,28 @@ Internal implementation of dirichlet BCs struct NeumannBCContainer{ IT <: Integer, IV <: AbstractArray{IT, 1}, - RV <: AbstractArray{<:Union{<:Number, <:SVector}, 2}, - RE <: ReferenceFE -} <: AbstractWeaklyEnforcedBCContainer{IT, IV, RV, RE} + RM <: AbstractArray{<:Union{<:Number, <:SVector}, 2} +} <: AbstractWeaklyEnforcedBCContainer{IT, IV, RM} element_conns::Connectivity{IT, IV} elements::IV sides::IV - ref_fe::RE - vals::RV + vals::RM + + function NeumannBCContainer(element_conns::Connectivity{IT, IV}, elements::IV, sides::IV, vals::RM) where {IT, IV, RM} + new{IT, IV, RM}(element_conns, elements, sides, vals) + end + + function NeumannBCContainer{IT, IV, RM}() where {IT, IV, RM} + new{IT, IV, RM}(Connectivity{IT, IV}(), IV(undef, 0), IV(undef, 0), RM(undef, 0)) + end +end + +function Adapt.adapt_structure(to, bc::NeumannBCContainer) + el_conns = adapt(to, bc.element_conns) + elements = adapt(to, bc.elements) + sides = adapt(to, bc.sides) + vals = adapt(to, bc.vals) + return NeumannBCContainer(el_conns, elements, sides, vals) end function _update_bc_values!(vals, func, ref_fe, conns, sides, X, t) @@ -51,12 +65,24 @@ function _update_bc_values!(vals, func, ref_fe, conns, sides, X, t) end end +struct NeumannBCFunction{F} <: AbstractBCFunction{F} + func::F +end + +function (f::NeumannBCFunction)(x, t) + return f.func(x, t) +end + struct NeumannBCs{ - BCCaches <: NamedTuple, - BCFuncs <: NamedTuple + BCFuncs, + IT <: Integer, + IV <: AbstractArray{IT, 1}, + RM <: AbstractArray{<:Union{<:Number, <:SVector}, 2}, } <: AbstractBCs{BCFuncs} - bc_caches::BCCaches + bc_caches::Vector{NeumannBCContainer{IT, IV, RM}} bc_funcs::BCFuncs + block_ids::Vector{Int} + block_names::Vector{String} end # note this method has the potential to make @@ -69,21 +95,50 @@ end # TODO below method also currently likely doesn't # handle blocks correclty function NeumannBCs(mesh, dof::DofManager, neumann_bcs::Vector{NeumannBC}) + bc_caches = NeumannBCContainer{Int, Vector{Int}, Matrix{Float64}}[] if length(neumann_bcs) == 0 - return NeumannBCs(NamedTuple(), NamedTuple()) + return NeumannBCs(bc_caches, NeumannBCFunction[], Int[], String[]) end - new_bcs, new_funcs = _setup_weakly_enforced_bc_container(mesh, dof, neumann_bcs, NeumannBCContainer) + new_bcs, new_funcs, block_ids, block_names = _setup_weakly_enforced_bc_container(mesh, dof, neumann_bcs, NeumannBCContainer) - syms = tuple(map(x -> Symbol("neumann_bc_$x"), 1:length(new_bcs))...) - new_bcs = NamedTuple{syms}(tuple(new_bcs...)) - new_funcs = NamedTuple{syms}(tuple(new_funcs...)) - return NeumannBCs(new_bcs, new_funcs) + # TODO fix this up eventually + new_bcs = convert(Vector{typeof(new_bcs[1])}, new_bcs) + new_funcs = NeumannBCFunction.(new_funcs) + # syms = tuple(map(x -> Symbol("neumann_bc_$x"), 1:length(new_bcs))...) + # new_bcs = NamedTuple{syms}(tuple(new_bcs...)) + # new_funcs = NamedTuple{syms}(tuple(new_funcs...)) + return NeumannBCs(new_bcs, new_funcs, block_ids, block_names) +end + +function Adapt.adapt_structure(to, bcs::NeumannBCs) + # NOTE + # below logic is needed due to improper + # adapt mapping for an empty array in julia 1.10/1.11 + # where Vector{T}(undef, 0) gets mappend to Vector{Any} + if length(bcs.bc_caches) > 0 + bc_caches = map(x -> adapt(to, x), bcs.bc_caches) + else + temp_int = adapt(to, zeros(Int, 0)) + temp_floats = adapt(to, zeros(Float64, 0, 0)) + bc_caches = NeumannBCContainer{Int, typeof(temp_int), typeof(temp_floats)}[] + end + return NeumannBCs( + bc_caches, + bcs.bc_funcs, + bcs.block_ids, + bcs.block_names + ) end -function update_bc_values!(bcs::NeumannBCs, X, t) - for (bc, func) in zip(bcs.bc_caches, bcs.bc_funcs) - _update_bc_values!(bc.vals, func, bc.ref_fe, bc.element_conns.data, bc.sides, X, t) +function update_bc_values!(bcs::NeumannBCs, assembler, X, t) + fspace = function_space(assembler.dof) + for n in axes(bcs.block_ids, 1) + bc = bcs.bc_caches[n] + block_id = bcs.block_ids[n] + func = bcs.bc_funcs[n] + ref_fe = block_reference_element(fspace, block_id) + _update_bc_values!(bc.vals, func, ref_fe, bc.element_conns.data, bc.sides, X, t) end return nothing end diff --git a/src/bcs/RobinBCs.jl b/src/bcs/RobinBCs.jl index 6b3850e..2086bd3 100644 --- a/src/bcs/RobinBCs.jl +++ b/src/bcs/RobinBCs.jl @@ -29,13 +29,11 @@ struct RobinBCContainer{ IT <: Integer, IV <: AbstractArray{IT, 1}, RV <: AbstractArray{<:Union{<:Number, <:SVector}, 2}, - dRV <: AbstractArray{<:Union{<:Number, <:SMatrix}, 2}, - RE <: ReferenceFE -} <: AbstractWeaklyEnforcedBCContainer{IT, IV, RV, RE} + dRV <: AbstractArray{<:Union{<:Number, <:SMatrix}, 2} +} <: AbstractWeaklyEnforcedBCContainer{IT, IV, RV} element_conns::Connectivity{IT, IV} elements::IV sides::IV - ref_fe::RE vals::RV dvalsdu::dRV end @@ -44,11 +42,9 @@ function Adapt.adapt_structure(to, bc::RobinBCContainer) el_conns = adapt(to, bc.element_conns) elements = adapt(to, bc.elements) sides = adapt(to, bc.sides) - ref_fe = adapt(to, bc.ref_fe) vals = adapt(to, bc.vals) dvalsdu = adapt(to, bc.dvalsdu) - type = eval(typeof(bc).name.name) - return type(el_conns, elements, sides, ref_fe, vals, dvalsdu) + return RobinBCContainer(el_conns, elements, sides, vals, dvalsdu) end struct RobinBCFunction{F1, F2} <: AbstractBCFunction{F1} @@ -83,24 +79,43 @@ struct RobinBCs{ } <: AbstractBCs{BCFuncs} bc_caches::BCCaches bc_funcs::BCFuncs + block_ids::Vector{Int} + block_names::Vector{String} end function RobinBCs(mesh, dof::DofManager, robin_bcs::Vector{RobinBC}) if length(robin_bcs) == 0 - return RobinBCs(NamedTuple(), NamedTuple()) + return RobinBCs(NamedTuple(), NamedTuple(), Int[], String[]) end - new_bcs, new_funcs = _setup_weakly_enforced_bc_container(mesh, dof, robin_bcs, RobinBCContainer) + new_bcs, new_funcs, block_ids, block_names = _setup_weakly_enforced_bc_container(mesh, dof, robin_bcs, RobinBCContainer) new_funcs = RobinBCFunction.(new_funcs) syms = tuple(map(x -> Symbol("robin_bc_$x"), 1:length(new_bcs))...) new_bcs = NamedTuple{syms}(tuple(new_bcs...)) new_funcs = NamedTuple{syms}(tuple(new_funcs...)) - return RobinBCs(new_bcs, new_funcs) + return RobinBCs(new_bcs, new_funcs, block_ids, block_names) end -function update_bc_values!(bcs::RobinBCs, X, t, U) - for (bc, func) in zip(bcs.bc_caches, bcs.bc_funcs) - _update_bc_values!(bc.vals, bc.dvalsdu, func, bc.ref_fe, bc.element_conns.data, bc.sides, X, t, U) +function Adapt.adapt_structure(to, bcs::RobinBCs) + return RobinBCs( + map(x -> adapt(to, x), bcs.bc_caches), + bcs.bc_funcs, + bcs.block_ids, + bcs.block_names + ) +end + +function update_bc_values!(bcs::RobinBCs, assembler, X, t, U) + # for (bc, func) in zip(bcs.bc_caches, bcs.bc_funcs) + # _update_bc_values!(bc.vals, bc.dvalsdu, func, bc.ref_fe, bc.element_conns.data, bc.sides, X, t, U) + # end + fspace = function_space(assembler.dof) + for n in axes(bcs.block_ids, 1) + bc = values(bcs.bc_caches)[n] + block_id = bcs.block_ids[n] + func = values(bcs.bc_funcs)[n] + ref_fe = block_reference_element(fspace, block_id) + _update_bc_values!(bc.vals, bc.dvalsdu, func, ref_fe, bc.element_conns.data, bc.sides, X, t, U) end return nothing end