Merge pull request #177 from Cthonios/sparse-more-testing

a bunch of re-work for poisson regression tests.

Author: cmhamel

src/DofManagers.jl

@@ -212,3 +212,9 @@ function update_field_unknowns!(
     _update_field_unknowns!(U, dof, Uu, backend)
     return nothing
 end
+
+function update_field_unknowns!(
+    U::F, dof::DofManager, Uu::F
+) where F <: AbstractField
+    update_field_unknowns!(U, dof, Uu.data)
+end

src/assemblers/Assemblers.jl

@@ -215,6 +215,13 @@ end
 #   end
 # end
 
+"""
+$(TYPEDSIGNATURES)
+"""
+function hessian(asm::AbstractAssembler)
+  return _hessian(asm, KA.get_backend(asm))
+end
+
 # new approach requiring access to the v that makes Hv
 """
 $(TYPEDSIGNATURES)

src/assemblers/SparseMatrixAssembler.jl

@@ -187,11 +187,21 @@ function _adjust_matrix_entries_for_constraints(
   return nothing
 end
 
+function _hessian(assembler::SparseMatrixAssembler, ::KA.CPU)
+  H = SparseArrays.sparse!(assembler.pattern, assembler.hessian_storage)
+
+  if _is_condensed(assembler.dof)
+    _adjust_matrix_entries_for_constraints!(H, assembler.constraint_storage, KA.get_backend(assembler))
+  end
+
+  return H
+end
+
 function _mass(assembler::SparseMatrixAssembler, ::KA.CPU)
   M = SparseArrays.sparse!(assembler.pattern, assembler.mass_storage)
 
   if _is_condensed(assembler.dof)
-    _adjust_matrix_entries_for_constraints!(M, assembler.constraint_storage, KA.get_getbackend(assembler))
+    _adjust_matrix_entries_for_constraints!(M, assembler.constraint_storage, KA.get_backend(assembler))
   end
 
   return M

src/bcs/BoundaryConditions.jl

@@ -91,6 +91,33 @@ function BCBookKeeping(
   )
 end
 
+function BCBookKeeping(mesh, sset_name::Symbol)
+    # 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)))
+      end
+    end
+
+    dofs = Vector{Int64}(undef, 0)
+      return BCBookKeeping(
+      blocks, dofs, elements, nodes, sides, side_nodes
+    )
+end
+
 """
 $(TYPEDSIGNATURES)
 """

src/bcs/DirichletBCs.jl

@@ -140,6 +140,14 @@ function _update_field_dirichlet_bcs!(U, bc::DirichletBCContainer, ::KA.CPU)
   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)
+    U[dof] = val
+    V[dof] = val_dot
+  end
+  return nothing
+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)
     U[dof] = val
@@ -188,6 +196,12 @@ function update_field_dirichlet_bcs!(U, bcs::NamedTuple)
   return nothing
 end
 
+function update_field_dirichlet_bcs!(U, V, bcs::NamedTuple)
+  for bc in values(bcs)
+    _update_field_dirichlet_bcs!(U, V, bc, KA.get_backend(bc))
+  end
+end
+
 function update_field_dirichlet_bcs!(U, V, A, bcs::NamedTuple)
   for bc in values(bcs)
     _update_field_dirichlet_bcs!(U, V, A, bc, KA.get_backend(bc))

src/fields/L2QuadratureField.jl

@@ -37,6 +37,11 @@ end
 #   setindex!(field.vals, v, (e - 1) * NQ + (q - 1) * NF + n)
 # end
 
+function Base.setindex!(field::L2QuadratureField{T, D, NF, NQ}, v, n::Int, q::Int, e::Int) where {T, D, NF, NQ}
+  setindex!(field.data, v, (e - 1) * NQ + (q - 1) * NF + n)
+  return nothing
+end
+
 function Base.size(field::L2QuadratureField{T, D, NF, NQ}) where {T, D, NF, NQ} 
   if NF == 0
     (NF, NQ, length(field.data) ÷ NQ)

test/poisson/TestPoisson.jl

@@ -18,7 +18,10 @@ include("TestPoissonCommon.jl")
 
 # read mesh and relevant quantities
 
-function test_poisson_direct(use_condensed)
+function test_poisson_dirichlet(
+  dev, use_condensed,
+  nsolver, lsolver
+)
   mesh = UnstructuredMesh(mesh_file)
   V = FunctionSpace(mesh, H1Field, Lagrange) 
   physics = Poisson()
@@ -34,29 +37,41 @@ function test_poisson_direct(use_condensed)
     DirichletBC(:u, :sset_4, bc_func),
   ]
 
-  # direct solver test
   # setup the parameters
-  @show p = create_parameters(mesh, asm, physics, props; dirichlet_bcs=dbcs)
+  p = create_parameters(mesh, asm, physics, props; dirichlet_bcs=dbcs)
+
+  if dev != cpu
+    p = p |> dev
+    asm = asm |> dev 
+  end
 
   # setup solver and integrator
-  solver = NewtonSolver(DirectLinearSolver(asm))
+  solver = nsolver(lsolver(asm))
   integrator = QuasiStaticIntegrator(solver)
   evolve!(integrator, p)
 
+  if dev != cpu
+    p = p |> cpu
+  end
+
+  U = p.h1_field
+
   pp = PostProcessor(mesh, output_file, u)
   write_times(pp, 1, 0.0)
-  write_field(pp, 1, ("u",), p.h1_field)
+  write_field(pp, 1, ("u",), U)
   close(pp)
 
   if !Sys.iswindows()
     @test exodiff(output_file, gold_file)
   end
   rm(output_file; force=true)
   display(solver.timer)
-
 end
 
-function test_poisson_direct_neumman(use_condensed)
+function test_poisson_neumann(
+  dev, use_condensed,
+  nsolver, lsolver
+)
   mesh = UnstructuredMesh(mesh_file)
   V = FunctionSpace(mesh, H1Field, Lagrange) 
   physics = Poisson()
@@ -77,123 +92,86 @@ function test_poisson_direct_neumman(use_condensed)
 
   # direct solver test
   # setup the parameters
-  @show p = create_parameters(
+  p = create_parameters(
     mesh, asm, physics, props; 
     dirichlet_bcs=dbcs,
     neumann_bcs=nbcs
   )
 
+  if dev != cpu
+    p = p |> dev
+    asm = asm |> dev 
+  end
+
   # setup solver and integrator
-  solver = NewtonSolver(DirectLinearSolver(asm))
+  solver = nsolver(lsolver(asm))
   integrator = QuasiStaticIntegrator(solver)
   evolve!(integrator, p)
 
-  pp = PostProcessor(mesh, output_file, u)
-  write_times(pp, 1, 0.0)
-  write_field(pp, 1, ("u",), p.h1_field)
-  close(pp)
-
-  if !Sys.iswindows()
-    @test exodiff(output_file, gold_file)
+  if dev != cpu
+    p = p |> cpu
   end
-  rm(output_file; force=true)
-  display(solver.timer)
-
-end
-
-function test_poisson_iterative(use_condensed)
-  mesh = UnstructuredMesh(mesh_file)
-  V = FunctionSpace(mesh, H1Field, Lagrange) 
-  physics = Poisson()
-  props = create_properties(physics)
-  u = ScalarFunction(V, :u)
-  asm = SparseMatrixAssembler(u; use_condensed)
-
-  # setup and update bcs
-  dbcs = DirichletBC[
-    DirichletBC(:u, :sset_1, bc_func),
-    DirichletBC(:u, :sset_2, bc_func),
-    DirichletBC(:u, :sset_3, bc_func),
-    DirichletBC(:u, :sset_4, bc_func),
-  ]
-
-  # iterative solver test
-  # setup the parameters
-  p = create_parameters(mesh, asm, physics, props; dirichlet_bcs=dbcs)
-
-  # setup solver and integrator
-  solver = NewtonSolver(IterativeLinearSolver(asm, :CgSolver))
-  integrator = QuasiStaticIntegrator(solver)
-  @time evolve!(integrator, p)
 
-  display(solver.timer)
+  # TODO make a neumann gold file
+  # U = p.h1_field
 
-  pp = PostProcessor(mesh, output_file, u)
-  write_times(pp, 1, 0.0)
-  write_field(pp, 1, ("u",), p.h1_field)
-  close(pp)
+  # pp = PostProcessor(mesh, output_file, u)
+  # write_times(pp, 1, 0.0)
+  # write_field(pp, 1, ("u",), U)
+  # close(pp)
 
-  if !Sys.iswindows()
-    @test exodiff(output_file, gold_file)
-  end
-  rm(output_file; force=true)
+  # if !Sys.iswindows()
+  #   @test exodiff(output_file, gold_file)
+  # end
+  # rm(output_file; force=true)
   display(solver.timer)
 end
 
-@time test_poisson_direct(false)
-@time test_poisson_direct(false)
-@time test_poisson_direct(true)
-@time test_poisson_direct(true)
-# @time test_poisson_direct_neumman()
-# @time test_poisson_direct_neumman()
-@time test_poisson_iterative(false)
-@time test_poisson_iterative(false)
-@time test_poisson_iterative(true)
-@time test_poisson_iterative(true)
-
-# # condensed test
-# mesh = UnstructuredMesh(mesh_file)
-# V = FunctionSpace(mesh, H1Field, Lagrange) 
-# physics = Poisson()
-# u = ScalarFunction(V, :u)
-# asm = SparseMatrixAssembler(H1Field, u)
-# # pp = PostProcessor(mesh, output_file, u)
-
-# # setup and update bcs
-# dbcs = DirichletBC[
-#   DirichletBC(asm.dof, :u, :sset_1, bc_func),
-#   DirichletBC(asm.dof, :u, :sset_2, bc_func),
-#   DirichletBC(asm.dof, :u, :sset_3, bc_func),
-#   DirichletBC(asm.dof, :u, :sset_4, bc_func),
-# ]
-# update_dofs!(asm, dbcs; use_condensed=true)
-# Uu = create_unknowns(asm)
-# Ubc = create_bcs(asm, H1Field)
-# U = create_field(asm, H1Field)
-# update_field!(U, asm, Uu, Ubc)
-# update_field_bcs!(U, asm.dof, dbcs, 0.)
-# assemble!(asm, physics, U, :residual_and_stiffness)
-# K = stiffness(asm)
-# G = constraint_matrix(asm)
-# # @time H = (G + I) * K
-# K[asm.dof.H1_bc_dofs, asm.dof.H1_bc_dofs] .= 1.
-# # R = G * residual(asm)
-# # R = G * asm.residual_storage.vals
-# R = asm.residual_storage
-# R[asm.dof.H1_bc_dofs] .= 0.
-# ΔUu = -K \ R.vals
-# U.vals .= U.vals .+ ΔUu
-# assemble!(asm, physics, U, :residual_and_stiffness)
-# K = stiffness(asm)
-# G = constraint_matrix(asm)
-# # @time H = (G + I) * K
-# K[asm.dof.H1_bc_dofs, asm.dof.H1_bc_dofs] .= 1.
-# # R = G * residual(asm)
-# # R = G * asm.residual_storage.vals
-# R = asm.residual_storage
-# R[asm.dof.H1_bc_dofs] .= 0.
-# ΔUu = -K \ R.vals
-# U.vals .= U.vals .+ ΔUu
-# U
-# # @time H = G * K + G * I
-# # @time H = G * K
+# function test_poisson_iterative(dev, use_condensed)
+#   mesh = UnstructuredMesh(mesh_file)
+#   V = FunctionSpace(mesh, H1Field, Lagrange) 
+#   physics = Poisson()
+#   props = create_properties(physics)
+#   u = ScalarFunction(V, :u)
+#   asm = SparseMatrixAssembler(u; use_condensed)
+
+#   # setup and update bcs
+#   dbcs = DirichletBC[
+#     DirichletBC(:u, :sset_1, bc_func),
+#     DirichletBC(:u, :sset_2, bc_func),
+#     DirichletBC(:u, :sset_3, bc_func),
+#     DirichletBC(:u, :sset_4, bc_func),
+#   ]
+
+#   # iterative solver test
+#   # setup the parameters
+#   p = create_parameters(mesh, asm, physics, props; dirichlet_bcs=dbcs)
+
+#   if dev != cpu
+#     p = p |> dev
+#     asm = asm |> dev 
+#   end
+
+#   # setup solver and integrator
+#   solver = NewtonSolver(IterativeLinearSolver(asm, :CgSolver))
+#   integrator = QuasiStaticIntegrator(solver)
+#   @time evolve!(integrator, p)
+
+#   display(solver.timer)
+
+#   if dev != cpu
+#     p = p |> cpu
+#   end
+
+#   U = p.h1_field
+#   pp = PostProcessor(mesh, output_file, u)
+#   write_times(pp, 1, 0.0)
+#   write_field(pp, 1, ("u",), U)
+#   close(pp)
+
+#   if !Sys.iswindows()
+#     @test exodiff(output_file, gold_file)
+#   end
+#   rm(output_file; force=true)
+#   display(solver.timer)
+# end