diff --git a/src/PotentialCalculation/PotentialCalculationSetup/BoundaryConditions/BoundaryConditions.jl b/src/PotentialCalculation/PotentialCalculationSetup/BoundaryConditions/BoundaryConditions.jl
index 46b681ce3..166de77d9 100644
--- a/src/PotentialCalculation/PotentialCalculationSetup/BoundaryConditions/BoundaryConditions.jl
+++ b/src/PotentialCalculation/PotentialCalculationSetup/BoundaryConditions/BoundaryConditions.jl
@@ -1,2 +1,3 @@
-include("BoundaryConditionsCylindrical.jl")
-include("BoundaryConditionsCartesian.jl")
+include("BoundaryConditionsCylindrical.jl")
+include("BoundaryConditionsCartesian.jl")
+include("BoundaryConditionsGPU.jl")
diff --git a/src/PotentialCalculation/PotentialCalculationSetup/BoundaryConditions/BoundaryConditionsCylindrical.jl b/src/PotentialCalculation/PotentialCalculationSetup/BoundaryConditions/BoundaryConditionsCylindrical.jl
index c9e4a4b83..6b3940e5b 100644
--- a/src/PotentialCalculation/PotentialCalculationSetup/BoundaryConditions/BoundaryConditionsCylindrical.jl
+++ b/src/PotentialCalculation/PotentialCalculationSetup/BoundaryConditions/BoundaryConditionsCylindrical.jl
@@ -30,8 +30,60 @@ function apply_boundary_conditions_on_r_axis!(  rbpot::AbstractArray{T, 4}, rbi:
     nothing
 end
 
+# One work-item per z-row: sum the row's φ columns of matching parity with the
+# precomputed weights and broadcast the result back to those columns. Replaces
+# a getindex + gemv + broadcast sequence (a temporary and several kernel
+# launches per SOR iteration) on GPU backends.
+@kernel function r0_bc_kernel!(rbpot, rbi::Int, nnz::Int, gw, even_points::Bool)
+    irbz = @index(Global, Linear)
+    rbz = irbz + 1
+    inφ_ub = size(rbpot, 2) - 1
+    @inbounds begin
+        last_odd = nnz - (iseven(nnz) ? 1 : 0)   # odd real-z rows
+        if rbz <= div(last_odd, 2) + mod(last_odd, 2) + 1
+            gwrow = even_points ? 6 : 5
+            φ0 = even_points ? 2 : 1
+            s = zero(eltype(rbpot))
+            j = 0
+            for inφ in φ0:2:inφ_ub
+                j += 1
+                s += rbpot[rbz, inφ + 1, 2, rbi] * gw[gwrow, j]
+            end
+            for inφ in φ0:2:inφ_ub
+                rbpot[rbz, inφ + 1, 2, rbi] = s
+            end
+        end
+        last_even = nnz - (isodd(nnz) ? 1 : 0)   # even real-z rows
+        if rbz <= div(last_even, 2) + mod(last_even, 2) + 1
+            gwrow = even_points ? 5 : 6
+            φ0 = even_points ? 1 : 2
+            s = zero(eltype(rbpot))
+            j = 0
+            for inφ in φ0:2:inφ_ub
+                j += 1
+                s += rbpot[rbz, inφ + 1, 2, rbi] * gw[gwrow, j]
+            end
+            for inφ in φ0:2:inφ_ub
+                rbpot[rbz, inφ + 1, 2, rbi] = s
+            end
+        end
+    end
+end
+
+function apply_boundary_conditions_at_r0!(
+    rbpot::GPUArrays.AbstractGPUArray{T, 4},
+    rbi::Int,
+    nnz::Int,
+    gw::AbstractArray{T},
+    update_even_points::Val{even_points}
+    ) where {T, even_points}
+    kern = r0_bc_kernel!(_ka_get_backend(rbpot))
+    kern(rbpot, rbi, nnz, gw, even_points, ndrange = rbidx(nnz) - 1)
+    nothing
+end
+
 function apply_boundary_conditions_at_r0!(
-    rbpot::AbstractArray{T, 4}, 
+    rbpot::AbstractArray{T, 4},
     rbi::Int,
     nnz::Int,
     gw::AbstractArray{T},
diff --git a/src/PotentialCalculation/PotentialCalculationSetup/BoundaryConditions/BoundaryConditionsGPU.jl b/src/PotentialCalculation/PotentialCalculationSetup/BoundaryConditions/BoundaryConditionsGPU.jl
new file mode 100644
index 000000000..2feeadefa
--- /dev/null
+++ b/src/PotentialCalculation/PotentialCalculationSetup/BoundaryConditions/BoundaryConditionsGPU.jl
@@ -0,0 +1,123 @@
+# GPU boundary-condition path: the CPU implementation applies the ghost-face
+# updates as ~10 view broadcasts per red-black color, whose per-call host
+# overhead dominates the GPU iteration loop. Here each face update is encoded
+# as (active, source index, factor) with ghost = factor * pot[source] and all
+# six faces are applied by one kernel launch per color.
+#
+# The op tables below are 1:1 transcriptions of the broadcast methods in
+# BoundaryConditionsCylindrical.jl / BoundaryConditionsCartesian.jl.
+# Ghost entries with more than one ghost coordinate are never read by the SOR
+# stencil, so the faces may be applied concurrently.
+
+const _GhostFaceOp{T} = Tuple{Bool, Int, T}
+_no_ghost_op(::Type{T}) where {T} = (false, 1, one(T))
+
+# faces of the compressed (red-black) dimension: neighbors sit at rows 2 / end-1
+_ghost_ops_compressed(::DiscreteAxis{T, :fixed, :fixed}, n::Int, gbf) where {T} =
+    (_no_ghost_op(T), _no_ghost_op(T))
+_ghost_ops_compressed(::DiscreteAxis{T, :infinite, :infinite}, n::Int, gbf) where {T} =
+    ((true, 2, gbf[1]), (true, n - 1, gbf[2]))
+_ghost_ops_compressed(::DiscreteAxis{T, :fixed, :infinite}, n::Int, gbf) where {T} =
+    (_no_ghost_op(T), (true, n - 1, gbf[2]))
+_ghost_ops_compressed(::DiscreteAxis{T, :infinite, :fixed}, n::Int, gbf) where {T} =
+    ((true, 2, gbf[1]), _no_ghost_op(T))
+_ghost_ops_compressed(::DiscreteAxis{T, :infinite, :reflecting}, n::Int, gbf) where {T} =
+    ((true, 2, gbf[1]), (true, n - 1, one(T)))
+_ghost_ops_compressed(::DiscreteAxis{T, :reflecting, :infinite}, n::Int, gbf) where {T} =
+    ((true, 2, one(T)), (true, n - 1, gbf[2]))
+_ghost_ops_compressed(::DiscreteAxis{T, :reflecting, :reflecting}, n::Int, gbf) where {T} =
+    ((true, 2, one(T)), (true, n - 1, one(T)))
+_ghost_ops_compressed(::DiscreteAxis{T, :fixed, :reflecting}, n::Int, gbf) where {T} =
+    (_no_ghost_op(T), (true, n - 1, one(T)))
+_ghost_ops_compressed(::DiscreteAxis{T, :reflecting, :fixed}, n::Int, gbf) where {T} =
+    ((true, 2, one(T)), _no_ghost_op(T))
+_ghost_ops_compressed(::DiscreteAxis{T, :periodic, :periodic}, n::Int, gbf) where {T} =
+    ((true, n - 1, one(T)), (true, 2, one(T)))
+
+# faces of the non-compressed dimensions: reflecting neighbors sit at rows 3 / end-2
+_ghost_ops(::DiscreteAxis{T, :fixed, :fixed}, n::Int, gbf) where {T} =
+    (_no_ghost_op(T), _no_ghost_op(T))
+_ghost_ops(::DiscreteAxis{T, :infinite, :infinite}, n::Int, gbf) where {T} =
+    ((true, 2, gbf[1]), (true, n - 1, gbf[2]))
+_ghost_ops(::DiscreteAxis{T, :fixed, :infinite}, n::Int, gbf) where {T} =
+    (_no_ghost_op(T), (true, n - 1, gbf[2]))
+_ghost_ops(::DiscreteAxis{T, :infinite, :fixed}, n::Int, gbf) where {T} =
+    ((true, 2, gbf[1]), _no_ghost_op(T))
+_ghost_ops(::DiscreteAxis{T, :infinite, :reflecting}, n::Int, gbf) where {T} =
+    ((true, 2, gbf[1]), (true, n - 2, one(T)))
+_ghost_ops(::DiscreteAxis{T, :reflecting, :infinite}, n::Int, gbf) where {T} =
+    ((true, 3, one(T)), (true, n - 1, gbf[2]))
+_ghost_ops(::DiscreteAxis{T, :reflecting, :reflecting}, n::Int, gbf) where {T} =
+    ((true, 3, one(T)), (true, n - 2, one(T)))
+_ghost_ops(::DiscreteAxis{T, :fixed, :reflecting}, n::Int, gbf) where {T} =
+    (_no_ghost_op(T), (true, n - 2, one(T)))
+_ghost_ops(::DiscreteAxis{T, :reflecting, :fixed}, n::Int, gbf) where {T} =
+    ((true, 3, one(T)), _no_ghost_op(T))
+_ghost_ops(::DiscreteAxis{T, :periodic, :periodic}, n::Int, gbf) where {T} =
+    ((true, n - 1, one(T)), (true, 2, one(T)))
+# the radial axis: no left face (r = 0 is handled separately), and the
+# infinite/reflecting right face reads row end-2 (as in the broadcast methods)
+_ghost_ops(::DiscreteAxis{T, :r0, :infinite}, n::Int, gbf) where {T} =
+    (_no_ghost_op(T), (true, n - 2, gbf[2]))
+_ghost_ops(::DiscreteAxis{T, :r0, :reflecting}, n::Int, gbf) where {T} =
+    (_no_ghost_op(T), (true, n - 2, one(T)))
+_ghost_ops(::DiscreteAxis{T, :r0, :fixed}, n::Int, gbf) where {T} =
+    (_no_ghost_op(T), _no_ghost_op(T))
+
+@kernel function faces_bc_kernel!(rbpot, rbi::Int,
+        a1l::Bool, s1l::Int, f1l, a1r::Bool, s1r::Int, f1r,
+        a2l::Bool, s2l::Int, f2l, a2r::Bool, s2r::Int, f2r,
+        a3l::Bool, s3l::Int, f3l, a3r::Bool, s3r::Int, f3r)
+    i, j = @index(Global, NTuple)
+    n1 = size(rbpot, 1); n2 = size(rbpot, 2); n3 = size(rbpot, 3)
+    @inbounds begin
+        if i <= n1 && j <= n2
+            if a3l rbpot[i, j,  1, rbi] = f3l * rbpot[i, j, s3l, rbi] end
+            if a3r rbpot[i, j, n3, rbi] = f3r * rbpot[i, j, s3r, rbi] end
+        end
+        if i <= n1 && j <= n3
+            if a2l rbpot[i,  1, j, rbi] = f2l * rbpot[i, s2l, j, rbi] end
+            if a2r rbpot[i, n2, j, rbi] = f2r * rbpot[i, s2r, j, rbi] end
+        end
+        if i <= n2 && j <= n3
+            if a1l rbpot[ 1, i, j, rbi] = f1l * rbpot[s1l, i, j, rbi] end
+            if a1r rbpot[n1, i, j, rbi] = f1r * rbpot[s1r, i, j, rbi] end
+        end
+    end
+end
+
+function _apply_faces_bc!(rbpot::GPUArrays.AbstractGPUArray{T, 4}, rbi::Int,
+        ops1::NTuple{2, _GhostFaceOp{T}}, ops2::NTuple{2, _GhostFaceOp{T}}, ops3::NTuple{2, _GhostFaceOp{T}}) where {T}
+    kern = faces_bc_kernel!(_ka_get_backend(rbpot))
+    n1, n2, n3 = size(rbpot, 1), size(rbpot, 2), size(rbpot, 3)
+    kern(rbpot, rbi,
+        ops1[1]..., ops1[2]..., ops2[1]..., ops2[2]..., ops3[1]..., ops3[2]...,
+        ndrange = (max(n1, n2), max(n2, n3)))
+    nothing
+end
+
+function apply_boundary_conditions!(pcs::PotentialCalculationSetup{T, Cylindrical, 3, <:GPUArrays.AbstractGPUArray},
+        update_even_points::Val{even_points}, only2d::Val{only_2d}) where {T, even_points, only_2d}
+    rbi::Int = even_points ? rb_even::Int : rb_odd::Int
+    rbpot = pcs.potential
+    if !only_2d
+        apply_boundary_conditions_at_r0!(rbpot, rbi, size(pcs.ϵ_r, 3) - 1, pcs.geom_weights[2], update_even_points)
+    end
+    # red-black dimension order is (z, φ, r)
+    ops_z = _ghost_ops_compressed(pcs.grid.axes[3], size(rbpot, 1), pcs.grid_boundary_factors[3])
+    ops_φ = only_2d ? (_no_ghost_op(T), _no_ghost_op(T)) : _ghost_ops(pcs.grid.axes[2], size(rbpot, 2), pcs.grid_boundary_factors[2])
+    ops_r = _ghost_ops(pcs.grid.axes[1], size(rbpot, 3), pcs.grid_boundary_factors[1])
+    _apply_faces_bc!(rbpot, rbi, ops_z, ops_φ, ops_r)
+    nothing
+end
+
+function apply_boundary_conditions!(pcs::PotentialCalculationSetup{T, Cartesian, 3, <:GPUArrays.AbstractGPUArray},
+        update_even_points::Val{even_points}, only2d::Val{only_2d}) where {T, even_points, only_2d}
+    rbi::Int = even_points ? rb_even::Int : rb_odd::Int
+    rbpot = pcs.potential
+    ops_x = _ghost_ops_compressed(pcs.grid.axes[1], size(rbpot, 1), pcs.grid_boundary_factors[1])
+    ops_y = _ghost_ops(pcs.grid.axes[2], size(rbpot, 2), pcs.grid_boundary_factors[2])
+    ops_z = _ghost_ops(pcs.grid.axes[3], size(rbpot, 3), pcs.grid_boundary_factors[3])
+    _apply_faces_bc!(rbpot, rbi, ops_x, ops_y, ops_z)
+    nothing
+end
diff --git a/src/PotentialCalculation/SuccessiveOverRelaxation/GPU_kernel.jl b/src/PotentialCalculation/SuccessiveOverRelaxation/GPU_kernel.jl
index b2c5c8f84..02f970681 100644
--- a/src/PotentialCalculation/SuccessiveOverRelaxation/GPU_kernel.jl
+++ b/src/PotentialCalculation/SuccessiveOverRelaxation/GPU_kernel.jl
@@ -164,6 +164,11 @@ Maybe this will change in the future.
 """
 get_sor_kernel(::Type{Cylindrical}, backend, CPU_via_KernelAbstractions, args...) = sor_cyl_gpu!(backend, args...)
 get_sor_kernel(::Type{Cartesian},   backend, CPU_via_KernelAbstractions, args...) = sor_car_gpu!(backend, args...)
+
+# Workgroup shapes tuned on an A100; both are 2.4-2.6x faster than the
+# KernelAbstractions default across 2D/3D cylindrical and Cartesian grids.
+# Dimension 1 is the memory-contiguous red-black dimension.
+_sor_kernel_workgroupsize(ndrange) = ndrange[2] == 1 ? (128, 1, 1) : (32, 4, 1)
 get_sor_kernel(::Type{Cylindrical}, ::CPU, CPU_via_KernelAbstractions::Val{false}) = nothing
 get_sor_kernel(::Type{Cartesian},   ::CPU, CPU_via_KernelAbstractions::Val{false}) = nothing
 get_sor_kernel(::Type{Cylindrical}, backend::CPU, CPU_via_KernelAbstractions::Val{true}, args...) = sor_cyl_gpu!(backend, args...)
diff --git a/src/PotentialCalculation/SuccessiveOverRelaxation/GPU_update.jl b/src/PotentialCalculation/SuccessiveOverRelaxation/GPU_update.jl
index 12dfd2f49..0a820e92c 100644
--- a/src/PotentialCalculation/SuccessiveOverRelaxation/GPU_update.jl
+++ b/src/PotentialCalculation/SuccessiveOverRelaxation/GPU_update.jl
@@ -8,14 +8,14 @@
     only2d::Val{only_2d} = Val{false}()
 )::Nothing where {T, only_2d, depletion_handling_enabled, _is_weighting_potential}
     update_even_points = true
-    _ka_synchronize(kernel, kernel( 
+    _ka_synchronize(kernel, kernel(
         pcs.potential, pcs.imp_scale, pcs.point_types, pcs.volume_weights, pcs.q_eff_imp, pcs.q_eff_fix, pcs.ϵ_r,
-        pcs.geom_weights, pcs.sor_const, update_even_points, depletion_handling_enabled, _is_weighting_potential, only_2d, 
+        pcs.geom_weights, pcs.sor_const, update_even_points, depletion_handling_enabled, _is_weighting_potential, only_2d,
         ndrange = ndrange
     ))
     apply_boundary_conditions!(pcs, Val(update_even_points), only2d)
     update_even_points = false
-    _ka_synchronize(kernel, kernel( 
+    _ka_synchronize(kernel, kernel(
         pcs.potential, pcs.imp_scale, pcs.point_types, pcs.volume_weights, pcs.q_eff_imp, pcs.q_eff_fix, pcs.ϵ_r,
         pcs.geom_weights, pcs.sor_const, update_even_points, depletion_handling_enabled, _is_weighting_potential, only_2d,
         ndrange = ndrange
diff --git a/src/PotentialCalculation/SuccessiveOverRelaxation/convergence.jl b/src/PotentialCalculation/SuccessiveOverRelaxation/convergence.jl
index 2ff8be8ab..3e3f06dca 100644
--- a/src/PotentialCalculation/SuccessiveOverRelaxation/convergence.jl
+++ b/src/PotentialCalculation/SuccessiveOverRelaxation/convergence.jl
@@ -11,7 +11,11 @@ function _update_till_convergence!( pcs::PotentialCalculationSetup{T, S, 3},
                                 ) where {T, S, depletion_handling_enabled, only_2d, _is_weighting_potential} #  <: GPUArrays.AbstractGPUArray
     backend = _ka_get_backend(pcs.potential)
     ndrange = size(pcs.potential)[1:3] .- 2
-    kernel = get_sor_kernel(S, backend, Val(via_KernelAbstractions))
+    kernel = if backend isa _KA_Backend && !(backend isa KernelAbstractions.CPU)
+        get_sor_kernel(S, backend, Val(via_KernelAbstractions), _sor_kernel_workgroupsize(ndrange))
+    else
+        get_sor_kernel(S, backend, Val(via_KernelAbstractions))
+    end
     c_limit = _is_weighting_potential ? convergence_limit : abs(convergence_limit * (iszero(pcs.bias_voltage) ? maximum(abs.(pcs.potential)) : pcs.bias_voltage))
     c = (one(c_limit) + c_limit) * 10 # Has to be larger than c_limit at the beginning
     n_performed_iterations = 0
diff --git a/src/ka_compat.jl b/src/ka_compat.jl
index d01f91927..6c9d7a3f0 100644
--- a/src/ka_compat.jl
+++ b/src/ka_compat.jl
@@ -4,7 +4,11 @@
 @static if isdefined(KernelAbstractions, :Backend)
     import KernelAbstractions.Backend as _KA_Backend
     import KernelAbstractions.get_backend as _ka_get_backend
-    _ka_synchronize(kernel, kernel_ret) = KernelAbstractions.synchronize(kernel.backend)
+    # GPU backends are stream-ordered: subsequent operations on the same queue see
+    # the kernel's results without a host-side synchronization. Only the CPU
+    # backend (task-based) needs an explicit wait.
+    _ka_synchronize(kernel, kernel_ret) = kernel.backend isa KernelAbstractions.CPU ?
+        KernelAbstractions.synchronize(kernel.backend) : nothing
 else
     # KernelAbstractions < v0.9:
     import KernelAbstractions.Device as _KA_Backend
diff --git a/test/Project.toml b/test/Project.toml
index 0b7a10ffc..a4aaa2952 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -3,6 +3,7 @@ ArraysOfArrays = "65a8f2f4-9b39-5baf-92e2-a9cc46fdf018"
 Geant4 = "559df036-b7a0-42fd-85df-7d5dd9d70f44"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 InverseFunctions = "3587e190-3f89-42d0-90ee-14403ec27112"
+JLArrays = "27aeb0d3-9eb9-45fb-866b-73c2ecf80fcb"
 LegendDataTypes = "99e09c13-5545-5ee2-bfa2-77f358fb75d8"
 LegendHDF5IO = "c9265ca6-b027-5446-b1a4-febfa8dd10b0"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -22,5 +23,6 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [compat]
 Geant4 = "0.3.1"
+JLArrays = "0.2, 0.3"
 LegendDataTypes = "0.1.3"
 LegendHDF5IO = "0.1"
diff --git a/test/runtests.jl b/test/runtests.jl
index 65668108b..5668aad8f 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -15,6 +15,10 @@ end
     include("test_gpu_backend.jl")
 end
 
+@timed_testset "SOR boundary conditions" begin
+    include("test_boundary_conditions.jl")
+end
+
 @timed_testset "ConstructiveSolidGeometry" begin
     include("ConstructiveSolidGeometry/CSG_IO.jl")
     include("ConstructiveSolidGeometry/CSG_coordinates.jl")
diff --git a/test/test_boundary_conditions.jl b/test/test_boundary_conditions.jl
new file mode 100644
index 000000000..eb8735397
--- /dev/null
+++ b/test/test_boundary_conditions.jl
@@ -0,0 +1,182 @@
+# This file is a part of SolidStateDetectors.jl, licensed under the MIT License (MIT).
+
+using Test
+using SolidStateDetectors
+
+import SolidStateDetectors: DiscreteAxis, Interval, KernelAbstractions, rbidx,
+    apply_boundary_conditions_on_x_axis!, apply_boundary_conditions_on_y_axis!,
+    apply_boundary_conditions_on_z_axis!, apply_boundary_conditions_on_r_axis!,
+    apply_boundary_conditions_on_φ_axis!, apply_boundary_conditions_at_r0!,
+    r0_bc_kernel!, _ghost_ops, _ghost_ops_compressed
+
+T = Float64
+
+# The red-black potential arrays carry one ghost row per face; the methods under
+# test rewrite ghost rows from interior rows of the same parity slab. Distinct
+# entries make every source row uniquely identifiable.
+fresh_rbpot() = reshape(T.(1:(6 * 7 * 8 * 2)), 6, 7, 8, 2)
+
+cc_ax(bl, br) = DiscreteAxis{T, bl, br}(Interval{:closed, :closed, T}(0, 1), T[0, 1])
+co_ax(bl, br) = DiscreteAxis{T, bl, br}(Interval{:closed, :open, T}(0, 1), T[0, 1])
+
+const GBF = (T(0.25), T(0.5))
+
+# expected ghost writes as (side, source_row, factor)
+lo(src, fac = one(T)) = (:lo, src, T(fac))
+hi(src, fac = one(T)) = (:hi, src, T(fac))
+
+function expected_result(A::Array{T, 4}, dim::Int, writes)
+    B = copy(A)
+    slab = view(B, :, :, :, 1)
+    for (side, src, fac) in writes
+        dst = side === :lo ? 1 : size(B, dim)
+        selectdim(slab, dim, dst) .= fac .* selectdim(slab, dim, src)
+    end
+    B
+end
+
+function apply_op_pair!(A::Array{T, 4}, dim::Int, ops)
+    slab = view(A, :, :, :, 1)
+    for (i, (active, src, fac)) in enumerate(ops)
+        active || continue
+        dst = i == 1 ? 1 : size(A, dim)
+        selectdim(slab, dim, dst) .= fac .* selectdim(slab, dim, src)
+    end
+    A
+end
+
+# Ghost-row semantics per (left, right) boundary type. The compressed
+# (red-black) dimension has its interior neighbors at rows 2 / end-1; the
+# non-compressed dimensions mirror at 3 / end-2 for reflecting faces.
+compressed_writes(n) = Dict(
+    (:fixed, :fixed)           => (),
+    (:infinite, :infinite)     => (lo(2, GBF[1]), hi(n - 1, GBF[2])),
+    (:infinite, :fixed)        => (lo(2, GBF[1]),),
+    (:fixed, :infinite)        => (hi(n - 1, GBF[2]),),
+    (:infinite, :reflecting)   => (lo(2, GBF[1]), hi(n - 1)),
+    (:reflecting, :infinite)   => (lo(2), hi(n - 1, GBF[2])),
+    (:reflecting, :reflecting) => (lo(2), hi(n - 1)),
+    (:fixed, :reflecting)      => (hi(n - 1),),
+    (:reflecting, :fixed)      => (lo(2),),
+    (:periodic, :periodic)     => (lo(n - 1), hi(2)),
+)
+noncompressed_writes(n) = Dict(
+    (:fixed, :fixed)           => (),
+    (:infinite, :infinite)     => (lo(2, GBF[1]), hi(n - 1, GBF[2])),
+    (:infinite, :fixed)        => (lo(2, GBF[1]),),
+    (:fixed, :infinite)        => (hi(n - 1, GBF[2]),),
+    (:infinite, :reflecting)   => (lo(2, GBF[1]), hi(n - 2)),
+    (:reflecting, :infinite)   => (lo(3), hi(n - 1, GBF[2])),
+    (:reflecting, :reflecting) => (lo(3), hi(n - 2)),
+    (:fixed, :reflecting)      => (hi(n - 2),),
+    (:reflecting, :fixed)      => (lo(3),),
+    (:periodic, :periodic)     => (lo(n - 1), hi(2)),
+)
+# r-axis outer face reads end-2 even for the infinite decay (whose factor is
+# defined from the end-1 tick); pins current behavior, see issue #617
+r_axis_writes(n) = Dict(
+    (:r0, :infinite)   => (hi(n - 2, GBF[2]),),
+    (:r0, :reflecting) => (hi(n - 2),),
+    (:r0, :fixed)      => (),
+)
+
+axis_for(bl, br) = bl === :periodic ? co_ax(bl, br) : cc_ax(bl, br)
+
+@testset "CPU ghost-row methods per axis and boundary types" begin
+    @testset "compressed dimension (Cartesian x / cylindrical z)" begin
+        for (key, writes) in compressed_writes(6)
+            A = fresh_rbpot()
+            ax = axis_for(key...)
+            apply_boundary_conditions_on_x_axis!(A, 1, ax, ax.interval, GBF)
+            @test A == expected_result(fresh_rbpot(), 1, writes)
+        end
+    end
+    @testset "non-compressed Cartesian y" begin
+        for (key, writes) in noncompressed_writes(7)
+            A = fresh_rbpot()
+            ax = axis_for(key...)
+            apply_boundary_conditions_on_y_axis!(A, 1, ax, ax.interval, GBF)
+            @test A == expected_result(fresh_rbpot(), 2, writes)
+        end
+    end
+    @testset "non-compressed Cartesian z" begin
+        # the mixed infinite/reflecting combinations have no z-axis methods
+        # (see issue #617), so only the implemented ones are pinned here
+        for (key, writes) in noncompressed_writes(8)
+            key in ((:infinite, :reflecting), (:reflecting, :infinite)) && continue
+            A = fresh_rbpot()
+            ax = axis_for(key...)
+            apply_boundary_conditions_on_z_axis!(A, 1, ax, ax.interval, GBF)
+            @test A == expected_result(fresh_rbpot(), 3, writes)
+        end
+    end
+    @testset "cylindrical r" begin
+        for (key, writes) in r_axis_writes(8)
+            A = fresh_rbpot()
+            ax = cc_ax(key...)
+            apply_boundary_conditions_on_r_axis!(A, 1, ax, ax.interval, GBF)
+            @test A == expected_result(fresh_rbpot(), 3, writes)
+        end
+    end
+    @testset "cylindrical φ" begin
+        for (key, writes) in ((:periodic, :periodic) => (lo(6), hi(2)),
+                              (:reflecting, :reflecting) => (lo(3), hi(5)))
+            A = fresh_rbpot()
+            ax = axis_for(key...)
+            apply_boundary_conditions_on_φ_axis!(A, 1, ax, ax.interval)
+            @test A == expected_result(fresh_rbpot(), 2, writes)
+        end
+    end
+end
+
+@testset "GPU ghost-face op tables match the CPU methods" begin
+    for (key, _) in compressed_writes(6)
+        A = fresh_rbpot()
+        ax = axis_for(key...)
+        apply_boundary_conditions_on_x_axis!(A, 1, ax, ax.interval, GBF)
+        B = apply_op_pair!(fresh_rbpot(), 1, _ghost_ops_compressed(ax, 6, GBF))
+        @test A == B
+    end
+    for (key, _) in noncompressed_writes(7)
+        A = fresh_rbpot()
+        ax = axis_for(key...)
+        apply_boundary_conditions_on_y_axis!(A, 1, ax, ax.interval, GBF)
+        B = apply_op_pair!(fresh_rbpot(), 2, _ghost_ops(ax, 7, GBF))
+        @test A == B
+    end
+    for (key, _) in r_axis_writes(8)
+        A = fresh_rbpot()
+        ax = cc_ax(key...)
+        apply_boundary_conditions_on_r_axis!(A, 1, ax, ax.interval, GBF)
+        B = apply_op_pair!(fresh_rbpot(), 3, _ghost_ops(ax, 8, GBF))
+        @test A == B
+    end
+end
+
+@testset "r0 axis kernel matches the CPU broadcast implementation" begin
+    # dim order (z, φ, r): 6 compressed z rows ⇒ up to 8 real z planes,
+    # φ columns 1..5 (plus ghosts), the axis ring lives at r-row 2
+    nnz = 8
+    gw = rand(T, 6, 4)
+    for even in (true, false), rbi in (1, 2)
+        A = rand(T, 6, 7, 8, 2)
+        B = copy(A)
+        apply_boundary_conditions_at_r0!(A, rbi, nnz, gw, Val(even))
+        backend = KernelAbstractions.CPU()
+        r0_bc_kernel!(backend)(B, rbi, nnz, gw, even, ndrange = rbidx(nnz) - 1)
+        KernelAbstractions.synchronize(backend)
+        @test A == B
+    end
+end
+
+@testset "Geometrical axis weights array interface" begin
+    w = rand(T, 4, 5)
+    for W in (SolidStateDetectors.GeometricalCartesianAxisWeights{T},
+              SolidStateDetectors.GeometricalAzimutalAxisWeights{T},
+              SolidStateDetectors.GeometricalRadialAxisWeights{T})
+        gw = W(w)
+        @test size(gw) == size(w)
+        @test gw[2, 3] == w[2, 3]
+        @test all(gw[i] == w[i] for i in eachindex(w))
+    end
+end
diff --git a/test/test_charge_drift_models.jl b/test/test_charge_drift_models.jl
index b4a401e7d..2eb370fe2 100644
--- a/test/test_charge_drift_models.jl
+++ b/test/test_charge_drift_models.jl
@@ -776,3 +776,73 @@ end
     @test isapprox(sum(nb_zero_cyl.energies), ustrip.(u"eV", Edep))
     @test all(x -> all(isfinite, (x.x, x.y, x.z)), nb_zero_cyl.locations)
 end
+
+@timed_testset "ADL2016 zero-field guard" begin
+    cdm2016 = ADL2016ChargeDriftModel(T = T)
+    @test getVe(zero(SVector{3, T}), cdm2016) == zero(SVector{3, T})
+    @test getVh(zero(SVector{3, T}), cdm2016) == zero(SVector{3, T})
+end
+
+@timed_testset "Simple drift models: numeric velocities" begin
+    fv = SVector{3, T}(100, -200, 300)
+    efcdm = ElectricFieldChargeDriftModel{T}()
+    @test getVe(fv, efcdm) == -fv
+    @test getVh(fv, efcdm) == fv
+    icdm = IsotropicChargeDriftModel{T}()
+    @test getVe(fv, icdm) ≈ -icdm.μ_e * fv
+    @test getVh(fv, icdm) ≈ icdm.μ_h * fv
+end
+
+@timed_testset "Temperature model selection and vacuum scaling" begin
+    @test SolidStateDetectors.VacuumTemperatureModel() isa SolidStateDetectors.VacuumTemperatureModel{Float32}
+    @test SolidStateDetectors.VacuumTemperatureModel{Float64}() isa SolidStateDetectors.VacuumTemperatureModel{Float64}
+    # no temperature_dependence entry => no rescaling (vacuum model)
+    @test SolidStateDetectors.TemperatureModel(T, Dict{String, Any}()) isa SolidStateDetectors.VacuumTemperatureModel{T}
+    # deprecated step-level scaling models fall back to vacuum with a warning
+    cfg_deprecated = Dict("temperature_dependence" => Dict("model" => "Boltzmann"))
+    @test (@test_logs (:warn, r"deprecated") SolidStateDetectors.TemperatureModel(T, cfg_deprecated)) isa SolidStateDetectors.VacuumTemperatureModel{T}
+    # unknown model names fall back to vacuum with an info message
+    cfg_unknown = Dict("temperature_dependence" => Dict("model" => "DoesNotExist"))
+    @test (@test_logs (:info, r"predefined temperature models") SolidStateDetectors.TemperatureModel(T, cfg_unknown)) isa SolidStateDetectors.VacuumTemperatureModel{T}
+    # a vacuum-model drift model is returned unchanged by temperature scaling
+    adl_cfg = SolidStateDetectors.parse_config_file(SolidStateDetectors.default_ADL_config_file)
+    delete!(adl_cfg["temperature_dependence"], "model")
+    cdm_vac = ADLChargeDriftModel(adl_cfg, T = T)
+    @test cdm_vac.temperaturemodel isa SolidStateDetectors.VacuumTemperatureModel{T}
+    @test scale_to_temperature(cdm_vac, T(300)) === cdm_vac
+    @test scale_to_temperature(cdm_vac, missing) === cdm_vac
+end
+
+# user hook for ArbitraryDriftModificationVolume, dispatched on the config id
+function SolidStateDetectors.modulate_driftvector(sv::CartesianVector{TV}, pt::CartesianPoint{TV},
+        tl::SolidStateDetectors.ArbitraryDriftModificationVolume{TV}, ::Type{Val{1042}}) where {TV <: SolidStateDetectors.SSDFloat}
+    sv / 2
+end
+
+@timed_testset "Virtual drift volumes" begin
+    cfg = SolidStateDetectors.parse_config_file(SSD_examples[:CGD])
+    cfg["detectors"][1]["virtual_drift_volumes"] = [
+        Dict("name" => "dead region", "model" => "dead",
+             "geometry" => Dict("box" => Dict("widths" => [2, 2, 2], "origin" => [0, 0, 0])))
+    ]
+    sim = Simulation{T}(cfg)
+    dv = sim.detector.virtual_drift_volumes[1]
+    @test dv isa SolidStateDetectors.DeadVolume{T}
+    origin = CartesianPoint{T}(0, 0, 0)
+    @test origin in dv
+    @test !(CartesianPoint{T}(0.1, 0, 0) in dv)
+    sv = CartesianVector{T}(1, 2, 3)
+    @test SolidStateDetectors.modulate_driftvector(sv, origin, dv) == CartesianVector{T}(0, 0, 0)
+
+    cfg["detectors"][1]["virtual_drift_volumes"] = [
+        Dict("model" => "arbitrary", "id" => 1042,
+             "geometry" => Dict("box" => Dict("widths" => [2, 2, 2], "origin" => [0, 0, 0])))
+    ]
+    sim2 = Simulation{T}(cfg)
+    av = sim2.detector.virtual_drift_volumes[1]
+    @test av isa SolidStateDetectors.ArbitraryDriftModificationVolume{T}
+    @test SolidStateDetectors.modulate_driftvector(sv, origin, av) == sv / 2
+    # without a user-defined method for the id, modulation must fail loudly
+    av_nohook = SolidStateDetectors.ArbitraryDriftModificationVolume{T}(av.name, 1043, av.geometry)
+    @test_throws ErrorException SolidStateDetectors.modulate_driftvector(sv, origin, av_nohook)
+end
diff --git a/test/test_gpu_backend.jl b/test/test_gpu_backend.jl
index 15cc2a2b8..4d64a7675 100644
--- a/test/test_gpu_backend.jl
+++ b/test/test_gpu_backend.jl
@@ -2,7 +2,16 @@
 
 using Test
 using SolidStateDetectors
+using Unitful
+using JLArrays
 
+T = Float32
+
+# pseudoGPUArray flips `via_KernelAbstractions` so the SOR kernel runs through
+# KernelAbstractions, but as an abstract type it survives Adapt unchanged: the
+# setup arrays stay `Array`s and none of the `AbstractGPUArray`-dispatched code
+# (fused boundary-condition kernels, r0 kernel, device convergence checks) runs.
+# The JLArray testsets below cover those paths with a real GPUArray on CPU memory.
 abstract type pseudoGPUArray <: SolidStateDetectors.GPUArrays.AbstractGPUArray{T where T, N where N} end
 
 @testset "Simulate example detector: Inverted Coax" begin
@@ -23,16 +32,40 @@ end
 @testset "Simulate example detector: CGD" begin
     sim = Simulation{T}(SSD_examples[:CGD])
 
-    calculate_electric_potential!( 
-        sim, 
-        device_array_type = pseudoGPUArray, 
+    calculate_electric_potential!(
+        sim,
+        device_array_type = pseudoGPUArray,
         convergence_limit = 1e-5,
         refinement_limits = [0.2],
         depletion_handling = true
-    )    
+    )
 
     W = SolidStateDetectors.calculate_stored_energy(sim)
     @test isapprox(ustrip(W), 4.39e-6, atol = 1e-5)        
 end
 
+@testset "JLArray backend reproduces the CPU result: Inverted Coax" begin
+    simC = Simulation{T}(SSD_examples[:InvertedCoax])
+    calculate_electric_potential!(simC, convergence_limit = 1e-5,
+        refinement_limits = [0.2, 0.1], depletion_handling = true, verbose = false)
+    simJ = Simulation{T}(SSD_examples[:InvertedCoax])
+    calculate_electric_potential!(simJ, device_array_type = JLArray, convergence_limit = 1e-5,
+        refinement_limits = [0.2, 0.1], depletion_handling = true, verbose = false)
+    @test simJ.electric_potential.grid == simC.electric_potential.grid
+    @test maximum(abs.(simJ.electric_potential.data .- simC.electric_potential.data)) < T(0.01)
+    @test count(simJ.point_types.data .!= simC.point_types.data) <= length(simC.point_types.data) ÷ 10^4
+end
+
+@testset "JLArray backend reproduces the CPU result: CGD" begin
+    simC = Simulation{T}(SSD_examples[:CGD])
+    calculate_electric_potential!(simC, convergence_limit = 1e-5,
+        refinement_limits = [0.2], depletion_handling = true, verbose = false)
+    simJ = Simulation{T}(SSD_examples[:CGD])
+    calculate_electric_potential!(simJ, device_array_type = JLArray, convergence_limit = 1e-5,
+        refinement_limits = [0.2], depletion_handling = true, verbose = false)
+    @test simJ.electric_potential.grid == simC.electric_potential.grid
+    @test maximum(abs.(simJ.electric_potential.data .- simC.electric_potential.data)) < T(0.01)
+    @test count(simJ.point_types.data .!= simC.point_types.data) <= length(simC.point_types.data) ÷ 10^4
+end
+
 
diff --git a/test/test_utils.jl b/test/test_utils.jl
index 8f9dcf88b..67c93dcac 100644
--- a/test/test_utils.jl
+++ b/test/test_utils.jl
@@ -1,9 +1,9 @@
 # This file is a part of SolidStateDetectors.jl, licensed under the MIT License (MIT).
 
 using TimerOutputs
-# using CUDAKernels # Uncomment this line in order to run all tests on (CUDA) GPU
+# using CUDA # Uncomment this line in order to run all tests on (CUDA) GPU
 
-device_array_type = (@isdefined CUDAKernels) ? CUDAKernels.CUDA.CuArray : Array
+device_array_type = (@isdefined CUDA) ? CUDA.CuArray : Array
 
 testtimer() = get_timer("_default_testtimer_")