Merge pull request #45 from Cthonios/add-residual-to-assembler

moving residual back to assembler.
Cthonios · Apr 3, 2024 · d0879e9 · d0879e9
2 parents 3f7476c + ca215ec
commit d0879e9
Show file tree

Hide file tree

Showing 5 changed files with 130 additions and 63 deletions.
diff --git a/Project.toml b/Project.toml
@@ -1,7 +1,7 @@
 name = "FiniteElementContainers"
 uuid = "d08262e4-672f-4e7f-a976-f2cea5767631"
 authors = ["Craig M. Hamel <cmhamel32@gmail.com> and contributors"]
-version = "0.5.1"
+version = "0.5.2"
 
 [deps]
 Atomix = "a9b6321e-bd34-4604-b9c9-b65b8de01458"
@@ -28,7 +28,7 @@ Adapt = "3"
 Aqua = "0.8"
 Atomix = "0.1.0"
 DocStringExtensions = "0.9"
-Exodus = "0.12"
+Exodus = "0.13"
 IterativeSolvers = "0.9"
 JET = "0.8"
 KernelAbstractions = "0.9"

diff --git a/src/assemblers/Assemblers.jl b/src/assemblers/Assemblers.jl
@@ -122,3 +122,54 @@ end
 include("DynamicAssemblers.jl")
 include("MatrixFreeStaticAssembler.jl")
 include("StaticAssembler.jl")
+
+"""
+$(TYPEDSIGNATURES)
+method that assumes first dof
+TODO move sorting of nodes up stream
+TODO remove other scratch unknowns and unknown_dofs arrays
+"""
+function update_unknown_dofs!(
+  assembler::Union{DynamicAssembler, StaticAssembler},
+  dof,
+  fspaces, 
+  nodes_in::V
+) where V <: AbstractVector{<:Integer}
+
+  # make this an assumption of the method
+  nodes = sort(nodes_in)
+
+  n_total_dofs = num_dofs_per_node(dof) * num_nodes(dof) - length(nodes)
+
+  # TODO change to a good sizehint!
+  resize!(assembler.Is, 0)
+  resize!(assembler.Js, 0)
+  resize!(assembler.unknown_dofs, 0)
+
+  n = 1
+  for fspace in fspaces
+    for e in 1:num_elements(fspace)
+      conn = dof_connectivity(fspace, e)
+      for temp in Iterators.product(conn, conn)
+        if insorted(temp[1], nodes) || insorted(temp[2], nodes)
+          # really do nothing here
+        else
+          push!(assembler.Is, temp[1] - count(x -> x < temp[1], nodes))
+          push!(assembler.Js, temp[2] - count(x -> x < temp[2], nodes))
+          push!(assembler.unknown_dofs, n)
+        end
+        n += 1
+      end
+    end
+  end
+
+  # resize cache arrays
+  resize!(assembler.klasttouch, n_total_dofs)
+  resize!(assembler.csrrowptr, n_total_dofs + 1)
+  resize!(assembler.csrcolval, length(assembler.Is))
+  resize!(assembler.csrnzval, length(assembler.Is))
+
+  # resize!(assembler.stiffnesses, length(assembler.Is))
+
+  return nothing
+end
diff --git a/src/assemblers/DynamicAssemblers.jl b/src/assemblers/DynamicAssemblers.jl
@@ -12,6 +12,7 @@ struct DynamicAssembler{
   U       <: AbstractArray{Itype, 1},
   Sizes   <: AbstractArray{Itype, 1},
   Offsets <: AbstractArray{Itype, 1},
+  R       <: NodalField,
   K       <: AbstractArray{Rtype, 1}, # should always be a vector type thing
   M       <: AbstractArray{Rtype, 1},
   # cache types
@@ -24,6 +25,7 @@ struct DynamicAssembler{
   unknown_dofs::U
   block_sizes::Sizes
   block_offsets::Offsets
+  residuals::R
   stiffnesses::K
   masses::M
   # cache arrays
@@ -78,6 +80,7 @@ function DynamicAssembler(dof::DofManager, fspaces::Fs) where Fs
     end
   end
 
+  residuals = create_fields(dof) #|> vec
   stiffnesses = zeros(Float64, size(Is))
   masses = zeros(Float64, size(Is))
 
@@ -94,13 +97,15 @@ function DynamicAssembler(dof::DofManager, fspaces::Fs) where Fs
   return DynamicAssembler{
     Float64, Int64, 
     typeof(Is), typeof(Js), typeof(unknown_dofs), typeof(block_sizes), typeof(block_offsets),
-    typeof(stiffnesses), typeof(masses),
+    typeof(residuals), typeof(stiffnesses), typeof(masses),
     # cache arrays
     typeof(klasttouch), typeof(csrrowptr), typeof(csrcolval), typeof(csrnzval),
     # additional cache arrays
     typeof(csccolptr), typeof(cscrowval), typeof(cscnzval)
   }(
-    Is, Js, unknown_dofs, block_sizes, block_offsets, stiffnesses, masses,
+    Is, Js, unknown_dofs, block_sizes, block_offsets, 
+    # fields
+    residuals, stiffnesses, masses,
     # cache arrays
     klasttouch, csrrowptr, csrcolval, csrnzval,
     # additional cache arrays
@@ -163,7 +168,7 @@ $(TYPEDSIGNATURES)
 Simple method for assembling in serial
 """
 function assemble!(
-  R,
+  # R,
   assembler::DynamicAssembler,
   dof::DofManager,
   fspace::FunctionSpace,
@@ -175,6 +180,8 @@ function assemble!(
   N    = num_nodes_per_element(fspace)
   NxNDof = N * NDof
 
+  R = assembler.residuals
+
   for e in 1:num_elements(fspace)
     U_el = element_level_fields(fspace, U, e)
     R_el = zeros(SVector{NxNDof, Float64})
@@ -203,19 +210,20 @@ end
 $(TYPEDSIGNATURES)
 """
 function assemble!(
-  R,
+  # R,
   assembler::DynamicAssembler,
   dof::DofManager,
   fspaces, X, U,
   residual_func, tangent_func, mass_func
 )
 
   # reset in some way
+  R = assembler.residuals
   R .= zero(eltype(R))
   assembler.stiffnesses .= zero(eltype(assembler.stiffnesses))
 
   for (block_id, fspace) in enumerate(fspaces)
-    assemble!(R, assembler, dof, fspace, X, U, block_id, residual_func, tangent_func, mass_func)
+    assemble!(assembler, dof, fspace, X, U, block_id, residual_func, tangent_func, mass_func)
   end
 
 end
diff --git a/src/assemblers/StaticAssembler.jl b/src/assemblers/StaticAssembler.jl
@@ -11,6 +11,7 @@ struct StaticAssembler{
   U       <: AbstractArray{Itype, 1},
   Sizes   <: AbstractArray{Itype, 1},
   Offsets <: AbstractArray{Itype, 1},
+  R       <: NodalField,
   K       <: AbstractArray{Rtype, 1}, # should always be a vector type thing
   # cache types
   C1, C2, C3, C4,
@@ -22,6 +23,7 @@ struct StaticAssembler{
   unknown_dofs::U
   block_sizes::Sizes
   block_offsets::Offsets
+  residuals::R
   stiffnesses::K
   # cache arrays
   klasttouch::C1
@@ -83,7 +85,7 @@ function StaticAssembler(dof::DofManager, fspaces::Fs) where Fs
     end
   end
 
-  # residuals = create_fields(dof)
+  residuals = create_fields(dof) #|> vec
   stiffnesses = zeros(Float64, size(Is))
 
   # create caches
@@ -99,13 +101,15 @@ function StaticAssembler(dof::DofManager, fspaces::Fs) where Fs
   return StaticAssembler{
     Float64, Int64, 
     typeof(Is), typeof(Js), typeof(unknown_dofs), typeof(block_sizes), typeof(block_offsets),
-    typeof(stiffnesses),
+    typeof(residuals), typeof(stiffnesses),
     # cache arrays
     typeof(klasttouch), typeof(csrrowptr), typeof(csrcolval), typeof(csrnzval),
     # additional cache arrays
     typeof(csccolptr), typeof(cscrowval), typeof(cscnzval)
   }(
-    Is, Js, unknown_dofs, block_sizes, block_offsets, stiffnesses,
+    Is, Js, unknown_dofs, block_sizes, block_offsets, 
+    # fields
+    residuals, stiffnesses,
     # cache arrays
     klasttouch, csrrowptr, csrcolval, csrnzval,
     # additional cache arrays
@@ -193,6 +197,8 @@ function assemble!(
   N    = num_nodes_per_element(fspace)
   NxNDof = N * NDof
 
+  R = assembler.residuals
+
   for e in 1:num_elements(fspace)
     U_el = element_level_fields(fspace, U, e)
     R_el = zeros(SVector{NxNDof, Float64})
@@ -221,7 +227,7 @@ Top level method using methods
 """
 # TODO figure this one out
 function assemble!(
-  R, 
+  # R, 
   assembler::StaticAssembler,
   dof::DofManager,
   fspaces, X, U,
@@ -230,6 +236,7 @@ function assemble!(
 
   # reset in some way
   # assembler.residuals .= 0.
+  R = assembler.residuals
   R .= zero(eltype(R))
   assembler.stiffnesses .= zero(eltype(assembler.stiffnesses))
 
@@ -239,53 +246,53 @@ function assemble!(
 
 end
 
-"""
-$(TYPEDSIGNATURES)
-method that assumes first dof
-TODO move sorting of nodes up stream
-TODO remove other scratch unknowns and unknown_dofs arrays
-"""
-function update_unknown_dofs!(
-  assembler::StaticAssembler,
-  dof,
-  fspaces, 
-  nodes_in::V
-) where V <: AbstractVector{<:Integer}
-
-  # make this an assumption of the method
-  nodes = sort(nodes_in)
-
-  n_total_dofs = num_dofs_per_node(dof) * num_nodes(dof) - length(nodes)
-
-  # TODO change to a good sizehint!
-  resize!(assembler.Is, 0)
-  resize!(assembler.Js, 0)
-  resize!(assembler.unknown_dofs, 0)
-
-  n = 1
-  for fspace in fspaces
-    for e in 1:num_elements(fspace)
-      conn = dof_connectivity(fspace, e)
-      for temp in Iterators.product(conn, conn)
-        if insorted(temp[1], nodes) || insorted(temp[2], nodes)
-          # really do nothing here
-        else
-          push!(assembler.Is, temp[1] - count(x -> x < temp[1], nodes))
-          push!(assembler.Js, temp[2] - count(x -> x < temp[2], nodes))
-          push!(assembler.unknown_dofs, n)
-        end
-        n += 1
-      end
-    end
-  end
-
-  # resize cache arrays
-  resize!(assembler.klasttouch, n_total_dofs)
-  resize!(assembler.csrrowptr, n_total_dofs + 1)
-  resize!(assembler.csrcolval, length(assembler.Is))
-  resize!(assembler.csrnzval, length(assembler.Is))
-
-  # resize!(assembler.stiffnesses, length(assembler.Is))
-
-  return nothing
-end
+# """
+# $(TYPEDSIGNATURES)
+# method that assumes first dof
+# TODO move sorting of nodes up stream
+# TODO remove other scratch unknowns and unknown_dofs arrays
+# """
+# function update_unknown_dofs!(
+#   assembler::StaticAssembler,
+#   dof,
+#   fspaces, 
+#   nodes_in::V
+# ) where V <: AbstractVector{<:Integer}
+
+#   # make this an assumption of the method
+#   nodes = sort(nodes_in)
+
+#   n_total_dofs = num_dofs_per_node(dof) * num_nodes(dof) - length(nodes)
+
+#   # TODO change to a good sizehint!
+#   resize!(assembler.Is, 0)
+#   resize!(assembler.Js, 0)
+#   resize!(assembler.unknown_dofs, 0)
+
+#   n = 1
+#   for fspace in fspaces
+#     for e in 1:num_elements(fspace)
+#       conn = dof_connectivity(fspace, e)
+#       for temp in Iterators.product(conn, conn)
+#         if insorted(temp[1], nodes) || insorted(temp[2], nodes)
+#           # really do nothing here
+#         else
+#           push!(assembler.Is, temp[1] - count(x -> x < temp[1], nodes))
+#           push!(assembler.Js, temp[2] - count(x -> x < temp[2], nodes))
+#           push!(assembler.unknown_dofs, n)
+#         end
+#         n += 1
+#       end
+#     end
+#   end
+
+#   # resize cache arrays
+#   resize!(assembler.klasttouch, n_total_dofs)
+#   resize!(assembler.csrrowptr, n_total_dofs + 1)
+#   resize!(assembler.csrcolval, length(assembler.Is))
+#   resize!(assembler.csrnzval, length(assembler.Is))
+
+#   # resize!(assembler.stiffnesses, length(assembler.Is))
+
+#   return nothing
+# end
diff --git a/test/poisson/TestPoisson.jl b/test/poisson/TestPoisson.jl
@@ -58,9 +58,10 @@ for type in types
   function solve(asm, dof, fspaces, X, U, Uu)
     for n in 1:10
       update_fields!(U, dof, Uu)
-      assemble!(R, asm, dof, fspaces, X, U, residual, tangent)
+      # assemble!(R, asm, dof, fspaces, X, U, residual, tangent)
+      assemble!(asm, dof, fspaces, X, U, residual, tangent)
       # R = asm.residuals[dof.unknown_dofs]
-      R_view = @views R[dof.unknown_dofs]
+      R_view = @views asm.residuals[dof.unknown_dofs]
       K = sparse(asm)
       cg!(ΔUu, -K, R_view)
       @show norm(ΔUu) norm(R_view)