diff --git a/include/gridtools/fn/backend/common.hpp b/include/gridtools/fn/backend/common.hpp
index 2626eb718..6616ee9ee 100644
--- a/include/gridtools/fn/backend/common.hpp
+++ b/include/gridtools/fn/backend/common.hpp
@@ -31,10 +31,30 @@ namespace gridtools::fn::backend {
                 meta::rename<tuple, Dims>());
         }
 
+        template <class Dims, class Sizes, class UnrollFactors>
+        constexpr GT_FUNCTION auto make_unrolled_loops(Sizes const &sizes, UnrollFactors) {
+            return tuple_util::host_device::fold(
+                [&](auto outer, auto dim) {
+                    using unroll_factor = std::remove_reference_t<decltype(host_device::at_key<decltype(dim)>(
+                        std::declval<UnrollFactors const &>()))>;
+                    return [outer = std::move(outer),
+                               inner = sid::make_unrolled_loop<decltype(dim), unroll_factor::value>(
+                                   host_device::at_key<decltype(dim)>(sizes))](
+                               auto &&...args) { return outer(inner(std::forward<decltype(args)>(args)...)); };
+                },
+                host_device::identity(),
+                meta::rename<tuple, Dims>());
+        }
+
         template <class Sizes>
         constexpr GT_FUNCTION auto make_loops(Sizes const &sizes) {
             return make_loops<get_keys<Sizes>>(sizes);
         }
+
+        template <class Sizes, class UnrollFactors>
+        constexpr GT_FUNCTION auto make_unrolled_loops(Sizes const &sizes, UnrollFactors unroll_factors) {
+            return make_unrolled_loops<get_keys<Sizes>>(sizes, unroll_factors);
+        }
     } // namespace common
 
     template <class T>
diff --git a/include/gridtools/fn/backend/gpu.hpp b/include/gridtools/fn/backend/gpu.hpp
index 5f1500b70..c771941ea 100644
--- a/include/gridtools/fn/backend/gpu.hpp
+++ b/include/gridtools/fn/backend/gpu.hpp
@@ -9,6 +9,7 @@
  */
 #pragma once
 
+#include <algorithm>
 #include <utility>
 
 #include "../../common/cuda_util.hpp"
@@ -23,108 +24,154 @@
 
 namespace gridtools::fn::backend {
     namespace gpu_impl_ {
+        template <class KeyValuePair, class = void>
+        struct is_valid_block_size_key_value_pair : std::false_type {};
+
+        template <template <class...> class List, class Key, class Value>
+        struct is_valid_block_size_key_value_pair<List<Key, Value>,
+            std::enable_if_t<is_integral_constant<Value>::value && (Value::value > 0)>> : std::true_type {};
+
+        template <class BlockSizes>
+        using is_valid_block_sizes =
+            std::bool_constant<meta::is_list<BlockSizes>::value &&
+                               meta::all<meta::transform<is_valid_block_size_key_value_pair, BlockSizes>>::value>;
+
         /*
-         * BlockSizes must be a meta map, mapping dimensions to integral constant block sizes.
+         * ThreadBlockSizes and LoopBlockSizes must be meta maps, mapping dimensions to integral constant block sizes.
          *
          * For example, meta::list<meta::list<dim::i, integral_constant<int, 32>>,
          *                         meta::list<dim::j, integral_constant<int, 8>>,
          *                         meta::list<dim::k, integral_constant<int, 1>>>;
          * When using a cartesian grid.
          */
-        template <class BlockSizes>
+        template <class ThreadBlockSizes, class LoopBlockSizes = meta::list<>>
         struct gpu {
-            using block_sizes_t = BlockSizes;
+            using thread_block_sizes_t = ThreadBlockSizes;
+            using loop_block_sizes_t = LoopBlockSizes;
             cudaStream_t stream = 0;
+
+            static_assert(is_valid_block_sizes<ThreadBlockSizes>::value, "invalid thread block sizes");
+            static_assert(is_valid_block_sizes<LoopBlockSizes>::value, "invalid loop block sizes");
         };
 
-        template <class BlockSizes, class Dims, int I>
-        using block_size_at_dim = meta::second<meta::mp_find<BlockSizes, meta::at_c<Dims, I>>>;
+        template <class BlockSizes>
+        struct block_size_at_dim {
+            template <class Dim>
+            using apply = std::conditional_t<std::is_void_v<meta::mp_find<BlockSizes, Dim, void>>,
+                meta::list<Dim, integral_constant<int, 1>>,
+                meta::mp_find<BlockSizes, Dim, void>>;
+        };
 
         template <class BlockSizes, class Sizes>
-        GT_FUNCTION_DEVICE auto global_thread_index() {
-            using all_keys_t = get_keys<Sizes>;
-            using ndims_t = meta::length<all_keys_t>;
-            using keys_t = meta::rename<hymap::keys, meta::take_c<std::min(3, (int)ndims_t::value), all_keys_t>>;
-            if constexpr (ndims_t::value == 0) {
-                return hymap::keys<>::values<>();
-            } else if constexpr (ndims_t::value == 1) {
-                using block_dim_x = block_size_at_dim<BlockSizes, keys_t, 0>;
-                using values_t = typename keys_t::template values<int>;
-                return values_t(blockIdx.x * block_dim_x::value + threadIdx.x);
-            } else if constexpr (ndims_t::value == 2) {
-                using block_dim_x = block_size_at_dim<BlockSizes, keys_t, 0>;
-                using block_dim_y = block_size_at_dim<BlockSizes, keys_t, 1>;
-                using values_t = typename keys_t::template values<int, int>;
-                return values_t(
-                    blockIdx.x * block_dim_x::value + threadIdx.x, blockIdx.y * block_dim_y::value + threadIdx.y);
-            } else {
-                using block_dim_x = block_size_at_dim<BlockSizes, keys_t, 0>;
-                using block_dim_y = block_size_at_dim<BlockSizes, keys_t, 1>;
-                using block_dim_z = block_size_at_dim<BlockSizes, keys_t, 2>;
-                using values_t = typename keys_t::template values<int, int, int>;
-                return values_t(blockIdx.x * block_dim_x::value + threadIdx.x,
-                    blockIdx.y * block_dim_y::value + threadIdx.y,
-                    blockIdx.z * block_dim_z::value + threadIdx.z);
+        using block_sizes_for_sizes =
+            hymap::from_meta_map<meta::transform<block_size_at_dim<BlockSizes>::template apply, get_keys<Sizes>>>;
+
+        struct global_thread_index_f {
+            template <std::size_t I, class Index>
+            GT_FUNCTION_DEVICE static constexpr int index_at_dim(Index const &idx) {
+                static_assert(I < 3);
+                if constexpr (I == 0)
+                    return idx.x;
+                if constexpr (I == 1)
+                    return idx.y;
+                return idx.z;
             }
-            // disable incorrect warning "missing return statement at end of non-void function"
-            GT_NVCC_DIAG_PUSH_SUPPRESS(940)
-        }
-        GT_NVCC_DIAG_POP_SUPPRESS(940)
 
-        template <class Key>
-        struct at_generator_f {
-            template <class Value>
-            GT_FUNCTION_DEVICE decltype(auto) operator()(Value &&value) const {
-                return device::at_key<Key>(std::forward<Value>(value));
+            template <std::size_t I, class ThreadBlockSize, class LoopBlockSize>
+            GT_FUNCTION_DEVICE constexpr auto operator()(ThreadBlockSize, LoopBlockSize) const {
+                if constexpr (I < 3) {
+                    return index_at_dim<I>(blockIdx) * (ThreadBlockSize::value * LoopBlockSize::value) +
+                           index_at_dim<I>(threadIdx) * LoopBlockSize::value;
+                } else {
+                    return integral_constant<int, 0>();
+                }
+                // disable incorrect warning "missing return statement at end of non-void function"
+                GT_NVCC_DIAG_PUSH_SUPPRESS(940)
+            }
+            GT_NVCC_DIAG_POP_SUPPRESS(940)
+        };
+
+        struct block_size_f {
+            template <std::size_t I, class GlobalThreadIndex, class LoopBlockSize, class Size>
+            GT_FUNCTION_DEVICE constexpr auto operator()(
+                GlobalThreadIndex global_thread_index, LoopBlockSize, Size size) const {
+                if constexpr (I < 3) {
+                    if constexpr (LoopBlockSize::value == 1)
+                        return integral_constant<int, 1>();
+                    else
+                        return std::clamp(size - global_thread_index, 0, int(LoopBlockSize::value));
+                } else {
+                    return size;
+                }
+                // disable incorrect warning "missing return statement at end of non-void function"
+                GT_NVCC_DIAG_PUSH_SUPPRESS(940)
             }
+            GT_NVCC_DIAG_POP_SUPPRESS(940)
         };
 
-        template <class Index, class Sizes>
-        GT_FUNCTION_DEVICE bool in_domain(Index const &index, Sizes const &sizes) {
-            using sizes_t = meta::rename<tuple, Index>;
-            using generators_t = meta::transform<at_generator_f, get_keys<Index>>;
-            auto indexed_sizes = tuple_util::device::generate<generators_t, sizes_t>(sizes);
-            return tuple_util::device::all_of(std::less(), index, indexed_sizes);
+        template <class ThreadBlockSizes, class LoopBlockSizes, class Sizes>
+        GT_FUNCTION_DEVICE constexpr auto global_thread_index(Sizes const &sizes) {
+            using keys_t = meta::rename<hymap::keys, get_keys<Sizes>>;
+
+            constexpr auto thread_block_sizes = block_sizes_for_sizes<ThreadBlockSizes, Sizes>();
+            constexpr auto loop_block_sizes = block_sizes_for_sizes<LoopBlockSizes, Sizes>();
+
+            auto global_thread_indices =
+                tuple_util::device::transform_index(global_thread_index_f{}, thread_block_sizes, loop_block_sizes);
+            auto block_sizes =
+                tuple_util::device::transform_index(block_size_f{}, global_thread_indices, loop_block_sizes, sizes);
+            return std::make_tuple(
+                std::move(global_thread_indices), std::move(block_sizes), std::move(loop_block_sizes));
+        }
+
+        template <class Int, Int... i>
+        constexpr int iseq_product(std::integer_sequence<Int, i...>) {
+            return (1 * ... * i);
         }
 
-        template <class BlockSizes,
+        template <class ThreadBlockSizes,
+            class LoopBlockSizes,
             class Sizes,
             class PtrHolder,
             class Strides,
             class Fun,
-            class NDims = tuple_util::size<Sizes>,
-            class SizeKeys = get_keys<Sizes>>
-        __global__ void kernel(Sizes sizes, PtrHolder ptr_holder, Strides strides, Fun fun) {
-            auto thread_idx = global_thread_index<BlockSizes, Sizes>();
-            if (!in_domain(thread_idx, sizes))
+            int NumThreads = iseq_product(meta::list_to_iseq<block_sizes_for_sizes<ThreadBlockSizes, Sizes>>())>
+        __global__ void __launch_bounds__(NumThreads)
+            kernel(Sizes sizes, PtrHolder ptr_holder, Strides strides, Fun fun) {
+            auto const [thread_idx, block_size, max_block_size] =
+                global_thread_index<ThreadBlockSizes, LoopBlockSizes>(sizes);
+            if (!tuple_util::device::all_of(std::less(), thread_idx, sizes))
                 return;
             auto ptr = ptr_holder();
             sid::multi_shift(ptr, strides, thread_idx);
-            if constexpr (NDims::value <= 3) {
-                fun(ptr, strides);
-            } else {
-                using loop_dims_t = meta::drop_front_c<3, SizeKeys>;
-                common::make_loops<loop_dims_t>(sizes)(std::move(fun))(ptr, strides);
-            }
+            common::make_unrolled_loops(block_size, max_block_size)(std::move(fun))(ptr, strides);
         }
 
-        template <class BlockSizes, class Sizes>
+        template <class ThreadBlockSizes, class LoopBlockSizes, class Sizes>
         std::tuple<dim3, dim3> blocks_and_threads(Sizes const &sizes) {
             using keys_t = get_keys<Sizes>;
             using ndims_t = meta::length<keys_t>;
+            [[maybe_unused]] constexpr auto thread_block_sizes = block_sizes_for_sizes<ThreadBlockSizes, Sizes>();
+            [[maybe_unused]] constexpr auto loop_block_sizes = block_sizes_for_sizes<LoopBlockSizes, Sizes>();
             dim3 blocks(1, 1, 1);
             dim3 threads(1, 1, 1);
             if constexpr (ndims_t::value >= 1) {
-                threads.x = block_size_at_dim<BlockSizes, keys_t, 0>();
-                blocks.x = (tuple_util::get<0>(sizes) + threads.x - 1) / threads.x;
+                threads.x = tuple_util::get<0>(thread_block_sizes);
+                constexpr int block_dim_x =
+                    tuple_util::get<0>(thread_block_sizes) * tuple_util::get<0>(loop_block_sizes);
+                blocks.x = (tuple_util::get<0>(sizes) + block_dim_x - 1) / block_dim_x;
             }
             if constexpr (ndims_t::value >= 2) {
-                threads.y = block_size_at_dim<BlockSizes, keys_t, 1>();
-                blocks.y = (tuple_util::get<1>(sizes) + threads.y - 1) / threads.y;
+                threads.y = tuple_util::get<1>(thread_block_sizes);
+                constexpr int block_dim_y =
+                    tuple_util::get<1>(thread_block_sizes) * tuple_util::get<1>(loop_block_sizes);
+                blocks.y = (tuple_util::get<1>(sizes) + block_dim_y - 1) / block_dim_y;
             }
             if constexpr (ndims_t::value >= 3) {
-                threads.z = block_size_at_dim<BlockSizes, keys_t, 2>();
-                blocks.z = (tuple_util::get<2>(sizes) + threads.z - 1) / threads.z;
+                threads.z = tuple_util::get<2>(thread_block_sizes);
+                constexpr int block_dim_z =
+                    tuple_util::get<2>(thread_block_sizes) * tuple_util::get<2>(loop_block_sizes);
+                blocks.z = (tuple_util::get<2>(sizes) + block_dim_z - 1) / block_dim_z;
             }
             return {blocks, threads};
         }
@@ -134,7 +181,7 @@ namespace gridtools::fn::backend {
             MakeIterator m_make_iterator;
 
             template <class Ptr, class Strides>
-            GT_FUNCTION_DEVICE void operator()(Ptr &ptr, Strides const &strides) const {
+            GT_FUNCTION_DEVICE constexpr void operator()(Ptr &ptr, Strides const &strides) const {
                 StencilStage()(m_make_iterator(), ptr, strides);
             }
         };
@@ -144,8 +191,13 @@ namespace gridtools::fn::backend {
             return tuple_util::host::apply([](auto... sizes) { return ((sizes == 0) || ...); }, sizes);
         }
 
-        template <class BlockSizes, class Sizes, class StencilStage, class MakeIterator, class Composite>
-        void apply_stencil_stage(gpu<BlockSizes> const &g,
+        template <class ThreadBlockSizes,
+            class LoopBlockSizes,
+            class Sizes,
+            class StencilStage,
+            class MakeIterator,
+            class Composite>
+        void apply_stencil_stage(gpu<ThreadBlockSizes, LoopBlockSizes> const &g,
             Sizes const &sizes,
             StencilStage,
             MakeIterator make_iterator,
@@ -158,13 +210,14 @@ namespace gridtools::fn::backend {
             auto ptr_holder = sid::get_origin(std::forward<Composite>(composite));
             auto strides = sid::get_strides(std::forward<Composite>(composite));
 
-            auto [blocks, threads] = blocks_and_threads<BlockSizes>(sizes);
+            auto [blocks, threads] = blocks_and_threads<ThreadBlockSizes, LoopBlockSizes>(sizes);
             assert(threads.x > 0 && threads.y > 0 && threads.z > 0);
             cuda_util::launch(blocks,
                 threads,
                 0,
                 g.stream,
-                kernel<BlockSizes,
+                kernel<ThreadBlockSizes,
+                    LoopBlockSizes,
                     Sizes,
                     decltype(ptr_holder),
                     decltype(strides),
@@ -187,14 +240,15 @@ namespace gridtools::fn::backend {
             }
         };
 
-        template <class BlockSizes,
+        template <class ThreadBlockSizes,
+            class LoopBlockSizes,
             class Sizes,
             class ColumnStage,
             class MakeIterator,
             class Composite,
             class Vertical,
             class Seed>
-        void apply_column_stage(gpu<BlockSizes> const &g,
+        void apply_column_stage(gpu<ThreadBlockSizes, LoopBlockSizes> const &g,
             Sizes const &sizes,
             ColumnStage,
             MakeIterator make_iterator,
@@ -211,13 +265,14 @@ namespace gridtools::fn::backend {
             auto h_sizes = hymap::canonicalize_and_remove_key<Vertical>(sizes);
             int v_size = at_key<Vertical>(sizes);
 
-            auto [blocks, threads] = blocks_and_threads<BlockSizes>(h_sizes);
+            auto [blocks, threads] = blocks_and_threads<ThreadBlockSizes, LoopBlockSizes>(h_sizes);
             assert(threads.x > 0 && threads.y > 0 && threads.z > 0);
             cuda_util::launch(blocks,
                 threads,
                 0,
                 g.stream,
-                kernel<BlockSizes,
+                kernel<ThreadBlockSizes,
+                    LoopBlockSizes,
                     decltype(h_sizes),
                     decltype(ptr_holder),
                     decltype(strides),
@@ -228,13 +283,14 @@ namespace gridtools::fn::backend {
                 column_fun_f<ColumnStage, MakeIterator, Seed>{std::move(make_iterator), std::move(seed), v_size});
         }
 
-        template <class BlockSizes>
-        auto tmp_allocator(gpu<BlockSizes> be) {
+        template <class ThreadBlockSizes, class LoopBlockSizes>
+        auto tmp_allocator(gpu<ThreadBlockSizes, LoopBlockSizes> be) {
             return std::make_tuple(be, sid::device::cached_allocator(&cuda_util::cuda_malloc<char[]>));
         }
 
-        template <class BlockSizes, class Allocator, class Sizes, class T>
-        auto allocate_global_tmp(std::tuple<gpu<BlockSizes>, Allocator> &alloc, Sizes const &sizes, data_type<T>) {
+        template <class ThreadBlockSizes, class LoopBlockSizes, class Allocator, class Sizes, class T>
+        auto allocate_global_tmp(
+            std::tuple<gpu<ThreadBlockSizes, LoopBlockSizes>, Allocator> &alloc, Sizes const &sizes, data_type<T>) {
             return sid::make_contiguous<T, int_t, sid::unknown_kind>(std::get<1>(alloc), sizes);
         }
     } // namespace gpu_impl_
diff --git a/include/gridtools/sid/loop.hpp b/include/gridtools/sid/loop.hpp
index 0fe7809ba..fc623cf39 100644
--- a/include/gridtools/sid/loop.hpp
+++ b/include/gridtools/sid/loop.hpp
@@ -15,6 +15,7 @@
 #include <utility>
 
 #include "../common/defs.hpp"
+#include "../common/for_each.hpp"
 #include "../common/functional.hpp"
 #include "../common/host_device.hpp"
 #include "../common/integral_constant.hpp"
@@ -637,6 +638,30 @@ namespace gridtools {
             return {};
         }
 
+        template <class Key, int UnrollFactor, class NumSteps, class Step = integral_constant<int, 1>>
+        constexpr GT_FUNCTION auto make_unrolled_loop(NumSteps num_steps, Step step = {}) {
+            using u = integral_constant<int, UnrollFactor>;
+            return [step,
+                       unrolled = make_loop<Key>(num_steps / u(), step * u()),
+                       epilogue = make_loop<Key>(num_steps % u(), step),
+                       epilogue_start = step * ((num_steps / u()) * u())](auto &&fun) {
+                return [unrolled = unrolled([step, fun=std::forward<decltype(fun)>(fun)](auto &&ptr, auto const strides) {
+                    ::gridtools::host_device::for_each<meta::make_indices_c<UnrollFactor>>([&](auto) {
+                        fun(std::forward<decltype(ptr)>(ptr), strides);
+                        shift(std::forward<decltype(ptr)>(ptr), get_stride<Key>(strides), step);
+                    });
+                    shift(std::forward<decltype(ptr)>(ptr), get_stride<Key>(strides), -step * u());
+                }),
+                           epilogue = epilogue(std::forward<decltype(fun)>(fun)),
+                           epilogue_start](auto &&ptr, auto const &strides) {
+                    unrolled(std::forward<decltype(ptr)>(ptr), strides);
+                    shift(std::forward<decltype(ptr)>(ptr), get_stride<Key>(strides), epilogue_start);
+                    epilogue(std::forward<decltype(ptr)>(ptr), strides);
+                    shift(std::forward<decltype(ptr)>(ptr), get_stride<Key>(strides), -epilogue_start);
+                };
+            };
+        }
+
         /**
          *   A helper that allows to use `SID`s with C++11 range based loop
          *
diff --git a/jenkins/envs/daint_nvcc_cray.sh b/jenkins/envs/daint_nvcc_cray.sh
index 609e81ab5..ae98ddd0e 100644
--- a/jenkins/envs/daint_nvcc_cray.sh
+++ b/jenkins/envs/daint_nvcc_cray.sh
@@ -7,6 +7,7 @@ module switch cudatoolkit/11.2.0_3.39-2.1__gf93aa1c
 module switch cce/10.0.2
 
 export GTCMAKE_GT_CLANG_CUDA_MODE=NVCC-CUDA
+export GTCMAKE_CMAKE_CUDA_FLAGS_RELEASE='-O3 -DNDEBUG --use_fast_math'
 
 export CTEST_PARALLEL_LEVEL=1
 
diff --git a/jenkins/envs/daint_nvcc_gcc.sh b/jenkins/envs/daint_nvcc_gcc.sh
index 713a86125..6598f232b 100644
--- a/jenkins/envs/daint_nvcc_gcc.sh
+++ b/jenkins/envs/daint_nvcc_gcc.sh
@@ -17,5 +17,6 @@ export CUDAHOSTCXX="$CXX"
 
 export GTCMAKE_CMAKE_CXX_FLAGS='-march=haswell'
 export GTCMAKE_CMAKE_CXX_FLAGS_RELEASE='-Ofast -DNDEBUG'
+export GTCMAKE_CMAKE_CUDA_FLAGS_RELEASE='-O3 -DNDEBUG --use_fast_math'
 
 export CTEST_PARALLEL_LEVEL=1
diff --git a/tests/include/fn_select.hpp b/tests/include/fn_select.hpp
index 68cc0ed0b..a207b6ee4 100644
--- a/tests/include/fn_select.hpp
+++ b/tests/include/fn_select.hpp
@@ -47,7 +47,7 @@ namespace {
                                  gridtools::integral_constant>,
             gridtools::meta::list<gridtools::integral_constant<int, sizes>...>>;
 
-    using fn_backend_t = gridtools::fn::backend::gpu<block_sizes_t<32, 8, 1>>;
+    using fn_backend_t = gridtools::fn::backend::gpu<block_sizes_t<32, 8, 1>, block_sizes_t<1, 1, 1>>;
 } // namespace
 #endif
 
@@ -68,17 +68,19 @@ namespace gridtools::fn::backend {
             return "naive";
         }
     } // namespace naive_impl_
+    using naive_impl_::naive_with_threadpool;
 
     namespace gpu_impl_ {
-        template <class>
+        template <class, class>
         struct gpu;
-        template <class BlockSizes>
-        storage::gpu backend_storage_traits(gpu<BlockSizes>);
-        template <class BlockSizes>
-        timer_cuda backend_timer_impl(gpu<BlockSizes>);
-        template <class BlockSizes>
-        inline char const *backend_name(gpu<BlockSizes> const &) {
+        template <class ThreadBlockSizes, class LoopBlockSizes>
+        storage::gpu backend_storage_traits(gpu<ThreadBlockSizes, LoopBlockSizes>);
+        template <class ThreadBlockSizes, class LoopBlockSizes>
+        timer_cuda backend_timer_impl(gpu<ThreadBlockSizes, LoopBlockSizes>);
+        template <class ThreadBlockSizes, class LoopBlockSizes>
+        inline char const *backend_name(gpu<ThreadBlockSizes, LoopBlockSizes> const &) {
             return "gpu";
         }
     } // namespace gpu_impl_
+    using gpu_impl_::gpu;
 } // namespace gridtools::fn::backend
diff --git a/tests/regression/fn/fn_unstructured_nabla.cpp b/tests/regression/fn/fn_unstructured_nabla.cpp
index 04bde8a6a..902da0bbb 100644
--- a/tests/regression/fn/fn_unstructured_nabla.cpp
+++ b/tests/regression/fn/fn_unstructured_nabla.cpp
@@ -184,6 +184,12 @@ namespace {
     static constexpr int vertex_field_id = 0;
     static constexpr int edge_field_id = 1;
 
+    using k_blocked_backend_t = meta::if_<meta::is_instantiation_of<backend::gpu, fn_backend_t>,
+        backend::gpu<meta::list<meta::list<integral_constant<int, 0>, integral_constant<int, 32>>,
+                         meta::list<integral_constant<int, 1>, integral_constant<int, 8>>>,
+            meta::list<meta::list<integral_constant<int, 1>, integral_constant<int, 5>>>>,
+        fn_backend_t>;
+
     constexpr inline auto make_comp = [](auto backend, auto const &mesh, auto &nabla) {
         using mesh_t = std::remove_reference_t<decltype(mesh)>;
         using float_t = typename mesh_t::float_t;
@@ -253,12 +259,12 @@ namespace {
         TypeParam::benchmark("fn_unstructured_nabla_field_of_tuples", comp);
     }
 
-    GT_REGRESSION_TEST(fn_unstructured_nabla_fused_field_of_tuples, test_environment<>, fn_backend_t) {
+    GT_REGRESSION_TEST(fn_unstructured_nabla_fused_field_of_tuples, test_environment<>, k_blocked_backend_t) {
         using float_t = typename TypeParam::float_t;
 
         auto mesh = TypeParam::fn_unstructured_mesh();
         auto nabla = mesh.template make_storage<tuple<float_t, float_t>>(mesh.nvertices(), mesh.nlevels());
-        auto comp = make_comp_fused(fn_backend_t(), mesh, nabla);
+        auto comp = make_comp_fused(k_blocked_backend_t(), mesh, nabla);
         comp();
         auto expected = make_expected(mesh);
         TypeParam::verify(expected, nabla);
@@ -266,6 +272,8 @@ namespace {
     }
 
     GT_REGRESSION_TEST(fn_unstructured_nabla_tuple_of_fields, test_environment<>, fn_backend_t) {
+        using float_t = typename TypeParam::float_t;
+
         auto mesh = TypeParam::fn_unstructured_mesh();
         auto nabla0 = mesh.template make_storage<float_t, vertex_field_id>(mesh.nvertices(), mesh.nlevels());
         auto nabla1 = mesh.template make_storage<float_t, vertex_field_id>(mesh.nvertices(), mesh.nlevels());
@@ -280,14 +288,16 @@ namespace {
         TypeParam::benchmark("fn_unstructured_nabla_tuple_of_fields", comp);
     }
 
-    GT_REGRESSION_TEST(fn_unstructured_nabla_fused_tuple_of_fields, test_environment<>, fn_backend_t) {
+    GT_REGRESSION_TEST(fn_unstructured_nabla_fused_tuple_of_fields, test_environment<>, k_blocked_backend_t) {
+        using float_t = typename TypeParam::float_t;
+
         auto mesh = TypeParam::fn_unstructured_mesh();
         auto nabla0 = mesh.template make_storage<float_t, vertex_field_id>(mesh.nvertices(), mesh.nlevels());
         auto nabla1 = mesh.template make_storage<float_t, vertex_field_id>(mesh.nvertices(), mesh.nlevels());
         auto nabla =
             sid::composite::keys<integral_constant<int, 0>, integral_constant<int, 1>>::make_values(nabla0, nabla1);
 
-        auto comp = make_comp_fused(fn_backend_t(), mesh, nabla);
+        auto comp = make_comp_fused(k_blocked_backend_t(), mesh, nabla);
         comp();
         auto expected = make_expected(mesh);
         TypeParam::verify([&](int vertex, int k) { return get<0>(expected(vertex, k)); }, nabla0);
diff --git a/tests/unit_tests/fn/test_fn_backend_gpu.cu b/tests/unit_tests/fn/test_fn_backend_gpu.cu
index 0f9c45a05..feb472ea2 100644
--- a/tests/unit_tests/fn/test_fn_backend_gpu.cu
+++ b/tests/unit_tests/fn/test_fn_backend_gpu.cu
@@ -63,10 +63,16 @@ namespace gridtools::fn::backend {
 
             column_stage<int_t<1>, sum_scan, 0, 1> cs;
 
-            using block_sizes_t = meta::list<meta::list<int_t<0>, int_t<4>>, meta::list<int_t<2>, int_t<2>>>;
+            using thread_block_sizes_t = meta::list<meta::list<int_t<0>, int_t<4>>, meta::list<int_t<2>, int_t<2>>>;
+            using loop_block_sizes_t = meta::list<meta::list<int_t<0>, int_t<1>>, meta::list<int_t<2>, int_t<2>>>;
 
-            apply_column_stage(
-                gpu<block_sizes_t>(), sizes, cs, make_iterator_mock(), composite, int_t<1>(), tuple(42, 1));
+            apply_column_stage(gpu<thread_block_sizes_t, loop_block_sizes_t>(),
+                sizes,
+                cs,
+                make_iterator_mock(),
+                composite,
+                int_t<1>(),
+                tuple(42, 1));
 
             cudaMemcpy(outh, out.get(), 5 * 7 * 3 * sizeof(int), cudaMemcpyDeviceToHost);
             for (int i = 0; i < 5; ++i)
@@ -100,10 +106,16 @@ namespace gridtools::fn::backend {
 
             column_stage<int_t<0>, sum_scan, 0, 1> cs;
 
-            using block_sizes_t = meta::list<meta::list<int_t<0>, int_t<4>>, meta::list<int_t<2>, int_t<2>>>;
+            using thread_block_sizes_t = meta::list<meta::list<int_t<0>, int_t<4>>, meta::list<int_t<2>, int_t<2>>>;
+            using loop_block_sizes_t = meta::list<meta::list<int_t<0>, int_t<2>>, meta::list<int_t<2>, int_t<2>>>;
 
-            apply_column_stage(
-                gpu<block_sizes_t>(), sizes, cs, make_iterator_mock(), composite, int_t<0>(), tuple(42, 1));
+            apply_column_stage(gpu<thread_block_sizes_t, loop_block_sizes_t>(),
+                sizes,
+                cs,
+                make_iterator_mock(),
+                composite,
+                int_t<0>(),
+                tuple(42, 1));
 
             cudaMemcpy(outh, out.get(), 5 * sizeof(int), cudaMemcpyDeviceToHost);
             int res = 42;
@@ -139,13 +151,22 @@ namespace gridtools::fn::backend {
 
             column_stage<int_t<1>, sum_scan, 0, 1> cs;
 
-            using block_sizes_t = meta::list<meta::list<int_t<0>, int_t<4>>,
+            using thread_block_sizes_t = meta::list<meta::list<int_t<0>, int_t<4>>,
                 meta::list<int_t<2>, int_t<2>>,
                 meta::list<int_t<3>, int_t<2>>,
                 meta::list<int_t<4>, int_t<1>>>;
+            using loop_block_sizes_t = meta::list<meta::list<int_t<0>, int_t<2>>,
+                meta::list<int_t<2>, int_t<2>>,
+                meta::list<int_t<3>, int_t<2>>,
+                meta::list<int_t<4>, int_t<2>>>;
 
-            apply_column_stage(
-                gpu<block_sizes_t>(), sizes, cs, make_iterator_mock(), composite, int_t<1>(), tuple(42, 1));
+            apply_column_stage(gpu<thread_block_sizes_t, loop_block_sizes_t>(),
+                sizes,
+                cs,
+                make_iterator_mock(),
+                composite,
+                int_t<1>(),
+                tuple(42, 1));
 
             cudaMemcpy(outh, out.get(), 5 * 7 * 3 * 2 * 3 * sizeof(int), cudaMemcpyDeviceToHost);
             for (int i = 0; i < 5; ++i)
@@ -197,8 +218,8 @@ namespace gridtools::fn::backend {
         };
 
         TEST(backend_gpu, global_tmp) {
-            using block_sizes_t = meta::list<meta::list<int_t<0>, int_t<4>>, meta::list<int_t<2>, int_t<2>>>;
-            auto alloc = tmp_allocator(gpu<block_sizes_t>());
+            using thread_block_sizes_t = meta::list<meta::list<int_t<0>, int_t<4>>, meta::list<int_t<2>, int_t<2>>>;
+            auto alloc = tmp_allocator(gpu<thread_block_sizes_t>());
             auto sizes = hymap::keys<int_t<0>, int_t<1>, int_t<2>>::values<int_t<5>, int_t<7>, int_t<3>>();
             auto tmp = allocate_global_tmp(alloc, sizes, data_type<int>());
             static_assert(sid::is_sid<decltype(tmp)>());
diff --git a/tests/unit_tests/sid/test_sid_loop.cpp b/tests/unit_tests/sid/test_sid_loop.cpp
index 1ace5b6a4..bd35c11ba 100644
--- a/tests/unit_tests/sid/test_sid_loop.cpp
+++ b/tests/unit_tests/sid/test_sid_loop.cpp
@@ -58,6 +58,35 @@ namespace gridtools {
                     EXPECT_EQ(88, data[i][j]) << " i:" << i << ", j:" << j;
         }
 
+        TEST(make_unrolled_loop, smoke) {
+            double data[10][10] = {};
+            auto strides = tuple(10_c, 1_c);
+
+            double *ptr = &data[0][0];
+            sid::make_unrolled_loop<i_t, 2>(5, 1)(assignment_f{42})(ptr, strides);
+            for (int i = 0; i < 5; ++i)
+                EXPECT_EQ(42, data[i][0]) << " i:" << i;
+
+            ptr = &data[2][3];
+            sid::make_unrolled_loop<j_t, 3>(4_c, -1_c)(assignment_f{5})(ptr, strides);
+            for (int i = 0; i < 4; ++i)
+                EXPECT_EQ(5, data[2][i]) << " i:" << i;
+
+            ptr = &data[0][0];
+            sid::make_unrolled_loop<i_t, 2>(10_c)(sid::make_unrolled_loop<j_t, 7>(10_c)(assignment_f{88}))(
+                ptr, strides);
+            for (int i = 0; i < 10; ++i)
+                for (int j = 0; j < 10; ++j)
+                    EXPECT_EQ(88, data[i][j]) << " i:" << i << ", j:" << j;
+
+            // pass ptr as r-value
+            sid::make_unrolled_loop<i_t, 3>(10_c)(sid::make_unrolled_loop<j_t, 1>(10_c)(assignment_f{88}))(
+                &data[0][0], strides);
+            for (int i = 0; i < 10; ++i)
+                for (int j = 0; j < 10; ++j)
+                    EXPECT_EQ(88, data[i][j]) << " i:" << i << ", j:" << j;
+        }
+
         TEST(nest_loops, smoke) {
             double data[10][10] = {};
             double *ptr = &data[0][0];