Cthonios · cmhamel · Aug 25, 2026 · Aug 25, 2026
diff --git a/Project.toml b/Project.toml
@@ -1,6 +1,6 @@
 name = "FiniteElementContainers"
 uuid = "d08262e4-672f-4e7f-a976-f2cea5767631"
-version = "0.15.2"
+version = "0.15.3"
 authors = ["Craig M. Hamel <cmhamel32@gmail.com> and contributors"]
 
 [deps]
@@ -16,6 +16,7 @@ GPUArrays = "0c68f7d7-f131-5f86-a1c3-88cf8149b2d7"
 KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
 Krylov = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Preferences = "21216c6a-2e73-6563-6e65-726566657250"
 ReferenceFiniteElements = "6dc62d09-f8eb-43fd-9672-074e490a997f"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SparseMatricesCSR = "a0a7dd2c-ebf4-11e9-1f05-cf50bc540ca1"
@@ -60,6 +61,7 @@ Krylov = "0.10"
 LinearAlgebra = "1"
 Metis = "1"
 PartitionedArrays = "0.5"
+Preferences = "1"
 ReferenceFiniteElements = "0.14"
 SparseArrays = "1"
 SparseMatricesCSR = "0.6"

diff --git a/src/FiniteElementContainers.jl b/src/FiniteElementContainers.jl
@@ -188,13 +188,41 @@ using ForwardDiff
 using GPUArrays
 using Krylov
 using LinearAlgebra
+using Preferences
 using ReferenceFiniteElements
 using SparseArrays
 using SparseMatricesCSR
 using StaticArrays
 using Tensors
 using TimerOutputs
 
+##################
+# preferences
+##################
+const ASSEMBLE_DIAGONAL_GPU_BLOCK_SIZE = @load_preference("assemble_diagonal_gpu_block_size", 256)
+const ASSEMBLE_LUMPED_MASS_GPU_BLOCK_SIZE = @load_preference("assemble_lumped_mass_gpu_block_size", 256)
+const ASSEMBLE_MATRIX_GPU_BLOCK_SIZE = @load_preference("assemble_matrix_gpu_block_size", 256)
+const ASSEMBLE_MATRIX_ACTION_GPU_BLOCK_SIZE = @load_preference("assemble_matrix_action_gpu_block_size", 256)
+const ASSEMBLE_MATRIX_FREE_ACTION_GPU_BLOCK_SIZE = @load_preference("assemble_matrix_free_action_gpu_block_size", 256)
+const ASSEMBLE_MATRIX_WEAKLY_ENFORCED_BC_GPU_BLOCK_SIZE = @load_preference("assemble_matrix_weakly_enforced_bc_gpu_block_size", 256)
+const ASSEMBLE_QUADRATURE_QUANTITY_GPU_BLOCK_SIZE = @load_preference("assemble_quadrature_quantity_gpu_block_size", 256)
+const ASSEMBLE_VECTOR_GPU_BLOCK_SIZE = @load_preference("assemble_vector_gpu_block_size", 256)
+const ASSEMBLE_VECTOR_SOURCE_GPU_BLOCK_SIZE = @load_preference("assemble_vector_source_gpu_block_size", 256)
+const ASSEMBLE_VECTOR_WEAKLY_ENFORCE_BC_GPU_BLOCK_SIZE = @load_preference("assemble_vector_weakly_enforced_bc_gpu_block_size", 256)
+
+function summarize_preferences()
+    println("GPU Preferences:")
+    println("  assemble_diagonal_gpu_block_size                  = ", ASSEMBLE_DIAGONAL_GPU_BLOCK_SIZE)
+    println("  assemble_lumped_mass_gpu_block_size               = ", ASSEMBLE_LUMPED_MASS_GPU_BLOCK_SIZE)
+    println("  assemble_matrix_gpu_block_size                    = ", ASSEMBLE_MATRIX_GPU_BLOCK_SIZE)
+    println("  assemble_matrix_action_gpu_block_size             = ", ASSEMBLE_MATRIX_ACTION_GPU_BLOCK_SIZE)
+    println("  assemble_matrix_free_action_gpu_block_size        = ", ASSEMBLE_MATRIX_FREE_ACTION_GPU_BLOCK_SIZE)
+    println("  assemble_matrix_weakly_enforced_bc_gpu_block_size = ", ASSEMBLE_MATRIX_WEAKLY_ENFORCED_BC_GPU_BLOCK_SIZE)
+    println("  assemble_quadrature_quantity_gpu_block_size       = ", ASSEMBLE_QUADRATURE_QUANTITY_GPU_BLOCK_SIZE)
+    println("  assemble_vector_gpu_block_size                    = ", ASSEMBLE_VECTOR_GPU_BLOCK_SIZE)
+    println("  assemble_vector_source_gpu_block_size             = ", ASSEMBLE_VECTOR_SOURCE_GPU_BLOCK_SIZE)
+    println("  assemble_vector_weakly_enforced_bc_gpu_block_size = ", ASSEMBLE_VECTOR_WEAKLY_ENFORCE_BC_GPU_BLOCK_SIZE)
+end
 ##################
 # exceptions
 ##################

diff --git a/src/assemblers/Assemblers.jl b/src/assemblers/Assemblers.jl
@@ -399,15 +399,16 @@ function _assemble_block!(
   X::AbstractField, t::T, dt::T,
   U::Solution, U_old::Solution, 
   state_old::StateVariableField, state_new::StateVariableField, props::PropertyField,
-  return_type::R
+  return_type::R;
+  gpu_block_size = 256
 ) where {
   T        <: Number,
   Solution <: AbstractField,
   R        <: AssembledReturnType
 }
   conns = conns_all.data
   coffset = conns_all.offsets[b]
-  foreach_element(conns_all, b) do e
+  foreach_element(conns_all, b; block_size = gpu_block_size) do e
     conn = connectivity(ref_fe, conns, e, coffset)
     x_el, u_el, u_el_old = element_level_fields(ref_fe, conn, X, U, U_old)
     props_el = properties(props, e, b)
@@ -431,14 +432,15 @@ function _assemble_block!(
   t::T, Δt::T,
   props::PropertyField, state_old::StateVariableField, state_new::StateVariableField,
   conns_all, ref_fe::ReferenceFE,
-  X::AbstractField, U::Solution, U_old::Solution
+  X::AbstractField, U::Solution, U_old::Solution;
+  gpu_block_size = 256
 ) where {
   T        <: Number,
   Solution <: AbstractField
 }
   conns = conns_all.data
   coffset = conns_all.offsets[b]
-  foreach_element(conns_all, b) do e
+  foreach_element(conns_all, b; block_size = gpu_block_size) do e
     conn = connectivity(ref_fe, conns, e, coffset)
     x_el, u_el, u_el_old = element_level_fields(ref_fe, conn, X, U, U_old)
     props_el = properties(props, e, b)

diff --git a/src/assemblers/Diagonal.jl b/src/assemblers/Diagonal.jl
@@ -47,7 +47,8 @@ function assemble_diagonal!(
         physics,
         t, Δt,
         p.properties, p.state_old, p.state_new,
-        conns, ref_fe, X, U, U_old
+        conns, ref_fe, X, U, U_old;
+        gpu_block_size = ASSEMBLE_DIAGONAL_GPU_BLOCK_SIZE
       )
     else
       _assemble_block!(
@@ -59,7 +60,8 @@ function assemble_diagonal!(
         X, t, Δt,
         U, U_old,
         p.state_old, p.state_new, p.properties,
-        return_type
+        return_type;
+        gpu_block_size = ASSEMBLE_DIAGONAL_GPU_BLOCK_SIZE
       )
     end
   end

diff --git a/src/assemblers/LumpedMass.jl b/src/assemblers/LumpedMass.jl
@@ -51,7 +51,8 @@ function assemble_lumped_mass!(
       X, t, Δt,
       U, U_old,
       p.state_old, p.state_new, p.properties,
-      return_type
+      return_type;
+      gpu_block_size = ASSEMBLE_LUMPED_MASS_GPU_BLOCK_SIZE
     )
   end
   return nothing

diff --git a/src/assemblers/Matrix.jl b/src/assemblers/Matrix.jl
@@ -56,7 +56,8 @@ function assemble_matrix!(
         t, dt,
         p.properties, p.state_old, p.state_new,
         conns,
-        ref_fe, X, U, U_old
+        ref_fe, X, U, U_old;
+        gpu_block_size = ASSEMBLE_MATRIX_GPU_BLOCK_SIZE
       )
     else
       _assemble_block!(
@@ -68,7 +69,8 @@ function assemble_matrix!(
         X, t, dt,
         U, U_old, 
         p.state_old, p.state_new, p.properties,
-        return_type
+        return_type;
+        gpu_block_size = ASSEMBLE_MATRIX_GPU_BLOCK_SIZE
       )
     end
   end

diff --git a/src/assemblers/MatrixAction.jl b/src/assemblers/MatrixAction.jl
@@ -42,7 +42,8 @@ function assemble_matrix_free_action!(
       physics, ref_fe,
       X, t, Δt,
       U, U_old, V,
-      p.state_old, p.state_new, p.properties
+      p.state_old, p.state_new, p.properties;
+      gpu_block_size = ASSEMBLE_MATRIX_FREE_ACTION_GPU_BLOCK_SIZE
     )
   end
 end
@@ -55,14 +56,15 @@ function _assemble_block_matrix_free_action!(
   physics::AbstractPhysics, ref_fe::ReferenceFE,
   X::AbstractField, t::T, Δt::T,
   U::Solution, U_old::Solution, V::Solution,
-  state_old::StateVariableField, state_new::StateVariableField, props::AbstractArray
+  state_old::StateVariableField, state_new::StateVariableField, props::AbstractArray;
+  gpu_block_size = 256
 ) where {
   T        <: Number,
   Solution <: AbstractField
 }
   conns = conns_all.data
   coffset = conns_all.offsets[b]
-  foreach_element(conns_all, b) do e
+  foreach_element(conns_all, b; block_size = gpu_block_size) do e
     conn = connectivity(ref_fe, conns, e, coffset)
     x_el, u_el, u_el_old, v_el = element_level_fields(ref_fe, conn, X, U, U_old, V)
     props_el = properties(props, e, b)
@@ -134,7 +136,8 @@ function assemble_matrix_free_action_full!(
       physics, ref_fe,
       X, t, Δt,
       U, U_old, V,
-      p.state_old, p.state_new, p.properties
+      p.state_old, p.state_new, p.properties;
+      gpu_block_size = ASSEMBLE_MATRIX_FREE_ACTION_GPU_BLOCK_SIZE
     )
   end
   # The free-DOF entry point above relies on `p.hvp_scratch_field`'s BC
@@ -193,7 +196,8 @@ function assemble_matrix_action!(
         t, Δt,
         p.properties, p.state_old, p.state_new,
         conns,
-        ref_fe, X, U, U_old, V
+        ref_fe, X, U, U_old, V;
+        gpu_block_size = ASSEMBLE_MATRIX_ACTION_GPU_BLOCK_SIZE
       )
     else
       _assemble_block_matrix_action!(
@@ -204,7 +208,8 @@ function assemble_matrix_action!(
         physics, ref_fe,
         X, t, Δt,
         U, U_old, V,
-        p.state_old, p.state_new, p.properties
+        p.state_old, p.state_new, p.properties;
+        gpu_block_size = ASSEMBLE_MATRIX_ACTION_GPU_BLOCK_SIZE
       )
     end
   end
@@ -218,14 +223,15 @@ function _assemble_block_matrix_action!(
   physics::AbstractPhysics, ref_fe::ReferenceFE,
   X::AbstractField, t::T, Δt::T,
   U::Solution, U_old::Solution, V::Solution,
-  state_old::StateVariableField, state_new::StateVariableField, props::AbstractArray
+  state_old::StateVariableField, state_new::StateVariableField, props::AbstractArray;
+  gpu_block_size = 256
 ) where {
   T        <: Number,
   Solution <: AbstractField
 }
   conns = conns_all.data
   coffset = conns_all.offsets[b]
-  foreach_element(conns_all, b) do e
+  foreach_element(conns_all, b; block_size = gpu_block_size) do e
     conn = connectivity(ref_fe, conns, e, coffset)
     x_el, u_el, u_el_old, v_el = element_level_fields(ref_fe, conn, X, U, U_old, V)
     props_el = properties(props, e, b)
@@ -251,14 +257,15 @@ function _assemble_block_matrix_action!(
   t::T, Δt::T,
   props::PropertyField, state_old::StateVariableField, state_new::StateVariableField,
   conns_all, ref_fe::ReferenceFE,
-  X::AbstractField, U::Solution, U_old::Solution, V::Solution
+  X::AbstractField, U::Solution, U_old::Solution, V::Solution;
+  gpu_block_size = 256
 ) where {
   T        <: Number,
   Solution <: AbstractField
 }
   conns = conns_all.data
   coffset = conns_all.offsets[b]
-  foreach_element(conns_all, b) do e
+  foreach_element(conns_all, b; block_size = gpu_block_size) do e
     conn = connectivity(ref_fe, conns, e, coffset)
     x_el, u_el, u_el_old, v_el = element_level_fields(ref_fe, conn, X, U, U_old, V)
     props_el = properties(props, e, b)

diff --git a/src/assemblers/QuadratureQuantity.jl b/src/assemblers/QuadratureQuantity.jl
@@ -39,13 +39,15 @@ function assemble_quadrature_quantity!(
       X, t, Δt,
       U, U_old,
       p.state_old, p.state_new, p.properties,
-      return_type
+      return_type;
+      gpu_block_size = ASSEMBLE_QUADRATURE_QUANTITY_GPU_BLOCK_SIZE
     )
   end
 end
 
 """
 $(TYPEDSIGNATURES)
+Consider deprecating
 """
 function assemble_quadrature_quantity!(
   storage::NamedTuple, pattern, dof,
@@ -76,7 +78,8 @@ function assemble_quadrature_quantity!(
       X, t, Δt,
       U, U_old,
       p.state_old, p.state_new, p.properties,
-      return_type
+      return_type;
+      gpu_block_size = ASSEMBLE_QUADRATURE_QUANTITY_GPU_BLOCK_SIZE
     )
   end
 end
diff --git a/src/assemblers/Source.jl b/src/assemblers/Source.jl
@@ -46,7 +46,7 @@ function _assemble_block_vector_source!(
   field::AbstractField, conns, coffset, ref_fe, X, U, vals
 )
   # need to maybe make this fec_foreach
-  fec_foraxes(vals, 2) do e
+  fec_foraxes(vals, 2; block_size = ASSEMBLE_VECTOR_SOURCE_GPU_BLOCK_SIZE) do e
     conn = connectivity(ref_fe, conns, e, coffset)
     X_el = _element_level_fields_flat(X, ref_fe, conn)
 

diff --git a/src/assemblers/Vector.jl b/src/assemblers/Vector.jl
@@ -54,7 +54,8 @@ function assemble_vector!(
         t, Δt,
         p.properties, p.state_old, p.state_new,
         conns,
-        ref_fe, X, U, U_old
+        ref_fe, X, U, U_old;
+        gpu_block_size = ASSEMBLE_VECTOR_GPU_BLOCK_SIZE
       )
     else
       _assemble_block!(
@@ -66,7 +67,8 @@ function assemble_vector!(
         X, t, Δt,
         U, U_old,
         p.state_old, p.state_new, p.properties,
-        return_type
+        return_type;
+        gpu_block_size = ASSEMBLE_VECTOR_GPU_BLOCK_SIZE
       )
     end
   end

diff --git a/src/assemblers/WeaklyEnforcedBCs.jl b/src/assemblers/WeaklyEnforcedBCs.jl
@@ -65,7 +65,7 @@ function _assemble_block_vector_weakly_enforced_bc!(
   U::AbstractField, X::AbstractField,
   conns, ref_fe, sides, vals
 )
-  fec_foreach(sides) do e
+  fec_foreach(sides; block_size = ASSEMBLE_VECTOR_WEAKLY_ENFORCE_BC_GPU_BLOCK_SIZE) do e
     side = sides[e]
     conn = connectivity(ref_fe, conns, e, 1) # 1 for coffset
     surf_conns = surface_connectivity(ref_fe, conns, side, e, 1) # 1 for coffset
@@ -121,7 +121,7 @@ function _assemble_block_matrix_weakly_enforced_bc!(
   ND   = size(dof, 1)
   NEPE = ReferenceFiniteElements.num_cell_dofs(ref_fe)
 
-  fec_foreach(sides) do e
+  fec_foreach(sides; block_size = ASSEMBLE_MATRIX_WEAKLY_ENFORCED_BC_GPU_BLOCK_SIZE) do e
     side  = sides[e]
     el_id = elements[e]   # local-to-block index — matches SparseMatrixPattern's element loop
     conn  = connectivity(ref_fe, conns, e, 1)

diff --git a/test/TestPreferences.jl b/test/TestPreferences.jl
@@ -0,0 +1,3 @@
+@testitem "test_summarize_prefs" begin
+    FiniteElementContainers.summarize_preferences()
+end