From 3f5333eb8952db9a43bcabf527ec16cade8587c9 Mon Sep 17 00:00:00 2001 From: "Craig M. Hamel" Date: Sun, 28 Sep 2025 19:02:47 -0400 Subject: [PATCH] some more testing. --- ext/FiniteElementContainersAdaptExt.jl | 11 ++- src/Parameters.jl | 6 ++ src/assemblers/NeumannBC.jl | 6 +- src/assemblers/SparseMatrixAssembler.jl | 2 +- src/bcs/BoundaryConditions.jl | 108 +++++++++++++++++------- src/bcs/DirichletBCs.jl | 35 ++++---- src/bcs/NeumannBCs.jl | 25 ++++-- src/ics/InitialConditions.jl | 11 +-- test/TestAssemblers.jl | 1 + test/runtests.jl | 14 +-- 10 files changed, 140 insertions(+), 79 deletions(-) diff --git a/ext/FiniteElementContainersAdaptExt.jl b/ext/FiniteElementContainersAdaptExt.jl index 83d2b11..e8ead92 100644 --- a/ext/FiniteElementContainersAdaptExt.jl +++ b/ext/FiniteElementContainersAdaptExt.jl @@ -60,20 +60,23 @@ function Adapt.adapt_structure(to, bk::FiniteElementContainers.BCBookKeeping{V}) end function Adapt.adapt_structure(to, bc::FiniteElementContainers.DirichletBCContainer) - bk = adapt(to, bc.bookkeeping) + dofs = adapt(to, bc.dofs) + nodes = adapt(to, bc.nodes) vals = adapt(to, bc.vals) vals_dot = adapt(to, bc.vals_dot) vals_dot_dot = adapt(to, bc.vals_dot_dot) - return FiniteElementContainers.DirichletBCContainer(bk, vals, vals_dot, vals_dot_dot) + return FiniteElementContainers.DirichletBCContainer(dofs, nodes, vals, vals_dot, vals_dot_dot) end function Adapt.adapt_structure(to, bc::FiniteElementContainers.NeumannBCContainer) - bk = adapt(to, bc.bookkeeping) el_conns = adapt(to, bc.element_conns) + elements = adapt(to, bc.elements) + side_nodes = adapt(to, bc.side_nodes) + sides = adapt(to, bc.sides) surf_conns = adapt(to, bc.surface_conns) ref_fe = adapt(to, bc.ref_fe) vals = adapt(to, bc.vals) - return FiniteElementContainers.NeumannBCContainer(bk, el_conns, surf_conns, ref_fe, vals) + return FiniteElementContainers.NeumannBCContainer(el_conns, elements, side_nodes, sides, surf_conns, ref_fe, vals) end # DofManagers diff --git a/src/Parameters.jl b/src/Parameters.jl index 661666f..ab7003a 100644 --- a/src/Parameters.jl +++ b/src/Parameters.jl @@ -196,6 +196,12 @@ function create_parameters( return Parameters(mesh, assembler, physics, props, ics, dirichlet_bcs, neumann_bcs, times) end +function initialize!(p::Parameters) + update_ic_values!(p.ics, p.ic_funcs, p.h1_coords) + update_field_ics!(p.h1_field, p.ics) + return nothing +end + """ $(TYPEDSIGNATURES) This method is used to update the stored bc values. diff --git a/src/assemblers/NeumannBC.jl b/src/assemblers/NeumannBC.jl index 266a818..0d26ac0 100644 --- a/src/assemblers/NeumannBC.jl +++ b/src/assemblers/NeumannBC.jl @@ -39,7 +39,7 @@ function _assemble_block_vector_neumann_bc!( for e in axes(conns, 2) x_el = _element_level_fields(X, ref_fe, conns, e) R_el = _element_scratch_vector(surface_element(ref_fe.element), U) - side = bc.bookkeeping.sides[e] + side = bc.sides[e] for q in 1:num_quadrature_points(surface_element(ref_fe.element)) interps = MappedSurfaceInterpolants(ref_fe, x_el, q, side) Nvec = interps.N_reduced @@ -54,7 +54,7 @@ function _assemble_block_vector_neumann_bc!( R_el = R_el + JxW * Nvec * f_val end block_id = 1 # doesn't matter for this method - el_id = bc.bookkeeping.elements[e] + el_id = bc.elements[e] @views _assemble_element!(field, R_el, bc.surface_conns[:, e], el_id, block_id) end end @@ -83,7 +83,7 @@ KA.@kernel function _assemble_block_vector_neumann_bc_kernel!( x_el = _element_level_fields(X, ref_fe, conns, E) R_el = _element_scratch_vector(surface_element(ref_fe.element), U) - side = bc.bookkeeping.sides[E] + side = bc.sides[E] for q in 1:num_quadrature_points(surface_element(ref_fe.element)) interps = MappedSurfaceInterpolants(ref_fe, x_el, q, side) diff --git a/src/assemblers/SparseMatrixAssembler.jl b/src/assemblers/SparseMatrixAssembler.jl index c5a0bc2..76b32dd 100644 --- a/src/assemblers/SparseMatrixAssembler.jl +++ b/src/assemblers/SparseMatrixAssembler.jl @@ -228,7 +228,7 @@ function update_dofs!(assembler::SparseMatrixAssembler, dirichlet_bcs) use_condensed = _is_condensed(assembler.dof) if length(dirichlet_bcs) > 0 - dirichlet_dofs = mapreduce(x -> x.bookkeeping.dofs, vcat, dirichlet_bcs) + dirichlet_dofs = mapreduce(x -> x.dofs, vcat, dirichlet_bcs) dirichlet_dofs = unique(sort(dirichlet_dofs)) else dirichlet_dofs = Vector{Int}(undef, 0) diff --git a/src/bcs/BoundaryConditions.jl b/src/bcs/BoundaryConditions.jl index c6a14e5..8bd867d 100644 --- a/src/bcs/BoundaryConditions.jl +++ b/src/bcs/BoundaryConditions.jl @@ -28,6 +28,8 @@ function _unique_sort_perm(array::AbstractArray{T, 1}) where T <: Number return [id_map[x] for x in sorted_unique] end +const SetName = Union{Nothing, Symbol} + """ $(TYPEDEF) $(TYPEDSIGNATURES) @@ -55,36 +57,81 @@ end $(TYPEDSIGNATURES) """ function BCBookKeeping( - mesh, dof::DofManager, var_name::Symbol, sset_name::Symbol + mesh, dof::DofManager, var_name::Symbol; #sset_name::Symbol + block_name::SetName = nothing, + nset_name::SetName = nothing, + sset_name::SetName = nothing ) + # check to ensure at least one name is supplied + if block_name === nothing && + nset_name === nothing && + sset_name === nothing + @assert false "Need to specify either a block, nodeset or sideset." + end + + # now check to make sure only one name is not nothing + if block_name !== nothing + @assert nset_name === nothing && sset_name === nothing + elseif nset_name !== nothing + @assert block_name === nothing && sset_name === nothing + elseif sset_name !== nothing + @assert block_name === nothing && nset_name === nothing + end + # get dof index associated with this var dof_index = _dof_index_from_var_name(dof, var_name) - fspace = function_space(dof) # get sset specific fields - elements = getproperty(mesh.sideset_elems, sset_name) - nodes = getproperty(mesh.sideset_nodes, sset_name) - sides = getproperty(mesh.sideset_sides, sset_name) - side_nodes = getproperty(mesh.sideset_side_nodes, sset_name) - - blocks = Vector{Int64}(undef, 0) - - # gather the blocks that are present in this sideset - # TODO this isn't quite right - for (n, val) in enumerate(values(mesh.element_id_maps)) - # note these are the local elem id to the block, e.g. starting from 1. - indices_in_sset = indexin(val, elements) - filter!(x -> x !== nothing, indices_in_sset) - - if length(indices_in_sset) > 0 - append!(blocks, repeat([n], length(indices_in_sset))) + all_dofs = reshape(1:length(dof), size(dof)) + + if block_name !== nothing + # for this case it is likely a DirichletBC or InitialCondition + # so we really only need the nodes/dofs although this might + # not be the case for say Hdiv or Hcurl fields... + # TODO eventually set the blocks, could be useful maybe? + blocks = Vector{Int64}(undef, 0) + conns = getproperty(mesh.element_conns, block_name) + nodes = sort(unique(conns.data)) + dofs = all_dofs[dof_index, nodes] + elements = getproperty(mesh.element_id_maps, block_name) + # below 2 don't make sense for other mesh entity types + sides = Vector{Int64}(undef, 0) + side_nodes = Matrix{Int64}(undef, 0, 0) + elseif nset_name !== nothing + # for this case we only setup the "nodes" and "dofs" + blocks = Vector{Int64}(undef, 0) # TODO we could eventually put the blocks present here + nodes = getproperty(mesh.nodeset_nodes, nset_name) + dofs = all_dofs[dof_index, nodes] + elements = Vector{Int64}(undef, 0) # TODO we could eventually put the elements present here + # below 2 don't make sense for other mesh entity types + sides = Vector{Int64}(undef, 0) + side_nodes = Matrix{Int64}(undef, 0, 0) + elseif sset_name !== nothing + + elements = getproperty(mesh.sideset_elems, sset_name) + nodes = getproperty(mesh.sideset_nodes, sset_name) + sides = getproperty(mesh.sideset_sides, sset_name) + side_nodes = getproperty(mesh.sideset_side_nodes, sset_name) + + blocks = Vector{Int64}(undef, 0) + + # gather the blocks that are present in this sideset + # TODO this isn't quite right + for (n, val) in enumerate(values(mesh.element_id_maps)) + # note these are the local elem id to the block, e.g. starting from 1. + indices_in_sset = indexin(val, elements) + filter!(x -> x !== nothing, indices_in_sset) + + if length(indices_in_sset) > 0 + append!(blocks, repeat([n], length(indices_in_sset))) + end end - end - # setup dofs local to this BC - all_dofs = reshape(1:length(dof), size(dof)) - dofs = all_dofs[dof_index, nodes] + # setup dofs local to this BC + # all_dofs = reshape(1:length(dof), size(dof)) + dofs = all_dofs[dof_index, nodes] + end return BCBookKeeping( blocks, dofs, elements, nodes, sides, side_nodes @@ -143,14 +190,15 @@ $(TYPEDEF) $(TYPEDSIGNATURES) $(TYPEDFIELDS) """ -abstract type AbstractBCContainer{ - IT <: Integer, - VT <: Union{<:Number, <:SVector}, - N, - IV <: AbstractArray{IT, 1}, - IM <: AbstractArray{IT, 2}, - VV <: AbstractArray{VT, N} -} end +# abstract type AbstractBCContainer{ +# IT <: Integer, +# VT <: Union{<:Number, <:SVector}, +# N, +# IV <: AbstractArray{IT, 1}, +# IM <: AbstractArray{IT, 2}, +# VV <: AbstractArray{VT, N} +# } end +abstract type AbstractBCContainer end KA.get_backend(x::AbstractBCContainer) = KA.get_backend(x.vals) diff --git a/src/bcs/DirichletBCs.jl b/src/bcs/DirichletBCs.jl index 8f29853..e7232cc 100644 --- a/src/bcs/DirichletBCs.jl +++ b/src/bcs/DirichletBCs.jl @@ -1,7 +1,7 @@ abstract type AbstractDirichletBC{F} <: AbstractBC{F} end -abstract type AbstractDirichletBCContainer{ - IT, VT, IV, IM, VV -} <: AbstractBCContainer{IT, VT, 1, IV, IM, VV} end +# abstract type AbstractDirichletBCContainer{ +# IT, VT, IV, IM, VV +# } <: AbstractBCContainer{IT, VT, 1, IV, IM, VV} end """ $(TYPEDEF) @@ -43,10 +43,10 @@ struct DirichletBCContainer{ IT <: Integer, VT <: Number, IV <: AbstractArray{IT, 1}, - IM <: AbstractArray{IT, 2}, VV <: AbstractArray{VT, 1} -} <: AbstractDirichletBCContainer{IT, VT, IV, IM, VV} - bookkeeping::BCBookKeeping{IT, IV, IM} +} <: AbstractBCContainer + dofs::IV + nodes::IV vals::VV vals_dot::VV vals_dot_dot::VV @@ -58,7 +58,7 @@ $(TYPEDSIGNATURES) $(TYPEDFIELDS) """ function DirichletBCContainer(mesh, dof::DofManager, dbc::DirichletBC) - bk = BCBookKeeping(mesh, dof, dbc.var_name, dbc.sset_name) + bk = BCBookKeeping(mesh, dof, dbc.var_name, sset_name=dbc.sset_name) # sort nodes and dofs for dirichlet bc dof_perm = _unique_sort_perm(bk.dofs) @@ -72,11 +72,11 @@ function DirichletBCContainer(mesh, dof::DofManager, dbc::DirichletBC) vals = zeros(length(bk.nodes)) vals_dot = zeros(length(bk.nodes)) vals_dot_dot = zeros(length(bk.nodes)) - return DirichletBCContainer(bk, vals, vals_dot, vals_dot_dot) + return DirichletBCContainer(bk.dofs, bk.nodes, vals, vals_dot, vals_dot_dot) end function Base.length(bc::DirichletBCContainer) - return length(bc.bookkeeping.dofs) + return length(bc.dofs) end # need checks on if field types are compatable @@ -87,7 +87,7 @@ based on the stored function """ function _update_bc_values!(bc::DirichletBCContainer, func, X, t, ::KA.CPU) ND = num_fields(X) - for (n, node) in enumerate(bc.bookkeeping.nodes) + for (n, node) in enumerate(bc.nodes) X_temp = @views SVector{ND, eltype(X)}(X[:, node]) bc.vals[n] = func.func(X_temp, t) bc.vals_dot[n] = func.func_dot(X_temp, t) @@ -104,7 +104,7 @@ GPU kernel for updating stored bc values based on the stored function KA.@kernel function _update_bc_values_kernel!(bc::DirichletBCContainer, func, X, t) I = KA.@index(Global) ND = num_fields(X) - node = bc.bookkeeping.nodes[I] + node = bc.nodes[I] # hacky for now, but it works # can't do X[:, node] on the GPU, this results in a dynamic @@ -134,14 +134,14 @@ end # TODO change below names to be more specific to dbcs function _update_field_dirichlet_bcs!(U, bc::DirichletBCContainer, ::KA.CPU) - for (dof, val) in zip(bc.bookkeeping.dofs, bc.vals) + for (dof, val) in zip(bc.dofs, bc.vals) U[dof] = val end return nothing end function _update_field_dirichlet_bcs!(U, V, bc::DirichletBCContainer, ::KA.CPU) - for (dof, val, val_dot) in zip(bc.bookkeeping.dofs, bc.vals, bc.vals_dot) + for (dof, val, val_dot) in zip(bc.dofs, bc.vals, bc.vals_dot) U[dof] = val V[dof] = val_dot end @@ -149,7 +149,7 @@ function _update_field_dirichlet_bcs!(U, V, bc::DirichletBCContainer, ::KA.CPU) end function _update_field_dirichlet_bcs!(U, V, A, bc::DirichletBCContainer, ::KA.CPU) - for (dof, val, val_dot, val_dot_dot) in zip(bc.bookkeeping.dofs, bc.vals, bc.vals_dot, bc.vals_dot_dot) + for (dof, val, val_dot, val_dot_dot) in zip(bc.dofs, bc.vals, bc.vals_dot, bc.vals_dot_dot) U[dof] = val V[dof] = val_dot A[dof] = val_dot_dot @@ -160,8 +160,7 @@ end # COV_EXCL_START KA.@kernel function _update_field_dirichlet_bcs_kernel!(U, bc::DirichletBCContainer) I = KA.@index(Global) - dof = bc.bookkeeping.dofs[I] - # val = bc.vals[I] + dof = bc.dofs[I] U[dof] = bc.vals[I] end # COV_EXCL_STOP @@ -169,7 +168,7 @@ end # COV_EXCL_START KA.@kernel function _update_field_dirichlet_bcs_kernel!(U, V, A, bc::DirichletBCContainer) I = KA.@index(Global) - dof = bc.bookkeeping.dofs[I] + dof = bc.dofs[I] val = bc.vals[I] U[dof] = bc.vals[I] V[dof] = bc.vals_dot[I] @@ -232,7 +231,7 @@ function create_dirichlet_bcs(mesh, dof::DofManager, dirichlet_bcs::Vector{<:Dir dirichlet_bcs = DirichletBCContainer.((mesh,), (dof,), dirichlet_bcs) if length(dirichlet_bcs) > 0 - temp_dofs = mapreduce(x -> x.bookkeeping.dofs, vcat, dirichlet_bcs) + temp_dofs = mapreduce(x -> x.dofs, vcat, dirichlet_bcs) temp_dofs = unique(sort(temp_dofs)) update_dofs!(dof, temp_dofs) end diff --git a/src/bcs/NeumannBCs.jl b/src/bcs/NeumannBCs.jl index ad83a32..1c3d2e3 100644 --- a/src/bcs/NeumannBCs.jl +++ b/src/bcs/NeumannBCs.jl @@ -44,9 +44,12 @@ struct NeumannBCContainer{ C1, # TODO specialize C2, # TODO specialize RE <: ReferenceFE -} <: AbstractBCContainer{IT, VT, 2, IV, IM, VV} - bookkeeping::BCBookKeeping{IT, IV, IM} +} <: AbstractBCContainer + # bookkeeping::BCBookKeeping{IT, IV, IM} element_conns::C1 + elements::IV + side_nodes::IM + sides::IV surface_conns::C2 ref_fe::RE vals::VV @@ -56,13 +59,13 @@ function _update_bc_values!(bc::NeumannBCContainer, func, X, t, ::KA.CPU) ND = size(X, 1) NN = num_vertices(bc.ref_fe) NNPS = num_vertices(surface_element(bc.ref_fe.element)) - for (n, e) in enumerate(bc.bookkeeping.elements) + for (n, e) in enumerate(bc.elements) conn = @views bc.element_conns[:, n] X_el = SVector{ND * NN, eltype(X)}(@views X[:, conn]) X_el = SMatrix{length(X_el) ÷ ND, ND, eltype(X_el), length(X_el)}(X_el...) for q in 1:num_quadrature_points(surface_element(bc.ref_fe.element)) - side = bc.bookkeeping.sides[n] + side = bc.sides[n] interps = MappedSurfaceInterpolants(bc.ref_fe, X_el, q, side) X_q = interps.X_q bc.vals[q, n] = func(X_q, t) @@ -77,14 +80,14 @@ KA.@kernel function _update_bc_values_kernel!(bc::NeumannBCContainer, func, X, t Q, E = KA.@index(Global, NTuple) # E = KA.@index(Global) - el_id = bc.bookkeeping.elements[E] + el_id = bc.elements[E] conn = @views bc.element_conns[:, E] X_el = SVector{ND * NN, eltype(X)}(@views X[:, conn]) X_el = SMatrix{length(X_el) ÷ ND, ND, eltype(X_el), length(X_el)}(X_el...) # for q in 1:num_quadrature_points(bc.ref_fe.surface_element) - side = bc.bookkeeping.sides[E] + side = bc.sides[E] interps = MappedSurfaceInterpolants(bc.ref_fe, X_el, Q, side) X_q = interps.X_q bc.vals[Q, E] = func(X_q, t) @@ -109,8 +112,11 @@ function create_neumann_bcs(mesh, dof::DofManager, neumann_bcs::Vector{NeumannBC sets = map(x -> x.sset_name, neumann_bcs) vars = map(x -> x.var_name, neumann_bcs) funcs = map(x -> x.func, neumann_bcs) - bks = BCBookKeeping.((mesh,), (dof,), vars, sets) - # bks = _split_bookkeeping_by_block(bks) + + # NOTE neumann bcs must be present on a sideset + # so that is the only mesh entity that will be + # supported for this BC type + bks = map((v, s) -> BCBookKeeping(mesh, dof, v; sset_name=s), vars, sets) fspace = function_space(dof) new_bcs = NeumannBCContainer[] new_funcs = Function[] @@ -178,7 +184,8 @@ function create_neumann_bcs(mesh, dof::DofManager, neumann_bcs::Vector{NeumannBC # typeof(new_bk), typeof(conns), typeof(surface_conns), typeof(ref_fe), eltype(vals), typeof(vals) # }( new_bc = NeumannBCContainer( - new_bk, conns, surface_conns, ref_fe, vals + # new_bk, conns, surface_conns, ref_fe, vals + conns, new_bk.elements, new_bk.side_nodes, new_bk.sides, surface_conns, ref_fe, vals ) push!(new_bcs, new_bc) push!(new_funcs, func) diff --git a/src/ics/InitialConditions.jl b/src/ics/InitialConditions.jl index 1d68c14..0493766 100644 --- a/src/ics/InitialConditions.jl +++ b/src/ics/InitialConditions.jl @@ -32,14 +32,9 @@ end # to "do the right thing" depending upon the field type # this is set up function InitialConditionContainer(mesh, dof, ic::InitialCondition) - block_conns = getproperty(mesh.element_conns, ic.block_name) - block_nodes = unique(sort(block_conns.data)) - dof_index = _dof_index_from_var_name(dof, ic.var_name) - - all_dofs = reshape(1:length(dof), size(dof)) - block_dofs = all_dofs[dof_index, block_nodes] - vals = zeros(length(block_dofs)) - return InitialConditionContainer(block_dofs, block_nodes, vals) + bk = BCBookKeeping(mesh, dof, ic.var_name; block_name=ic.block_name) + vals = zeros(length(bk.dofs)) + return InitialConditionContainer(bk.dofs, bk.nodes, vals) end function Base.length(ic::InitialConditionContainer) diff --git a/test/TestAssemblers.jl b/test/TestAssemblers.jl index dfbc4c6..14df7a9 100644 --- a/test/TestAssemblers.jl +++ b/test/TestAssemblers.jl @@ -21,6 +21,7 @@ include("poisson/TestPoissonCommon.jl") DirichletBC(:u, :sset_4, bc_func), ] p = create_parameters(mesh, asm, physics, props; dirichlet_bcs=dbcs) + FiniteElementContainers.initialize!(p) Uu = create_unknowns(asm) Vu = create_unknowns(asm) diff --git a/test/runtests.jl b/test/runtests.jl index 9f1a231..019b3ad 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -56,12 +56,14 @@ include("TestPhysics.jl") end end - # if CUDA.functional() - # test_poisson_dirichlet(cuda, false, NewtonSolver, cg_solver) - # test_poisson_dirichlet(cuda, true, NewtonSolver, cg_solver) - # test_poisson_neumann(cuda, false, NewtonSolver, cg_solver) - # test_poisson_neumann(cuda, true, NewtonSolver, cg_solver) - # end + if CUDA.functional() + for cond in condensed + test_poisson_dirichlet(cuda, cond, NewtonSolver, cg_solver) + test_poisson_dirichlet_multi_block_quad4_quad4(cuda, cond, NewtonSolver, cg_solver) + test_poisson_dirichlet_multi_block_quad4_tri3(cuda, cond, NewtonSolver, cg_solver) + test_poisson_neumann(cuda, cond, NewtonSolver, cg_solver) + end + end end @testset ExtendedTestSet "Mechanics Problem" begin