lots of clean up on periodic bcs. Now working with inplace and sparse…

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,27 +26,27 @@ 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")
     PeriodicBC("u", "y", zero_func, "sset_4", "sset_2")
 ]
 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)
 
 # display(Uu)
 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)

src/DofManagers.jl

@@ -52,19 +52,22 @@ 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
@@ -81,6 +84,7 @@ function Adapt.adapt_structure(to, dof::DofManager)
     )
 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 +94,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 +174,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 +186,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 +200,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

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, 

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

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

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