From f138f544604f89102e3635bc90f62e1c07b29214 Mon Sep 17 00:00:00 2001 From: "Craig M. Hamel" Date: Fri, 12 Jun 2026 09:20:47 -0400 Subject: [PATCH 1/4] lots of clean up on periodic bcs. Now working with inplace and sparse vector stuff. --- examples/poisson/periodic_bc_v2.jl | 18 ++-- src/DofManagers.jl | 109 +++++++++++++++++------- src/assemblers/Assemblers.jl | 7 ++ src/assemblers/SparseMatrixAssembler.jl | 3 +- src/bcs/PeriodicBCs.jl | 102 +++++++++++----------- test/poisson/TestPoissonPBCs.jl | 40 +++++++++ 6 files changed, 185 insertions(+), 94 deletions(-) diff --git a/examples/poisson/periodic_bc_v2.jl b/examples/poisson/periodic_bc_v2.jl index 544abc7a..7f598e00 100644 --- a/examples/poisson/periodic_bc_v2.jl +++ b/examples/poisson/periodic_bc_v2.jl @@ -4,12 +4,6 @@ using Krylov using LinearAlgebra using SparseArrays - -# THINGS TO DO -# 1. Need to propogate dof_to_unknown map to assemblers -# so we can use the new dof_to_unknown_index method -# or at least that behavior so that PBC assembly works correctly - include("../../test/poisson/TestPoissonCommon.jl") # function simulation() @@ -32,7 +26,7 @@ physics = Poisson(f) props = create_properties(physics) u = ScalarFunction(V, "u") dof = DofManager(u) -asm = SparseMatrixAssembler(dof) +asm = SparseMatrixAssembler(dof; use_inplace_methods = false, use_sparse_vector = false) pbcs = PeriodicBC[ PeriodicBC("u", "x", zero_func, "sset_1", "sset_3") @@ -40,7 +34,7 @@ pbcs = PeriodicBC[ ] p = create_parameters(mesh, asm, physics, props; periodic_bcs = pbcs) # U = create_unknowns(dof) -solver = NewtonSolver(DirectLinearSolver(asm)) +solver = NewtonSolver(IterativeLinearSolver(asm, :cg)) integrator = QuasiStaticIntegrator(solver) evolve!(integrator, p) @@ -48,11 +42,11 @@ evolve!(integrator, p) U = p.field # U_an = similar(U) # X = p.coords -u_an = u_analytic.(eachcol(mesh.nodal_coords)) +u_an = H1Field{Float64, Vector{Float64}, 1}(u_analytic.(eachcol(mesh.nodal_coords))) -# for n in axes(mesh.nodal_coords, 2) -# @assert isapprox(U[n], u_an[n], atol=5e-2) -# end +for n in axes(mesh.nodal_coords, 2) + @assert isapprox(U[n], u_an[n], atol=5e-4) +end # pp = PostProcessor(mesh, output_file, u; copy_mesh_file = false) diff --git a/src/DofManagers.jl b/src/DofManagers.jl index a05c636b..81f4768a 100644 --- a/src/DofManagers.jl +++ b/src/DofManagers.jl @@ -13,6 +13,46 @@ abstract type AbstractDofManager{ IDs <: AbstractArray{IT, 1} } end +# struct DofManagerCache{ +# IT <: Integer, +# IDs <: AbstractVector{IT} +# } <: AbstractDofManager{IT, IDs} +# dirichlet_dofs::IDs +# dof_to_unknown::IDs +# periodic_side_a_dofs::IDs +# periodic_side_b_dofs::IDs +# periodic_side_b_to_side_a_unknown::IDs +# unknown_dofs::IDs +# end + +# function Adapt.adapt_structure(to, cache::DofManagerCache) +# return DofManagerCache( +# adapt(to, cache.dirichlet_dofs), +# adapt(to, cache.dof_to_unknown), +# adapt(to, cache.periodic_side_a_dofs), +# adapt(to, cache.periodic_side_b_dofs), +# adapt(to, cache.periodic_side_b_to_side_a_unknown), +# adapt(to, cache.unknown_dofs) +# ) +# end + +# Base.eltype(::DofManagerCache{IT, IDs}) where {IT, IDs} = IT +# _ids_type(::DofManagerCache{IT, IDs}) where {IT, IDs} = IDs + +# @inline function dof_to_unknown_index(dof::DofManagerCache, n::Int) +# dtu = dof.dof_to_unknown[n] +# if dtu == DIRICHLET_DOF +# return DIRICHLET_DOF +# elseif dtu == PERIODIC_SIDE_B_DOF +# side_a_unknown = dof.periodic_side_b_to_side_a_unknown[n] +# # @assert side_a_unknown != PERIODIC_SIDE_B_DOF +# # TODO need some kind of check here +# return side_a_unknown +# else +# return dtu +# end +# end + """ $(TYPEDEF) $(TYPEDSIGNATURES) @@ -52,35 +92,40 @@ struct DofManager{ return DofManager{Condensed}( dirichlet_dofs, dof_to_unknown, periodic_side_a_dofs, periodic_side_b_dofs, - periodic_side_b_to_side_a_unknown, unknown_dofs, var + periodic_side_b_to_side_a_unknown, + unknown_dofs, var ) end function DofManager{Condensed}( - dirichlet_dofs, dof_to_unknown, + dirichlet_dofs::IDs, dof_to_unknown, periodic_side_a_dofs, periodic_side_b_dofs, - periodic_side_b_to_side_a_unknown, unknown_dofs, var - ) where Condensed - new{Condensed, eltype(dirichlet_dofs), typeof(dirichlet_dofs), typeof(var)}( + periodic_side_b_to_side_a_unknown, + unknown_dofs, var::V + ) where {Condensed, IDs, V} + new{Condensed, eltype(dirichlet_dofs), IDs, V}( dirichlet_dofs, dof_to_unknown, - periodic_side_a_dofs, periodic_side_b_dofs, - periodic_side_b_to_side_a_unknown, unknown_dofs, var + periodic_side_a_dofs, periodic_side_b_dofs, + periodic_side_b_to_side_a_unknown, + unknown_dofs, var ) end end function Adapt.adapt_structure(to, dof::DofManager) return DofManager{_is_condensed(dof)}( - adapt(to, dof.dirichlet_dofs), - adapt(to, dof.dof_to_unknown), - adapt(to, dof.periodic_side_a_dofs), - adapt(to, dof.periodic_side_b_dofs), - adapt(to, dof.periodic_side_b_to_side_a_unknown), - adapt(to, dof.unknown_dofs), + adapt(to, dof), + adapt(to, cache.dirichlet_dofs), + adapt(to, cache.dof_to_unknown), + adapt(to, cache.periodic_side_a_dofs), + adapt(to, cache.periodic_side_b_dofs), + adapt(to, cache.periodic_side_b_to_side_a_unknown), + adapt(to, cache.unknown_dofs), adapt(to, dof.var) ) end +_ids_type(::DofManager{C, IT, IDs, V}) where {C, IT, IDs, V} = IDs _is_condensed(::DofManager{C, IT, IDs, V}) where {C, IT, IDs, V} = C """ @@ -90,7 +135,10 @@ Return the total lenth of dofs in the problem. E.g. for an H1 space this will the number of nodes times the number of degrees of freedom per node. """ -Base.length(dof::DofManager) = length(dof.dirichlet_dofs) + length(dof.periodic_side_b_dofs) + length(dof.unknown_dofs) +function Base.length(dof::DofManager) + cache = dof + return length(cache.dirichlet_dofs) + length(cache.periodic_side_b_dofs) + length(cache.unknown_dofs) +end """ $(TYPEDSIGNATURES) @@ -167,22 +215,6 @@ function create_unknowns(dof::DofManager{Flag, IT, IDs, Var}, float_type = Float end end -function extract_field_unknowns!( - Uu::V, - dof::DofManager, - U::AbstractField -) where V <: AbstractVector{<:Number} - unknown_dofs = dof.unknown_dofs - fec_foreach(Uu) do n - Uu[n] = U[unknown_dofs[n]] - end - return nothing -end - -function function_space(dof::DofManager) - return dof.var.fspace -end - """ Map a single global DOF index to its reduced unknown index. @@ -195,6 +227,7 @@ by global DOF, so this is a pure array lookup with no allocation — safe to call inside GPU kernels. """ @inline function dof_to_unknown_index(dof::DofManager, n::Int) + # return dof_to_unknown_index(dof, n) dtu = dof.dof_to_unknown[n] if dtu == DIRICHLET_DOF return DIRICHLET_DOF @@ -208,6 +241,22 @@ call inside GPU kernels. end end +function extract_field_unknowns!( + Uu::V, + dof::DofManager, + U::AbstractField +) where V <: AbstractVector{<:Number} + unknown_dofs = dof.unknown_dofs + fec_foreach(Uu) do n + Uu[n] = U[unknown_dofs[n]] + end + return nothing +end + +function function_space(dof::DofManager) + return dof.var.fspace +end + """ $(TYPEDSIGNATURES) Takes in a list of dof ids associated with dirichlet bcs diff --git a/src/assemblers/Assemblers.jl b/src/assemblers/Assemblers.jl index 7b1cf62b..dd009d58 100644 --- a/src/assemblers/Assemblers.jl +++ b/src/assemblers/Assemblers.jl @@ -87,6 +87,7 @@ function _assemble_element!( end # sparse vector attempt +# this one is good for pbcs function _assemble_element!( storage::AbstractVector, R_el::SVector{NDOF, T}, conns, @@ -336,6 +337,12 @@ function residual(asm::AbstractAssembler) ) return asm.residual_storage.data else + # need to adjust for PBCs + dof = asm.dof + for (side_a_dof, side_b_dof) in zip(dof.periodic_side_a_dofs, dof.periodic_side_b_dofs) + asm.residual_storage[side_a_dof] += asm.residual_storage[side_b_dof] + end + extract_field_unknowns!( asm.residual_unknowns, asm.dof, diff --git a/src/assemblers/SparseMatrixAssembler.jl b/src/assemblers/SparseMatrixAssembler.jl index d1e845aa..14ff32a9 100644 --- a/src/assemblers/SparseMatrixAssembler.jl +++ b/src/assemblers/SparseMatrixAssembler.jl @@ -40,7 +40,7 @@ struct SparseMatrixAssembler{ ) where {Condensed, SparseMatrixType, UseInPlaceMethods, UseSparseVec} new{ Condensed, SparseMatrixType, UseInPlaceMethods, UseSparseVec, - typeof(dof.unknown_dofs), typeof(residual_storage.data), + _ids_type(dof), typeof(residual_storage.data), typeof(dof.var), typeof(stiffness_action_storage) }( dof, matrix_pattern, vector_pattern, @@ -267,5 +267,6 @@ function update_dofs!(assembler::AbstractAssembler, dirichlet_bcs::DirichletBCs, _, n_entries = _setup_block_sizes(assembler.dof, 1) resize!(assembler.residual_unknowns, n_entries) end + return nothing end diff --git a/src/bcs/PeriodicBCs.jl b/src/bcs/PeriodicBCs.jl index d4c490a0..798b1b6b 100644 --- a/src/bcs/PeriodicBCs.jl +++ b/src/bcs/PeriodicBCs.jl @@ -1,9 +1,3 @@ -# For a periodic pair in direction dir_id, match on the OTHER coordinates -# function _transverse_key(coords, node, dir_id, tolerance) -# ndim = size(coords, 1) -# return tuple([round(Int, coords[d, node] / tolerance) for d in 1:ndim if d != dir_id]...) -# end - struct PeriodicBC{F} <: AbstractBC{F} direction::String func::F @@ -135,7 +129,7 @@ function _update_bc_values!(bc::PeriodicBCContainer, func, X, t) elseif ND == 3 X_temp = SVector{ND, eltype(X)}(X[1, node], X[2, node], X[3, node]) end - bc.vals[n] = func.func(X_temp, t) + bc.vals[n] = func(X_temp, t) end return nothing end @@ -152,6 +146,10 @@ struct PeriodicBCFunction{F} <: AbstractBCFunction{F} end end +function (func::PeriodicBCFunction)(X, t) + return func.func(X, t) +end + struct PeriodicBCs{ BCFuncs, IV <: AbstractVector{<:Integer}, @@ -220,7 +218,53 @@ function Adapt.adapt_structure(to, pbcs::PeriodicBCs) ) end +Base.length(bcs::PeriodicBCs) = length(bcs.bc_funcs) + +function Base.show(io::IO, bcs::PeriodicBCs) + # for (n, (cache, func)) in enumerate(zip(bcs.bc_cache, bcs.bc_funcs)) + show(io, "PeriodicBC:") + show(io, bcs.bc_cache) + # show(io, bcs.bc_lengths) + # show(io, func) + show(io, "\n") + # end +end + +function periodic_dofs(bcs::PeriodicBCs) + if length(bcs) > 0 + side_a_dofs = mapreduce(x -> x.side_a_dofs, vcat, bcs.bc_caches) + side_b_dofs = mapreduce(x -> x.side_b_dofs, vcat, bcs.bc_caches) + perm = _unique_sort_perm(side_b_dofs) # sort by side b, side a will seem unsorted + side_a_dofs = side_a_dofs[perm] + side_b_dofs = side_b_dofs[perm] + else + return Vector{Int}(undef, 0), Vector{Int}(undef, 0) + end + return side_a_dofs, side_b_dofs +end + +function update_bc_values!(bcs::PeriodicBCs, X, t) + for n in axes(bcs.bc_caches, 1) + _update_bc_values!(bcs.bc_caches[n], bcs.bc_funcs[n], X, t) + end + return nothing +end + +function update_field_periodic_bcs!(U, bcs::PeriodicBCs) + for n in axes(bcs.bc_caches, 1) + cache = bcs.bc_caches[n] + fec_foreach(cache.side_b_dofs) do I + side_a_dof = cache.side_a_dofs[I] + side_b_dof = cache.side_b_dofs[I] + U[side_b_dof] = U[side_a_dof] + cache.vals[I] + end + end +end + +##################################################### +# old likely unused below # below is for lagrange multiplier approach +##################################################### function _constraint_matrix(dof::DofManager, pbcs::PeriodicBCs) Is, Js, Vs = Int[], Int[], Float64[] n = 1 @@ -273,47 +317,3 @@ function _create_constraint_field(dof::DofManager, pbcs::PeriodicBCs) n_constraints = mapreduce(x -> length(x.side_a_dofs), +, pbcs.bc_caches) return zeros(n_constraints) end - -Base.length(bcs::PeriodicBCs) = length(bcs.bc_funcs) - -function Base.show(io::IO, bcs::PeriodicBCs) - # for (n, (cache, func)) in enumerate(zip(bcs.bc_cache, bcs.bc_funcs)) - show(io, "PeriodicBC:") - show(io, bcs.bc_cache) - # show(io, bcs.bc_lengths) - # show(io, func) - show(io, "\n") - # end -end - -function periodic_dofs(bcs::PeriodicBCs) - if length(bcs) > 0 - side_a_dofs = mapreduce(x -> x.side_a_dofs, vcat, bcs.bc_caches) - side_b_dofs = mapreduce(x -> x.side_b_dofs, vcat, bcs.bc_caches) - perm = _unique_sort_perm(side_b_dofs) # sort by side b, side a will seem unsorted - side_a_dofs = side_a_dofs[perm] - side_b_dofs = side_b_dofs[perm] - else - return Vector{Int}(undef, 0), Vector{Int}(undef, 0) - end - return side_a_dofs, side_b_dofs -end - -function update_bc_values!(bcs::PeriodicBCs, X, t) - # for bc_cache, bc_func in zbcs.bc_caches - for n in axes(bcs.bc_caches, 1) - _update_bc_values!(bcs.bc_caches[n], bcs.bc_funcs[n], X, t) - end - return nothing -end - -function update_field_periodic_bcs!(U, bcs::PeriodicBCs) - for n in axes(bcs.bc_caches, 1) - cache = bcs.bc_caches[n] - fec_foreach(cache.side_b_dofs) do I - side_a_dof = cache.side_a_dofs[I] - side_b_dof = cache.side_b_dofs[I] - U[side_b_dof] = U[side_a_dof] + cache.vals[I] - end - end -end diff --git a/test/poisson/TestPoissonPBCs.jl b/test/poisson/TestPoissonPBCs.jl index 3ccc2498..f6af2112 100644 --- a/test/poisson/TestPoissonPBCs.jl +++ b/test/poisson/TestPoissonPBCs.jl @@ -78,3 +78,43 @@ end @testitem "Regression test - test_poisson_pbcs" setup=[PoissonPBCHelper] begin test_poisson_pbcs() end + +@testitem "Regression test - test_poisson_periodic" setup=[PoissonRegressionHelper] begin + output_file = "poisson_periodic_test.e" + f_pbc(X, _) = begin + x, y = X + (2π)^2 * cos(2π * x) + + 0.5 * (4π)^2 * cos(4π * y) + + 0.25 * ((2π)^2 + (4π)^2) * sin(2π * x) * sin(4π * y) + end + u_analytic(x) = cos(2π*x[1]) + 0.5*cos(4π*x[2]) + 0.25*sin(2π*x[1])*sin(4π*x[2]) + zero_func(_, _) = 0.0 + + physics_pbcs = Poisson(f_pbc) + for dev in backends + for use_inplace_methods in [false, true] + asm = SparseMatrixAssembler( + u; + sparse_matrix_type = :csc, + use_condensed = false, + use_inplace_methods = use_inplace_methods + ) + + # setup bcs + pbcs = PeriodicBC[ + PeriodicBC("u", "x", zero_func, "sset_1", "sset_3") + PeriodicBC("u", "y", zero_func, "sset_4", "sset_2") + ] + p = create_parameters(mesh, asm, physics_pbcs, props; periodic_bcs = pbcs) + solver = NewtonSolver(IterativeLinearSolver(asm, :cg)) + integrator = QuasiStaticIntegrator(solver) + evolve!(integrator, p) + U = p.field + U_an = H1Field{Float64, Vector{Float64}, 1}(u_analytic.(eachcol(mesh.nodal_coords))) + + for n in axes(mesh.nodal_coords, 2) + @test isapprox(U[n], U_an[n], atol=5e-4) + end + end + end +end \ No newline at end of file From aa57d67fba4158a1d90c2bc9f4beb01226123316 Mon Sep 17 00:00:00 2001 From: "Craig M. Hamel" Date: Fri, 12 Jun 2026 09:22:59 -0400 Subject: [PATCH 2/4] some cleanup in dof managers. --- src/DofManagers.jl | 40 ---------------------------------------- 1 file changed, 40 deletions(-) diff --git a/src/DofManagers.jl b/src/DofManagers.jl index 81f4768a..222f66cb 100644 --- a/src/DofManagers.jl +++ b/src/DofManagers.jl @@ -13,46 +13,6 @@ abstract type AbstractDofManager{ IDs <: AbstractArray{IT, 1} } end -# struct DofManagerCache{ -# IT <: Integer, -# IDs <: AbstractVector{IT} -# } <: AbstractDofManager{IT, IDs} -# dirichlet_dofs::IDs -# dof_to_unknown::IDs -# periodic_side_a_dofs::IDs -# periodic_side_b_dofs::IDs -# periodic_side_b_to_side_a_unknown::IDs -# unknown_dofs::IDs -# end - -# function Adapt.adapt_structure(to, cache::DofManagerCache) -# return DofManagerCache( -# adapt(to, cache.dirichlet_dofs), -# adapt(to, cache.dof_to_unknown), -# adapt(to, cache.periodic_side_a_dofs), -# adapt(to, cache.periodic_side_b_dofs), -# adapt(to, cache.periodic_side_b_to_side_a_unknown), -# adapt(to, cache.unknown_dofs) -# ) -# end - -# Base.eltype(::DofManagerCache{IT, IDs}) where {IT, IDs} = IT -# _ids_type(::DofManagerCache{IT, IDs}) where {IT, IDs} = IDs - -# @inline function dof_to_unknown_index(dof::DofManagerCache, n::Int) -# dtu = dof.dof_to_unknown[n] -# if dtu == DIRICHLET_DOF -# return DIRICHLET_DOF -# elseif dtu == PERIODIC_SIDE_B_DOF -# side_a_unknown = dof.periodic_side_b_to_side_a_unknown[n] -# # @assert side_a_unknown != PERIODIC_SIDE_B_DOF -# # TODO need some kind of check here -# return side_a_unknown -# else -# return dtu -# end -# end - """ $(TYPEDEF) $(TYPEDSIGNATURES) From dc75ab7d74929be67309b9546a94827a7ddcead3 Mon Sep 17 00:00:00 2001 From: "Craig M. Hamel" Date: Fri, 12 Jun 2026 09:25:12 -0400 Subject: [PATCH 3/4] bug in adapt bixed. --- src/DofManagers.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/DofManagers.jl b/src/DofManagers.jl index 222f66cb..a594643d 100644 --- a/src/DofManagers.jl +++ b/src/DofManagers.jl @@ -74,7 +74,6 @@ end function Adapt.adapt_structure(to, dof::DofManager) return DofManager{_is_condensed(dof)}( - adapt(to, dof), adapt(to, cache.dirichlet_dofs), adapt(to, cache.dof_to_unknown), adapt(to, cache.periodic_side_a_dofs), From 26f541add0bef5f952ddc2e436cc0c6304f14949 Mon Sep 17 00:00:00 2001 From: "Craig M. Hamel" Date: Fri, 12 Jun 2026 09:27:27 -0400 Subject: [PATCH 4/4] another adapt bugfix. --- src/DofManagers.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/DofManagers.jl b/src/DofManagers.jl index a594643d..7507f929 100644 --- a/src/DofManagers.jl +++ b/src/DofManagers.jl @@ -74,12 +74,12 @@ end function Adapt.adapt_structure(to, dof::DofManager) return DofManager{_is_condensed(dof)}( - adapt(to, cache.dirichlet_dofs), - adapt(to, cache.dof_to_unknown), - adapt(to, cache.periodic_side_a_dofs), - adapt(to, cache.periodic_side_b_dofs), - adapt(to, cache.periodic_side_b_to_side_a_unknown), - adapt(to, cache.unknown_dofs), + adapt(to, dof.dirichlet_dofs), + adapt(to, dof.dof_to_unknown), + adapt(to, dof.periodic_side_a_dofs), + adapt(to, dof.periodic_side_b_dofs), + adapt(to, dof.periodic_side_b_to_side_a_unknown), + adapt(to, dof.unknown_dofs), adapt(to, dof.var) ) end