Merge pull request #180 from Cthonios/multi-block-fixes

fixes to multi-block stiffness assembly methods and adding tests for …

Author: cmhamel

ext/FiniteElementContainersAdaptExt.jl

@@ -28,8 +28,8 @@ function Adapt.adapt_structure(to, asm::FiniteElementContainers.SparsityPattern)
   Is = adapt(to, asm.Is)
   Js = adapt(to, asm.Js)
   unknown_dofs = adapt(to, asm.unknown_dofs)
-  block_sizes = adapt(to, asm.block_sizes)
-  block_offsets = adapt(to, asm.block_offsets)
+  block_start_indices = adapt(to, asm.block_start_indices)
+  block_el_level_sizes = adapt(to, asm.block_el_level_sizes)
   #
   klasttouch = adapt(to, asm.klasttouch)
   csrrowptr = adapt(to, asm.csrrowptr)
@@ -41,7 +41,7 @@ function Adapt.adapt_structure(to, asm::FiniteElementContainers.SparsityPattern)
   cscnzval = adapt(to, asm.cscnzval)
   return FiniteElementContainers.SparsityPattern(
     Is, Js,
-    unknown_dofs, block_sizes, block_offsets,
+    unknown_dofs, block_start_indices, block_el_level_sizes,
     klasttouch, csrrowptr, csrcolval, csrnzval,
     csccolptr, cscrowval, cscnzval
   )

src/assemblers/Matrix.jl

@@ -137,11 +137,11 @@ KA.@kernel function _assemble_block_matrix_kernel!(
     end
   end
 
-  block_size = values(pattern.block_sizes)[block_id]
-  block_offset = values(pattern.block_offsets)[block_id]
-  start_id = (block_id - 1) * block_size + 
-             (E - 1) * block_offset + 1
-  end_id = start_id + block_offset - 1
+  block_start_index = values(pattern.block_start_indices)[block_id]
+  block_el_level_size = values(pattern.block_el_level_sizes)[block_id]
+  start_id = block_start_index + 
+             (E - 1) * block_el_level_size
+  end_id = start_id + block_el_level_size - 1
   ids = start_id:end_id
   for (i, id) in enumerate(ids)
     Atomix.@atomic field[id] += K_el.data[i]

src/assemblers/SparseMatrixAssembler.jl

@@ -97,12 +97,11 @@ function _assemble_element!(
   pattern::SparsityPattern, storage, K_el::SMatrix, el_id::Int, block_id::Int
 )
   # figure out ids needed to update
-  block_size = values(pattern.block_sizes)[block_id]
-  block_offset = values(pattern.block_offsets)[block_id]
-  # get range of ids
-  start_id = (block_id - 1) * block_size + 
-             (el_id - 1) * block_offset + 1
-  end_id = start_id + block_offset - 1
+  block_start_index = values(pattern.block_start_indices)[block_id]
+  block_el_level_size = values(pattern.block_el_level_sizes)[block_id]
+  start_id = block_start_index + 
+             (el_id - 1) * block_el_level_size
+  end_id = start_id + block_el_level_size - 1
   ids = start_id:end_id
 
   # get appropriate storage and update values

src/assemblers/SparsityPattern.jl

@@ -13,8 +13,8 @@ struct SparsityPattern{
   Is::I
   Js::I
   unknown_dofs::I
-  block_sizes::B
-  block_offsets::B
+  block_start_indices::B
+  block_el_level_sizes::B
   # cache arrays
   klasttouch::I
   csrrowptr::I
@@ -46,23 +46,33 @@ function SparsityPattern(dof::DofManager)
 
   # first get total number of entries in a stupid manner
   n_entries = 0
-  block_sizes = Vector{Int64}(undef, n_blocks)
-  block_offsets = Vector{Int64}(undef, n_blocks)
+  block_start_indices = Vector{Int64}(undef, n_blocks)
+  block_el_level_sizes = Vector{Int64}(undef, n_blocks)
 
   # for (n, conn) in enumerate(values(dof.H1_vars[1].fspace.elem_conns))
+  start_carry = 1
   for (n, conn) in enumerate(values(fspace.elem_conns))
     ids = reshape(1:n_total_dofs, ND, NN)
     # TODO do we need this operation?
     conn = reshape(ids[:, conn], ND * size(conn, 1), size(conn, 2))
     n_entries += size(conn, 1)^2 * size(conn, 2)
-    block_sizes[n] = size(conn, 2)
-    block_offsets[n] = size(conn, 1)^2
+
+    # block_start_indices[n] = size(conn, 2)
+    # if n == 1
+    #   block_start_indices[n] = 1
+    # else
+    #   block_start_indices[n] = carry
+    # end
+    block_start_indices[n] = start_carry
+  
+    start_carry = start_carry + size(conn, 1)^2 * size(conn, 2)
+    block_el_level_sizes[n] = size(conn, 1)^2
   end
 
   # convert to NamedTuples so it's easy to index
   block_syms = keys(fspace.ref_fes)
-  block_sizes = NamedTuple{block_syms}(tuple(block_sizes)...)
-  block_offsets = NamedTuple{block_syms}(tuple(block_offsets)...)
+  block_start_indices = NamedTuple{block_syms}(tuple(block_start_indices)...)
+  block_el_level_sizes = NamedTuple{block_syms}(tuple(block_el_level_sizes)...)
 
   # setup pre-allocated arrays based on number of entries found above
   Is = Vector{Int64}(undef, n_entries)
@@ -100,7 +110,7 @@ function SparsityPattern(dof::DofManager)
   return SparsityPattern(
     Is, Js, 
     unknown_dofs, 
-    block_sizes, block_offsets, 
+    block_start_indices, block_el_level_sizes, 
     # cache arrays
     klasttouch, csrrowptr, csrcolval, csrnzval,
     # additional cache arrays

test/poisson/TestPoisson.jl

@@ -1,3 +1,5 @@
+using Adapt
+using AMDGPU
 using Exodus
 using FiniteElementContainers
 using StaticArrays
@@ -126,3 +128,97 @@ function test_poisson_neumann(
   # rm(output_file; force=true)
   display(solver.timer)
 end
+
+function test_poisson_dirichlet_multi_block_quad4_quad4(
+  dev, use_condensed,
+  nsolver, lsolver
+)
+  mesh = UnstructuredMesh(Base.source_dir() * "/poisson/multi_block_mesh_quad4_quad4.g")
+  V = FunctionSpace(mesh, H1Field, Lagrange) 
+  physics = Poisson()
+  props = create_properties(physics)
+  u = ScalarFunction(V, :u)
+  asm = SparseMatrixAssembler(u; use_condensed=use_condensed)
+
+  # setup and update bcs
+  dbcs = DirichletBC[
+    DirichletBC(:u, :boundary, bc_func)
+  ]
+
+  # 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 = 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",), 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_dirichlet_multi_block_quad4_tri3(
+  dev, use_condensed,
+  nsolver, lsolver
+)
+  mesh = UnstructuredMesh(Base.source_dir() * "/poisson/multi_block_mesh_quad4_tri3.g")
+  V = FunctionSpace(mesh, H1Field, Lagrange) 
+  physics = Poisson()
+  props = create_properties(physics)
+  u = ScalarFunction(V, :u)
+  asm = SparseMatrixAssembler(u; use_condensed=use_condensed)
+
+  # setup and update bcs
+  dbcs = DirichletBC[
+    DirichletBC(:u, :boundary, bc_func)
+  ]
+
+  # 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 = 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",), U)
+  close(pp)
+
+  # if !Sys.iswindows()
+  #   @test exodiff(output_file, gold_file)
+  # end
+  rm(output_file; force=true)
+  display(solver.timer)
+end

test/poisson/multi_block_mesh_quad4_quad4.g

@@ -0,0 +1,38 @@
+reset
+
+create surface rectangle width 1 height 1 zplane
+Surface 1 Move X 0.5 Y 0.5
+
+create surface circle radius 0.25 zplane
+Surface 2 Move X 0.5 Y 0.5
+
+subtract body 2 from body 1 keep_tool
+
+imprint body all
+merge tol 1e-6
+merge all
+
+surface 2 scheme trimesh
+surface 2 size 0.05
+surface 3 size 0.05
+
+mesh surface 3
+mesh surface 2
+
+block 1 add surface 3
+block 2 add surface 2
+
+block 1 name "block_1"
+block 2 name "block_2"
+
+block 1 element type quad4
+block 2 element type tri3
+#block 2 element type quad4
+
+nodeset 1 add curve 1 2 3 4
+nodeset 1 name "boundary"
+sideset 1 add curve 1 2 3 4
+sideset 1 name "boundary"
+
+#export genesis "multi_block_mesh_quad4_quad4.g" overwrite
+export genesis "multi_block_mesh_quad4_tri3.g" overwrite

test/poisson/multi_block_mesh_quad4_tri3.g

@@ -35,30 +35,37 @@ include("TestPhysics.jl")
   include("poisson/TestPoisson.jl")
 
   cg_solver = x -> IterativeLinearSolver(x, :CgSolver)
-
+  condensed = [false, true]
+  lsolvers = [cg_solver, DirectLinearSolver]
   # cpu tests
-  test_poisson_dirichlet(cpu, false, NewtonSolver, DirectLinearSolver)
-  test_poisson_dirichlet(cpu, true, NewtonSolver, DirectLinearSolver)
-  test_poisson_dirichlet(cpu, false, NewtonSolver, cg_solver)
-  test_poisson_dirichlet(cpu, true, NewtonSolver, cg_solver)
-  test_poisson_neumann(cpu, false, NewtonSolver, DirectLinearSolver)
-  test_poisson_neumann(cpu, true, NewtonSolver, DirectLinearSolver)
-  test_poisson_neumann(cpu, false, NewtonSolver, cg_solver)
-  test_poisson_neumann(cpu, true, NewtonSolver, cg_solver)
-  
-  if AMDGPU.functional()
-    test_poisson_dirichlet(rocm, false, NewtonSolver, cg_solver)
-    test_poisson_dirichlet(rocm, true, NewtonSolver, cg_solver)
-    test_poisson_neumann(rocm, false, NewtonSolver, cg_solver)
-    test_poisson_neumann(rocm, true, NewtonSolver, cg_solver)
+  for cond in condensed
+    for lsolver in lsolvers
+      test_poisson_dirichlet(cpu, cond, NewtonSolver, lsolver)
+      test_poisson_dirichlet_multi_block_quad4_quad4(cpu, cond, NewtonSolver, lsolver)
+      test_poisson_dirichlet_multi_block_quad4_tri3(cpu, cond, NewtonSolver, lsolver)
+      test_poisson_neumann(cpu, cond, NewtonSolver, lsolver)
+    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)
+  if AMDGPU.functional()
+    for cond in condensed
+      test_poisson_dirichlet(rocm, cond, NewtonSolver, cg_solver)
+      test_poisson_dirichlet_multi_block_quad4_quad4(rocm, cond, NewtonSolver, cg_solver)
+      test_poisson_dirichlet_multi_block_quad4_tri3(rocm, cond, NewtonSolver, cg_solver)
+      test_poisson_neumann(rocm, cond, NewtonSolver, cg_solver)
+    end
+  #   test_poisson_dirichlet(rocm, false, NewtonSolver, cg_solver)
+  #   test_poisson_dirichlet(rocm, true, NewtonSolver, cg_solver)
+  #   test_poisson_neumann(rocm, false, NewtonSolver, cg_solver)
+  #   test_poisson_neumann(rocm, true, NewtonSolver, cg_solver)
   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
 end
 
 @testset ExtendedTestSet "Mechanics Problem" begin