Repository: GridTools/gt4py
Pull request #1120 (2023-01-09 16:57:56)

feat[next]: use context vars instead of global state in embedded iterator execution
Changed:
- In the iterator embedded execution, both `column_range` and `offset_provider` global variables have been converted to context variables and fencil closures are now executed in isolated contexts with properly initialized values, to avoid using global state.

havogt: Opening as a reminder, instead of #896 
egparedes: I can clean up the code and provide a final version of this prototype without too much effort. @tehrengruber  do you also agree with using this implementation strategy to fix the problem?
tehrengruber: Yes, given the model itself this is much better than the global variable.
egparedes: I've requested a formal review from @tehrengruber because apparently @havogt can't approve it since he's the creator of the PR, but I think both of you should take a look.

Patch:
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index 9d58e2de..031b5e1e 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -130,7 +130,14 @@ repos:
       tests/cartesian_tests/integration_tests/multi_feature_tests/utils.py |
       tests/cartesian_tests/integration_tests/feature_tests/test_call_interface.py |
       tests/cartesian_tests/unit_tests/frontend_tests/test_gtscript_frontend.py |
-      tests/next_tests/.*
+      tests/next_tests/unit_tests/.* |
+      tests/next_tests/integration_tests/multi_feature_tests/.* |
+      tests/next_tests/integration_tests/feature_tests/ffront_tests/.* |
+      tests/next_tests/integration_tests/feature_tests/iterator_tests/.* |
+      tests/next_tests/integration_tests/feature_tests/otf_tests/.* |
+      tests/next_tests/integration_tests/feature_tests/math_builtin_test_data.py |
+      tests/next_tests/past_common_fixtures.py |
+      tests/next_tests/toy_connectivity.py |
       )$
 
 - repo: https://github.com/pre-commit/mirrors-mypy
@@ -187,5 +194,15 @@ repos:
       build/.* |
       ci/.* |
       docs/.* |
+      tests/conftest.py |
+      tests/cartesian_tests/.* |
+      tests/next_tests/unit_tests/.* |
+      tests/next_tests/integration_tests/multi_feature_tests/.* |
+      tests/next_tests/integration_tests/feature_tests/ffront_tests/.* |
+      tests/next_tests/integration_tests/feature_tests/iterator_tests/.* |
+      tests/next_tests/integration_tests/feature_tests/otf_tests/.* |
+      tests/next_tests/integration_tests/feature_tests/math_builtin_test_data.py |
+      tests/next_tests/past_common_fixtures.py |
+      tests/next_tests/toy_connectivity.py |
       tests/.*
       )$
diff --git a/CODING_GUIDELINES.md b/CODING_GUIDELINES.md
index 60dd7c65..a5b766a8 100644
--- a/CODING_GUIDELINES.md
+++ b/CODING_GUIDELINES.md
@@ -132,6 +132,8 @@ Each test suite should follow the following structure:
 
 Temporarily it may be allowed to split unit tests for a module into multiple `test_<module>_<feature>.py` test modules. This should be taken as a reason to think about splitting the source module though.
 
+Temporarily, tests for testing utilities can be placed next to the module containing them, with the name `test_util_<module>.py`. This should be taken as a hint that the tested utils should be moved into the library.
+
 TODO: add missing test conventions.
 
 <!--
diff --git a/pyproject.toml b/pyproject.toml
index 98734e39..117e62ee 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -224,7 +224,22 @@ skip_glob = ['*.venv/**', '_local/**']
 # -- mypy  --
 [tool.mypy]
 disallow_incomplete_defs = true
-exclude = ['^setup\.py$', '^build/.*py$', '^ci/*.py$', '^tests/*.py$']
+exclude = [
+  '^setup\.py$',
+  '^build/.*py$',
+  '^ci/*.py$',
+  '^tests/cartesian_tests/*.py$',
+  '^tests/storage_tests/*.py$',
+  '^tests/eve_tests/*.py$',
+  '^tests/next_tests/unit_tests/*.py$',
+  '^tests/next_tests/integration_tests/multi_feature_tests/*.py$',
+  '^tests/next_tests/integration_tests/feature_tests/ffront_tests/*.py$',
+  '^tests/next_tests/integration_tests/feature_tests/iterator_tests/*.py$',
+  '^tests/next_tests/integration_tests/feature_tests/otf_tests/*.py$',
+  '^tests/next_tests/integration_tests/feature_tests/math_builtin_test_data.py$',
+  '^tests/next_tests/past_common_fixtures.py$',
+  '^tests/next_tests/toy_connectivity.py$'
+]
 ignore_missing_imports = true
 install_types = true
 namespace_packages = false
diff --git a/src/gt4py/next/common.py b/src/gt4py/next/common.py
index 056330ac..5ece3a23 100644
--- a/src/gt4py/next/common.py
+++ b/src/gt4py/next/common.py
@@ -37,7 +37,7 @@ from gt4py.eve.type_definitions import StrEnum
 
 DimT = TypeVar("DimT", bound="Dimension")
 DimsT = TypeVar("DimsT", bound=Sequence["Dimension"])
-DT = TypeVar("DT", bound="DType")
+DT = TypeVar("DT", bound="Scalar")
 
 Scalar: TypeAlias = SupportsInt | SupportsFloat | np.int32 | np.int64 | np.float32 | np.float64
 
diff --git a/src/gt4py/next/ffront/decorator.py b/src/gt4py/next/ffront/decorator.py
index 632681fc..25d194e1 100644
--- a/src/gt4py/next/ffront/decorator.py
+++ b/src/gt4py/next/ffront/decorator.py
@@ -238,7 +238,7 @@ class Program:
                 f"The following closure variables are undefined: {', '.join(undefined_symbols)}"
             )
 
-    def with_backend(self, backend: ppi.ProgramExecutor) -> "Program":
+    def with_backend(self, backend: ppi.ProgramExecutor) -> Program:
         return Program(
             past_node=self.past_node,
             closure_vars=self.closure_vars,
@@ -246,6 +246,15 @@ class Program:
             definition=self.definition,  # type: ignore[arg-type]  # mypy wrongly deduces definition as method here
         )
 
+    def with_grid_type(self, grid_type: GridType) -> Program:
+        return Program(
+            past_node=self.past_node,
+            closure_vars=self.closure_vars,
+            backend=self.backend,
+            definition=self.definition,
+            grid_type=grid_type,
+        )
+
     @functools.cached_property
     def _all_closure_vars(self) -> dict[str, Any]:
         return _get_closure_vars_recursively(self.closure_vars)
diff --git a/src/gt4py/next/iterator/embedded.py b/src/gt4py/next/iterator/embedded.py
index e3dd18a8..c5befa58 100644
--- a/src/gt4py/next/iterator/embedded.py
+++ b/src/gt4py/next/iterator/embedded.py
@@ -17,6 +17,7 @@
 from __future__ import annotations
 
 import abc
+import contextvars as cvars
 import copy
 import dataclasses
 import itertools
@@ -182,9 +183,10 @@ class MutableLocatedField(LocatedField, Protocol):
         ...
 
 
-# TODO see https://github.com/GridTools/gt4py/pull/1120
-_column_range: Optional[range] = None
-_offset_provider: Optional[OffsetProvider] = None
+#: Column range used in column mode (`column_axis != None`) in the current closure execution context.
+column_range_cvar: cvars.ContextVar[range] = cvars.ContextVar("column_range")
+#: Offset provider dict in the current closure execution context.
+offset_provider_cvar: cvars.ContextVar[OffsetProvider] = cvars.ContextVar("offset_provider")
 
 
 class Column(np.lib.mixins.NDArrayOperatorsMixin):
@@ -197,7 +199,8 @@ class Column(np.lib.mixins.NDArrayOperatorsMixin):
     def __init__(self, kstart: int, data: np.ndarray | Scalar) -> None:
         self.kstart = kstart
         assert isinstance(data, (np.ndarray, Scalar))  # type: ignore # mypy bug
-        self.data = data if isinstance(data, np.ndarray) else np.full(len(_column_range), data)  # type: ignore[arg-type]
+        column_range = column_range_cvar.get()
+        self.data = data if isinstance(data, np.ndarray) else np.full(len(column_range), data)
 
     def __getitem__(self, i: int) -> Any:
         result = self.data[i - self.kstart]
@@ -720,22 +723,23 @@ def _make_tuple(
     *,
     column_axis: Optional[Tag] = None,
 ) -> Column | npt.DTypeLike | tuple[tuple | Column | npt.DTypeLike, ...]:
+    column_range = column_range_cvar.get()
     if isinstance(field_or_tuple, tuple):
         if column_axis is not None:
-            assert _column_range
+            assert column_range
             # construct a Column of tuples
             column_axis_idx = _axis_idx(_get_axes(field_or_tuple), column_axis)
             if column_axis_idx is None:
                 column_axis_idx = -1  # field doesn't have the column index, e.g. ContantField
             first = tuple(
-                _make_tuple(f, _single_vertical_idx(indices, column_axis_idx, _column_range.start))
+                _make_tuple(f, _single_vertical_idx(indices, column_axis_idx, column_range.start))
                 for f in field_or_tuple
             )
             col = Column(
-                _column_range.start, np.zeros(len(_column_range), dtype=_column_dtype(first))
+                column_range.start, np.zeros(len(column_range), dtype=_column_dtype(first))
             )
             col[0] = first
-            for i in _column_range[1:]:
+            for i in column_range[1:]:
                 col[i] = tuple(
                     _make_tuple(f, _single_vertical_idx(indices, column_axis_idx, i))
                     for f in field_or_tuple
@@ -747,8 +751,8 @@ def _make_tuple(
         data = field_or_tuple.field_getitem(indices)
         if column_axis is not None:
             # wraps a vertical slice of an input field into a `Column`
-            assert _column_range is not None
-            return Column(_column_range.start, data)
+            assert column_range is not None
+            return Column(column_range.start, data)
         else:
             return data
 
@@ -768,10 +772,11 @@ class MDIterator:
 
     def shift(self, *offsets: OffsetPart) -> MDIterator:
         complete_offsets = group_offsets(*offsets)
-        assert _offset_provider is not None
+        offset_provider = offset_provider_cvar.get()
+        assert offset_provider is not None
         return MDIterator(
             self.field,
-            shift_position(self.pos, *complete_offsets, offset_provider=_offset_provider),
+            shift_position(self.pos, *complete_offsets, offset_provider=offset_provider),
             column_axis=self.column_axis,
         )
 
@@ -792,13 +797,14 @@ class MDIterator:
             if not all(axis.value in shifted_pos.keys() for axis in axes if axis is not None):
                 raise IndexError("Iterator position doesn't point to valid location for its field.")
         slice_column = dict[Tag, range]()
+        column_range = column_range_cvar.get()
         if self.column_axis is not None:
-            assert _column_range is not None
+            assert column_range is not None
             k_pos = shifted_pos.pop(self.column_axis)
             assert isinstance(k_pos, int)
             # the following range describes a range in the field
             # (negative values are relative to the origin, not relative to the size)
-            slice_column[self.column_axis] = range(k_pos, k_pos + len(_column_range))
+            slice_column[self.column_axis] = range(k_pos, k_pos + len(column_range))
 
         assert _is_concrete_position(shifted_pos)
         ordered_indices = get_ordered_indices(
@@ -834,9 +840,10 @@ def make_in_iterator(
         init = [None] * sparse_dimensions.count(sparse_dim)
         new_pos[sparse_dim] = init  # type: ignore[assignment] # looks like mypy is confused
     if column_axis is not None:
+        column_range = column_range_cvar.get()
         # if we deal with column stencil the column position is just an offset by which the whole column needs to be shifted
-        assert _column_range is not None
-        new_pos[column_axis] = _column_range.start
+        assert column_range is not None
+        new_pos[column_axis] = column_range.start
     it = MDIterator(
         inp,
         new_pos,
@@ -1083,8 +1090,9 @@ class _ConstList(Generic[DT]):
 def neighbors(offset: runtime.Offset, it: ItIterator) -> _List:
     offset_str = offset.value if isinstance(offset, runtime.Offset) else offset
     assert isinstance(offset_str, str)
-    assert _offset_provider is not None
-    connectivity = _offset_provider[offset_str]
+    offset_provider = offset_provider_cvar.get()
+    assert offset_provider is not None
+    connectivity = offset_provider[offset_str]
     assert isinstance(connectivity, common.Connectivity)
     return _List(
         shifted.deref()
@@ -1141,8 +1149,9 @@ class SparseListIterator:
     offsets: Sequence[OffsetPart] = dataclasses.field(default_factory=list, kw_only=True)
 
     def deref(self) -> Any:
-        assert _offset_provider is not None
-        connectivity = _offset_provider[self.list_offset]
+        offset_provider = offset_provider_cvar.get()
+        assert offset_provider is not None
+        connectivity = offset_provider[self.list_offset]
         assert isinstance(connectivity, common.Connectivity)
         return _List(
             shifted.deref()
@@ -1290,13 +1299,14 @@ def _column_dtype(elem: Any) -> np.dtype:
 @builtins.scan.register(EMBEDDED)
 def scan(scan_pass, is_forward: bool, init):
     def impl(*iters: ItIterator):
-        if _column_range is None:
+        column_range = column_range_cvar.get()
+        if column_range is None:
             raise RuntimeError("Column range is not defined, cannot scan.")
 
-        column_range = _column_range if is_forward else reversed(_column_range)
+        sorted_column_range = column_range if is_forward else reversed(column_range)
         state = init
-        col = Column(_column_range.start, np.zeros(len(_column_range), dtype=_column_dtype(init)))
-        for i in column_range:
+        col = Column(column_range.start, np.zeros(len(column_range), dtype=_column_dtype(init)))
+        for i in sorted_column_range:
             state = scan_pass(state, *map(shifted_scan_arg(i), iters))
             col[i] = state
 
@@ -1321,8 +1331,7 @@ def fendef_embedded(fun: Callable[..., None], *args: Any, **kwargs: Any):
     if "offset_provider" not in kwargs:
         raise RuntimeError("offset_provider not provided")
 
-    global _offset_provider
-    _offset_provider = kwargs["offset_provider"]
+    offset_provider = kwargs["offset_provider"]
 
     @runtime.closure.register(EMBEDDED)
     def closure(
@@ -1336,42 +1345,48 @@ def fendef_embedded(fun: Callable[..., None], *args: Any, **kwargs: Any):
         if not (is_located_field(out) or can_be_tuple_field(out)):
             raise TypeError("Out needs to be a located field.")
 
-        global _column_range
+        column_range = None
         column: Optional[ColumnDescriptor] = None
         if kwargs.get("column_axis") and kwargs["column_axis"].value in domain:
             column_axis = kwargs["column_axis"]
             column = ColumnDescriptor(column_axis.value, domain[column_axis.value])
             del domain[column_axis.value]
 
-            _column_range = column.col_range
+            column_range = column.col_range
 
         out = as_tuple_field(out) if can_be_tuple_field(out) else out
 
-        for pos in _domain_iterator(domain):
-            promoted_ins = [promote_scalars(inp) for inp in ins]
-            ins_iters = list(
-                make_in_iterator(
-                    inp,
-                    pos,
-                    column_axis=column.axis if column else None,
+        def _closure_runner():
+            # Set context variables before executing the closure
+            column_range_cvar.set(column_range)
+            offset_provider_cvar.set(offset_provider)
+
+            for pos in _domain_iterator(domain):
+                promoted_ins = [promote_scalars(inp) for inp in ins]
+                ins_iters = list(
+                    make_in_iterator(
+                        inp,
+                        pos,
+                        column_axis=column.axis if column else None,
+                    )
+                    for inp in promoted_ins
                 )
-                for inp in promoted_ins
-            )
-            res = sten(*ins_iters)
+                res = sten(*ins_iters)
 
-            if column is None:
-                assert _is_concrete_position(pos)
-                ordered_indices = get_ordered_indices(out.axes, pos)
-                out.field_setitem(ordered_indices, res)
-            else:
-                col_pos = pos.copy()
-                for k in column.col_range:
-                    col_pos[column.axis] = k
-                    assert _is_concrete_position(col_pos)
-                    ordered_indices = get_ordered_indices(out.axes, col_pos)
-                    out.field_setitem(ordered_indices, res[k])
-
-        _column_range = None
+                if column is None:
+                    assert _is_concrete_position(pos)
+                    ordered_indices = get_ordered_indices(out.axes, pos)
+                    out.field_setitem(ordered_indices, res)
+                else:
+                    col_pos = pos.copy()
+                    for k in column.col_range:
+                        col_pos[column.axis] = k
+                        assert _is_concrete_position(col_pos)
+                        ordered_indices = get_ordered_indices(out.axes, col_pos)
+                        out.field_setitem(ordered_indices, res[k])
+
+        ctx = cvars.copy_context()
+        ctx.run(_closure_runner)
 
     fun(*args)
 


Test patch:
diff --git a/tests/next_tests/integration_tests/feature_tests/cases.py b/tests/next_tests/integration_tests/feature_tests/cases.py
new file mode 100644
index 00000000..61afc376
--- /dev/null
+++ b/tests/next_tests/integration_tests/feature_tests/cases.py
@@ -0,0 +1,564 @@
+# GT4Py - GridTools Framework
+#
+# Copyright (c) 2014-2023, ETH Zurich
+# All rights reserved.
+#
+# This file is part of the GT4Py project and the GridTools framework.
+# GT4Py is free software: you can redistribute it and/or modify it under
+# the terms of the GNU General Public License as published by the
+# Free Software Foundation, either version 3 of the License, or any later
+# version. See the LICENSE.txt file at the top-level directory of this
+# distribution for a copy of the license or check <https://www.gnu.org/licenses/>.
+#
+# SPDX-License-Identifier: GPL-3.0-or-later
+
+from __future__ import annotations
+
+import dataclasses
+import functools
+import inspect
+import types
+import typing
+from typing import Any, Callable, Literal, Optional, Protocol, TypeAlias
+
+import numpy as np
+import pytest
+
+from gt4py.eve import extended_typing as xtyping
+from gt4py.eve.extended_typing import Self
+from gt4py.next import common
+from gt4py.next.ffront import decorator, fbuiltins
+from gt4py.next.iterator import embedded, ir as itir
+from gt4py.next.program_processors import processor_interface as ppi
+from gt4py.next.type_system import type_specifications as ts, type_translation
+
+from next_tests.integration_tests.feature_tests.ffront_tests.ffront_test_utils import (  # noqa: F401 #  fixture and aliases
+    Edge,
+    IDim,
+    Ioff,
+    JDim,
+    Joff,
+    KDim,
+    Koff,
+    Vertex,
+    fieldview_backend,
+    reduction_setup,
+)
+
+
+# mypy does not accept [IDim, ...] as a type
+IField: TypeAlias = common.Field[[IDim], np.int64]  # type: ignore [valid-type]
+IJKField: TypeAlias = common.Field[[IDim, JDim, KDim], np.int64]  # type: ignore [valid-type]
+IJKFloatField: TypeAlias = common.Field[[IDim, JDim, KDim], np.float64]  # type: ignore [valid-type]
+VField: TypeAlias = common.Field[[Vertex], np.int64]  # type: ignore [valid-type]
+EField: TypeAlias = common.Field[[Edge], np.int64]  # type: ignore [valid-type]
+
+# TODO(ricoh): unify the following with the `ffront_test_utils.reduction_setup`
+#   fixture if `ffront_test_utils.reduction_setup` is not completely superseded
+#   by `unstructured_case`.
+V2EDim = common.Dimension("V2E", kind=common.DimensionKind.LOCAL)
+E2VDim = common.Dimension("E2V", kind=common.DimensionKind.LOCAL)
+V2E = fbuiltins.FieldOffset("V2E", source=Edge, target=(Vertex, V2EDim))
+E2V = fbuiltins.FieldOffset("E2V", source=Vertex, target=(Edge, E2VDim))
+
+ScalarValue: TypeAlias = np.int32 | np.int64 | np.float32 | np.float64 | np.generic
+FieldValue: TypeAlias = common.Field | embedded.LocatedFieldImpl
+FieldViewArg: TypeAlias = FieldValue | ScalarValue | tuple["FieldViewArg", ...]
+FieldViewInout: TypeAlias = FieldValue | tuple["FieldViewInout", ...]
+ReferenceValue: TypeAlias = (
+    common.Field | np.typing.NDArray[ScalarValue] | tuple["ReferenceValue", ...]
+)
+OffsetProvider: TypeAlias = dict[str, common.Connectivity | common.Dimension]
+
+
+#: To allocate the return value of a field operator, we must pass
+#: something that is not an argument name. Currently this is  the
+#: literal string "return" (because it is read from annotations).
+#: This could change if implemented differently. RETURN acts as a
+#: proxy to avoid propagating the change to client code.
+RETURN = "return"
+
+
+class DataInitializer(Protocol):
+    @property
+    def scalar_value(self) -> ScalarValue:
+        ...
+
+    def scalar(self, dtype: np.typing.DTypeLike) -> ScalarValue:
+        # some unlikely numpy dtypes are picky about arguments
+        return np.dtype(dtype).type(self.scalar_value)  # type: ignore [arg-type]
+
+    def field(
+        self,
+        backend: ppi.ProgramProcessor,
+        sizes: dict[common.Dimension, int],
+        dtype: np.typing.DTypeLike,
+    ) -> FieldValue:
+        ...
+
+    def from_case(
+        self: Self,
+        case: Case,
+        fieldview_prog: decorator.FieldOperator | decorator.Program,
+        arg_name: str,
+    ) -> Self:
+        return self
+
+
+@dataclasses.dataclass
+class ConstInitializer(DataInitializer):
+    """Initialize with a given value across the coordinate space."""
+
+    value: ScalarValue
+
+    @property
+    def scalar_value(self) -> ScalarValue:
+        return self.value
+
+    def field(
+        self,
+        backend: ppi.ProgramProcessor,
+        sizes: dict[common.Dimension, int],
+        dtype: np.typing.DTypeLike,
+    ) -> FieldValue:
+        return embedded.np_as_located_field(*sizes.keys())(
+            np.full(tuple(sizes.values()), self.value, dtype=dtype)
+        )
+
+
+@dataclasses.dataclass(init=False)
+class ZeroInitializer(ConstInitializer):
+    """Initialize with zeros."""
+
+    def __init__(self):
+        self.value = 0
+
+
+@dataclasses.dataclass
+class UniqueInitializer(DataInitializer):
+    """
+    Initialize with a unique value in each coordinate point.
+
+    Data initialized with the same instance will also have unique values across
+    data containers.
+    """
+
+    start: int = 0
+
+    @property
+    def scalar_value(self) -> ScalarValue:
+        start = self.start
+        self.start += 1
+        return np.int64(start)
+
+    def field(
+        self,
+        backend: ppi.ProgramProcessor,
+        sizes: dict[common.Dimension, int],
+        dtype: np.typing.DTypeLike,
+    ) -> FieldValue:
+        start = self.start
+        svals = tuple(sizes.values())
+        n_data = int(np.prod(svals))
+        self.start += n_data
+        return embedded.np_as_located_field(*sizes.keys())(
+            np.arange(start, start + n_data, dtype=dtype).reshape(svals)
+        )
+
+    def from_case(
+        self: Self,
+        case: Case,
+        fieldview_prog: decorator.FieldOperator | decorator.Program,
+        arg_name: str,
+    ) -> Self:
+        param_types = get_param_types(fieldview_prog)
+        param_sizes = [
+            get_param_size(param_type, sizes=case.default_sizes)
+            for param_type in param_types.values()
+        ]
+        param_index = list(param_types.keys()).index(arg_name)
+        return self.__class__(start=self.start + sum(param_sizes[:param_index]))
+
+
+@dataclasses.dataclass(frozen=True)
+class Builder:
+    partial: functools.partial
+
+    def __call__(self, *args: Any, **kwargs: Any) -> Any:
+        return self.partial(*args, **kwargs)
+
+    def __getattr__(self, name: str) -> Any:
+        raise AttributeError(f"No setter for argument {name}.")
+
+
+@typing.overload
+def make_builder(*args: Callable) -> Callable[..., Builder]:
+    ...
+
+
+@typing.overload
+def make_builder(
+    *args: Literal[None], **kwargs: dict[str, Any]
+) -> Callable[[Callable], Callable[..., Builder]]:
+    ...
+
+
+@typing.overload
+def make_builder(
+    *args: Optional[Callable], **kwargs: dict[str, Any]
+) -> Callable[[Callable], Callable[..., Builder]] | Callable[..., Builder]:
+    ...
+
+
+# TODO(ricoh): Think about improving the type hints using `typing.ParamSpec`.
+def make_builder(
+    *args: Optional[Callable], **kwargs: dict[str, Any]
+) -> Callable[[Callable], Callable[..., Builder]] | Callable[..., Builder]:
+    """Create a fluid interface for a function with many arguments."""
+
+    def make_builder_inner(func: Callable) -> Callable[..., Builder]:
+        def make_setter(argname: str) -> Callable[[Builder, Any], Builder]:
+            def setter(self: Builder, arg: Any) -> Builder:
+                return self.__class__(
+                    partial=functools.partial(
+                        self.partial.func,
+                        *self.partial.args,
+                        **(self.partial.keywords | {argname: arg}),
+                    )
+                )
+
+            setter.__name__ = argname
+            return setter
+
+        def make_flag_setter(
+            flag_name: str, flag_kwargs: dict[str, Any]
+        ) -> Callable[[Builder], Builder]:
+            def setter(self: Builder) -> Builder:
+                return self.__class__(
+                    partial=functools.partial(
+                        self.partial.func,
+                        *self.partial.args,
+                        **(self.partial.keywords | flag_kwargs),
+                    )
+                )
+
+            setter.__name__ = flag_name
+            return setter
+
+        argspec = inspect.getfullargspec(func)
+
+        @dataclasses.dataclass(frozen=True)
+        class NewBuilder(Builder):
+            ...
+
+        for argname in argspec.args + argspec.kwonlyargs:
+            setattr(NewBuilder, argname, make_setter(argname))
+
+        for flag, flag_kwargs in kwargs.items():
+            setattr(NewBuilder, flag, make_flag_setter(flag, flag_kwargs))
+
+        func_snake_words = func.__name__.split("_")
+        func_camel_name = "".join(word.capitalize() for word in func_snake_words)
+
+        NewBuilder.__name__ = f"{func_camel_name}Builder"
+
+        return lambda *args, **kwargs: NewBuilder(functools.partial(func, *args, **kwargs))
+
+    if 0 < len(args) <= 1 and args[0] is not None:
+        return make_builder_inner(args[0])
+    if len(args) > 1:
+        raise ValueError(f"make_builder takes only one positional argument, {len(args)} received!")
+    return make_builder_inner
+
+
+@make_builder(zeros={"strategy": ZeroInitializer()}, unique={"strategy": UniqueInitializer()})
+def allocate(
+    case: Case,
+    fieldview_prog: decorator.FieldOperator | decorator.Program,
+    name: str,
+    *,
+    sizes: Optional[dict[common.Dimension, int]] = None,
+    strategy: Optional[DataInitializer] = None,
+    dtype: Optional[np.typing.DTypeLike] = None,
+    extend: Optional[dict[common.Dimension, tuple[int, int]]] = None,
+) -> FieldViewArg:
+    """
+    Allocate a parameter or return value from a fieldview code object.
+
+    Args:
+        case: The test case.
+        fieldview_prog: The field operator or program to be verified.
+        name: The name of the input argument to allocate, or ``RETURN``
+            for the return value of a field operator.
+    Keyword Args:
+        sizes: Override for the test case dimension sizes.
+            Use with caution.
+        strategy: How to initialize the data.
+        dtype: Override for the dtype in the argument's type hint.
+        extend: Lower and upper size extension per dimension.
+            Useful for shifted fields, which must start off bigger
+            than the output field in the shifted dimension.
+    """
+    sizes = extend_sizes(case.default_sizes | (sizes or {}), extend)
+    arg_type = get_param_types(fieldview_prog)[name]
+    if strategy is None:
+        if name in ["out", RETURN]:
+            strategy = ZeroInitializer()
+        else:
+            strategy = UniqueInitializer()
+    return _allocate_from_type(
+        case=case,
+        arg_type=arg_type,
+        sizes=sizes,
+        dtype=dtype,
+        strategy=strategy.from_case(case=case, fieldview_prog=fieldview_prog, arg_name=name),
+    )
+
+
+def run(
+    case: Case,
+    fieldview_prog: decorator.FieldOperator | decorator.Program,
+    *args: FieldViewArg,
+    **kwargs: Any,
+) -> None:
+    """Run fieldview code in the context of a given test case."""
+    if kwargs.get("offset_provider", None) is None:
+        kwargs["offset_provider"] = case.offset_provider
+    fieldview_prog.with_grid_type(case.grid_type).with_backend(case.backend)(*args, **kwargs)
+
+
+def verify(
+    case: Case,
+    fieldview_prog: decorator.FieldOperator | decorator.Program,
+    *args: FieldViewArg,
+    ref: ReferenceValue,
+    out: Optional[FieldViewInout] = None,
+    inout: Optional[FieldViewInout] = None,
+    offset_provider: Optional[OffsetProvider] = None,
+    comparison: Callable[[Any, Any], bool] = np.allclose,
+) -> None:
+    """
+    Check the result of executing a fieldview program or operator against ref.
+
+    Args:
+        case: The test case.
+        fieldview_prog: The field operator or program to be verified.
+        *args: positional input arguments to the fieldview code.
+    Keyword Args:
+        ref: A field or array which will be compared to the results of the
+            fieldview code.
+        out: If given will be passed to the fieldview code as ``out=`` keyword
+            argument. This will hold the results and be used to compare
+            to ``ref``.
+        inout: If ``out`` is not passed it is assumed that the fieldview code
+            does not take an ``out`` keyword argument, so some of the other
+            arguments must be in/out parameters. Pass the according field
+            or tuple of fields here and they will be compared to ``ref`` under
+            the assumption that the fieldview code stores it's results in
+            them.
+        offset_provied: An override for the test case's offset_provider.
+            Use with care!
+        comparison: A comparison function, which will be called as
+            ``comparison(ref, <out | inout>)`` and should return a boolean.
+
+    One of ``out`` or ``inout`` must be passed. If ``out`` is passed it will be
+    used as an argument to the fieldview program and compared against ``ref``.
+    Else, ``inout`` will not be passed and compared to ``ref``.
+    """
+    if out:
+        run(
+            case,
+            fieldview_prog,
+            *args,
+            out=out,
+            offset_provider=offset_provider,
+        )
+    else:
+        run(case, fieldview_prog, *args, offset_provider=offset_provider)
+
+    out_comp = out or inout
+    out_comp_str = str(out_comp)
+    assert out_comp is not None
+    if hasattr(out_comp, "array"):
+        out_comp_str = str(out_comp.array())
+    assert comparison(ref, out_comp), (
+        f"Verification failed:\n"
+        f"\tcomparison={comparison.__name__}(ref, out)\n"
+        f"\tref = {ref}\n\tout = {out_comp_str}"
+    )
+
+
+def verify_with_default_data(
+    case: Case,
+    fieldop: decorator.FieldOperator,
+    ref: Callable,
+    comparison: Callable[[Any, Any], bool] = np.allclose,
+) -> None:
+    """
+    Check the fieldview code against a reference calculation.
+
+    This is a convenience function to reduce boiler plate
+    and is not meant to hide logic for complex cases. The
+    ``verify`` function allows more fine grained control for such tests.
+
+    Args:
+        case: The test case.
+        fieldview_prog: The field operator or program to be verified.
+        ref: A callable which will be called with all the input arguments
+            of the fieldview code, after applying ``.array()`` on the fields.
+        comparison: A comparison function, which will be called as
+            ``comparison(ref, <out | inout>)`` and should return a boolean.
+    """
+    inps, kwfields = get_default_data(case, fieldop)
+    ref_args = tuple(i.array() if hasattr(i, "array") else i for i in inps)
+    verify(
+        case,
+        fieldop,
+        *inps,
+        **kwfields,
+        ref=ref(*ref_args),
+        offset_provider=case.offset_provider,
+        comparison=comparison,
+    )
+
+
+@pytest.fixture
+def no_default_backend():
+    """Temporarily switch off default backend for feature tests."""
+    backup_backend = decorator.DEFAULT_BACKEND
+    decorator.DEFAULT_BACKEND = no_backend
+    yield
+    decorator.DEFAULT_BACKEND = backup_backend
+
+
+@pytest.fixture
+def cartesian_case(no_default_backend, fieldview_backend):  # noqa: F811 # fixtures
+    yield Case(
+        fieldview_backend,
+        offset_provider={"Ioff": IDim, "Joff": JDim, "Koff": KDim},
+        default_sizes={IDim: 10, JDim: 10, KDim: 10},
+        grid_type=common.GridType.CARTESIAN,
+    )
+
+
+@pytest.fixture
+def unstructured_case(
+    no_default_backend, reduction_setup, fieldview_backend  # noqa: F811 # fixtures
+):
+    yield Case(
+        fieldview_backend,
+        offset_provider=reduction_setup.offset_provider,
+        default_sizes={
+            Vertex: reduction_setup.num_vertices,
+            Edge: reduction_setup.num_edges,
+        },
+        grid_type=common.GridType.UNSTRUCTURED,
+    )
+
+
+def _allocate_from_type(
+    case: Case,
+    arg_type: ts.TypeSpec,
+    sizes: dict[common.Dimension, int],
+    strategy: DataInitializer,
+    dtype: Optional[np.typing.DTypeLike] = None,
+    tuple_start: Optional[int] = None,
+) -> FieldViewArg:
+    """Allocate data based on the type or a (nested) tuple thereof."""
+    match arg_type:
+        case ts.FieldType(dims=dims, dtype=arg_dtype):
+            return strategy.field(
+                backend=case.backend,
+                sizes={dim: sizes[dim] for dim in dims},
+                dtype=dtype or arg_dtype.kind.name.lower(),
+            )
+        case ts.ScalarType(kind=kind):
+            return strategy.scalar(dtype=dtype or kind.name.lower())
+        case ts.TupleType(types=types):
+            return tuple(
+                (
+                    _allocate_from_type(
+                        case=case, arg_type=t, sizes=sizes, dtype=dtype, strategy=strategy
+                    )
+                    for t in types
+                )
+            )
+        case _:
+            raise TypeError(
+                f"Can not allocate for type {arg_type} with initializer {strategy or 'default'}"
+            )
+
+
+def get_param_types(
+    fieldview_prog: decorator.FieldOperator | decorator.Program,
+) -> dict[str, ts.TypeSpec]:
+    if fieldview_prog.definition is None:
+        raise ValueError(
+            f"test cases do not support {type(fieldview_prog)} with empty .definition attribute (as you would get from .as_program())!"
+        )
+    annotations = xtyping.get_type_hints(fieldview_prog.definition)
+    return {
+        name: type_translation.from_type_hint(type_hint) for name, type_hint in annotations.items()
+    }
+
+
+def get_param_size(param_type: ts.TypeSpec, sizes: dict[common.Dimension, int]) -> int:
+    match param_type:
+        case ts.FieldType(dims=dims):
+            return int(np.prod([sizes[dim] for dim in sizes if dim in dims]))
+        case ts.ScalarType(shape=shape):
+            return int(np.prod(shape)) if shape else 1
+        case ts.TupleType(types):
+            return sum([get_param_size(t, sizes=sizes) for t in types])
+        case _:
+            raise TypeError(f"Can not get size for parameter of type {param_type}")
+
+
+def extend_sizes(
+    sizes: dict[common.Dimension, int],
+    extend: Optional[dict[common.Dimension, tuple[int, int]]] = None,
+) -> dict[common.Dimension, int]:
+    """Calculate the sizes per dimension given a set of extensions."""
+    sizes = sizes.copy()
+    if extend:
+        for dim, (lower, upper) in extend.items():
+            sizes[dim] += upper - lower
+    return sizes
+
+
+def get_default_data(
+    case: Case,
+    fieldview_prog: decorator.FieldOperator | decorator.Program,
+) -> tuple[tuple[common.Field | ScalarValue | tuple, ...], dict[str, common.Field]]:
+    """
+    Allocate default data for a fieldview code object given a test case.
+
+    Meant to reduce boiler plate for simple cases, everything else
+    should rely on ``allocate()``.
+    """
+    param_types = get_param_types(fieldview_prog)
+    kwfields: dict[str, Any] = {}
+    if RETURN in param_types:
+        if not isinstance(param_types[RETURN], types.NoneType):
+            kwfields = {"out": allocate(case, fieldview_prog, RETURN).zeros()()}
+        param_types.pop(RETURN)
+    if "out" in param_types:
+        kwfields = {"out": allocate(case, fieldview_prog, "out").zeros()()}
+        param_types.pop("out")
+    inps = tuple(allocate(case, fieldview_prog, name)() for name in param_types)
+    return inps, kwfields
+
+
+def no_backend(program: itir.FencilDefinition, *args: Any, **kwargs: Any) -> None:
+    """Temporary default backend to not accidentally test the wrong backend."""
+    raise ValueError("No backend selected! Backend selection is mandatory in tests.")
+
+
+@dataclasses.dataclass
+class Case:
+    """Parametrizable components for single feature integration tests."""
+
+    backend: ppi.ProgramProcessor
+    offset_provider: dict[str, common.Connectivity | common.Dimension]
+    default_sizes: dict[common.Dimension, int]
+    grid_type: common.GridType
diff --git a/tests/next_tests/integration_tests/feature_tests/ffront_tests/test_execution.py b/tests/next_tests/integration_tests/feature_tests/ffront_tests/test_execution.py
index 6fb9d880..821953cd 100644
--- a/tests/next_tests/integration_tests/feature_tests/ffront_tests/test_execution.py
+++ b/tests/next_tests/integration_tests/feature_tests/ffront_tests/test_execution.py
@@ -38,75 +38,72 @@ from gt4py.next.ffront.foast_passes.type_deduction import FieldOperatorTypeDeduc
 from gt4py.next.iterator.embedded import index_field, np_as_located_field
 from gt4py.next.program_processors.runners import gtfn_cpu
 
-from next_tests.integration_tests.feature_tests.ffront_tests.ffront_test_utils import *
+from next_tests.integration_tests.feature_tests import cases
+from next_tests.integration_tests.feature_tests.cases import (
+    E2V,
+    V2E,
+    E2VDim,
+    Edge,
+    IDim,
+    Ioff,
+    JDim,
+    Joff,
+    KDim,
+    Koff,
+    V2EDim,
+    Vertex,
+    cartesian_case,
+    no_default_backend,
+    unstructured_case,
+)
+from next_tests.integration_tests.feature_tests.ffront_tests.ffront_test_utils import (
+    Cell,
+    fieldview_backend,
+    reduction_setup,
+    size,
+)
 
 
-def test_copy(fieldview_backend):
-    a_I_float = np_as_located_field(IDim)(np.random.randn(size).astype("float64"))
-    b_I_float = np_as_located_field(IDim)(np.random.randn(size).astype("float64"))
-
-    @field_operator(backend=fieldview_backend)
-    def copy(inp: Field[[IDim], float64]) -> Field[[IDim], float64]:
-        field_tuple = (inp, inp)
+def test_copy(cartesian_case):  # noqa: F811 # fixtures
+    @field_operator
+    def testee(a: cases.IJKField) -> cases.IJKField:
+        field_tuple = (a, a)
         field_0 = field_tuple[0]
         field_1 = field_tuple[1]
         return field_0
 
-    copy(a_I_float, out=b_I_float, offset_provider={})
-
-    assert np.allclose(a_I_float, b_I_float)
-
-
-def test_multicopy(fieldview_backend):
-    inp0 = np_as_located_field(IDim)(np.random.randn(size).astype("float64"))
-    inp1 = np_as_located_field(IDim)(np.random.randn(size).astype("float32"))
-    out0 = np_as_located_field(IDim)(np.zeros((size), dtype=float64))
-    out1 = np_as_located_field(IDim)(np.zeros((size), dtype=float32))
-
-    @field_operator(backend=fieldview_backend)
-    def multicopy(
-        inp1: Field[[IDim], float64], inp2: Field[[IDim], float32]
-    ) -> tuple[Field[[IDim], float64], Field[[IDim], float32]]:
-        return inp1, inp2
-
-    multicopy(inp0, inp1, out=(out0, out1), offset_provider={})
+    cases.verify_with_default_data(cartesian_case, testee, ref=lambda a: a)
 
-    assert np.allclose(inp0, out0)
-    assert np.allclose(inp1, out1)
 
+def test_multicopy(cartesian_case):  # noqa: F811 # fixtures
+    @field_operator
+    def testee(a: cases.IJKField, b: cases.IJKField) -> tuple[cases.IJKField, cases.IJKField]:
+        return a, b
 
-def test_cartesian_shift(fieldview_backend):
-    a = np_as_located_field(IDim)(np.arange(size + 1, dtype=np.float64))
-    out_I_float = np_as_located_field(IDim)(np.zeros((size,), dtype=float64))
-
-    @field_operator(backend=fieldview_backend)
-    def shift_by_one(inp: Field[[IDim], float64]) -> Field[[IDim], float64]:
-        return inp(Ioff[1])
-
-    @program(backend=fieldview_backend)
-    def fencil(inp: Field[[IDim], float64], out: Field[[IDim], float64]) -> None:
-        shift_by_one(inp, out=out)
-
-    fencil(a, out_I_float, offset_provider={"Ioff": IDim})
-
-    assert np.allclose(out_I_float.array(), np.arange(1, 11))
+    cases.verify_with_default_data(cartesian_case, testee, ref=lambda a, b: (a, b))
 
 
-def test_unstructured_shift(reduction_setup, fieldview_backend):
-    E2V = reduction_setup.E2V
+def test_cartesian_shift(cartesian_case):  # noqa: F811 # fixtures
+    @field_operator
+    def testee(a: cases.IJKField) -> cases.IJKField:
+        return a(Ioff[1])
 
-    a = np_as_located_field(Vertex)(np.arange(0, reduction_setup.num_vertices, dtype=np.float64))
-    b = np_as_located_field(Edge)(np.zeros(reduction_setup.num_edges))
+    a = cases.allocate(cartesian_case, testee, "a").extend({IDim: (0, 1)})()
+    out = cases.allocate(cartesian_case, testee, cases.RETURN)()
 
-    @field_operator(backend=fieldview_backend)
-    def shift_unstructured(inp: Field[[Vertex], float64]) -> Field[[Edge], float64]:
-        return inp(E2V[0])
+    cases.verify(cartesian_case, testee, a, out=out, ref=a[1:])
 
-    shift_unstructured(a, out=b, offset_provider={"E2V": reduction_setup.offset_provider["E2V"]})
 
-    ref = np.asarray(a)[reduction_setup.offset_provider["E2V"].table[slice(0, None), 0]]
+def test_unstructured_shift(unstructured_case):  # noqa: F811 # fixtures
+    @field_operator
+    def testee(a: cases.VField) -> cases.EField:
+        return a(E2V[0])
 
-    assert np.allclose(b, ref)
+    cases.verify_with_default_data(
+        unstructured_case,
+        testee,
+        ref=lambda a: a[unstructured_case.offset_provider["E2V"].table[:, 0]],
+    )
 
 
 def test_composed_unstructured_shift(reduction_setup, fieldview_backend):
@@ -149,223 +146,187 @@ def test_composed_unstructured_shift(reduction_setup, fieldview_backend):
         assert np.allclose(b, ref)
 
 
-def test_fold_shifts(fieldview_backend):
+def test_fold_shifts(cartesian_case):  # noqa: F811 # fixtures
     """Shifting the result of an addition should work."""
-    a = np_as_located_field(IDim)(np.arange(size + 1, dtype=np.float64))
-    b = np_as_located_field(IDim)(np.ones((size + 2)) * 2)
-    out_I_float = np_as_located_field(IDim)(np.zeros((size,), dtype=float64))
 
-    @field_operator(backend=fieldview_backend)
-    def auto_lift(
-        inp1: Field[[IDim], float64], inp2: Field[[IDim], float64]
-    ) -> Field[[IDim], float64]:
-        tmp = inp1 + inp2(Ioff[1])
+    @field_operator
+    def testee(a: cases.IJKField, b: cases.IJKField) -> cases.IJKField:
+        tmp = a + b(Ioff[1])
         return tmp(Ioff[1])
 
-    @program(backend=fieldview_backend)
-    def fencil(
-        inp1: Field[[IDim], float64], inp2: Field[[IDim], float64], out: Field[[IDim], float64]
-    ) -> None:
-        auto_lift(inp1, inp2, out=out)
-
-    fencil(a, b, out_I_float, offset_provider={"Ioff": IDim})
-
-    assert np.allclose(a[1:] + b[2:], out_I_float)
+    a = cases.allocate(cartesian_case, testee, "a").extend({cases.IDim: (0, 1)})()
+    b = cases.allocate(cartesian_case, testee, "b").extend({cases.IDim: (0, 2)})()
+    out = cases.allocate(cartesian_case, testee, cases.RETURN)()
 
+    cases.verify(cartesian_case, testee, a, b, out=out, ref=a[1:] + b[2:])
 
-def test_tuples(fieldview_backend):
-    a_I_float = np_as_located_field(IDim)(np.random.randn(size).astype("float64"))
-    b_I_float = np_as_located_field(IDim)(np.random.randn(size).astype("float64"))
-    out_I_float = np_as_located_field(IDim)(np.zeros((size,), dtype=float64))
 
-    @field_operator(backend=fieldview_backend)
-    def tuples(
-        inp1: Field[[IDim], float64], inp2: Field[[IDim], float64]
-    ) -> Field[[IDim], float64]:
-        inps = inp1, inp2
+def test_tuples(cartesian_case):  # noqa: F811 # fixtures
+    @field_operator
+    def testee(a: cases.IJKFloatField, b: cases.IJKFloatField) -> cases.IJKFloatField:
+        inps = a, b
         scalars = 1.3, float64(5.0), float64("3.4")
         return (inps[0] * scalars[0] + inps[1] * scalars[1]) * scalars[2]
 
-    @program(backend=fieldview_backend)
-    def fencil(
-        inp1: Field[[IDim], float64], inp2: Field[[IDim], float64], out: Field[[IDim], float64]
-    ) -> None:
-        tuples(inp1, inp2, out=out)
-
-    fencil(a_I_float, b_I_float, out_I_float, offset_provider={})
-
-    assert np.allclose((a_I_float.array() * 1.3 + b_I_float.array() * 5.0) * 3.4, out_I_float)
+    cases.verify_with_default_data(
+        cartesian_case, testee, ref=lambda a, b: (a * 1.3 + b * 5.0) * 3.4
+    )
 
 
-def test_scalar_arg(fieldview_backend):
+def test_scalar_arg(unstructured_case):  # noqa: F811 # fixtures
     """Test scalar argument being turned into 0-dim field."""
-    inp = 5.0
-    out = np_as_located_field(Vertex)(np.zeros([size]))
-
-    @field_operator(backend=fieldview_backend)
-    def scalar_arg(scalar_inp: float64) -> Field[[Vertex], float64]:
-        return broadcast(scalar_inp + 1.0, (Vertex,))
-
-    scalar_arg(inp, out=out, offset_provider={})
-
-    ref = np.full([size], 6.0)
-    assert np.allclose(ref, out.array())
-
-
-def test_nested_scalar_arg(fieldview_backend):
-    inp = 5.0
-    out = np_as_located_field(Vertex)(np.zeros([size]))
 
-    @field_operator(backend=fieldview_backend)
-    def scalar_arg_inner(scalar_inp: float64) -> Field[[Vertex], float64]:
-        return broadcast(scalar_inp + 1.0, (Vertex,))
-
-    @field_operator(backend=fieldview_backend)
-    def scalar_arg(scalar_inp: float64) -> Field[[Vertex], float64]:
-        return scalar_arg_inner(scalar_inp + 1.0)
-
-    scalar_arg(inp, out=out, offset_provider={})
+    @field_operator
+    def testee(a: int) -> cases.VField:
+        return broadcast(a + 1, (Vertex,))
+
+    cases.verify_with_default_data(
+        unstructured_case,
+        testee,
+        ref=lambda a: np.full(
+            [unstructured_case.default_sizes[Vertex]],
+            a + 1,
+            dtype=int,
+        ),
+        comparison=lambda a, b: np.all(a == b),
+    )
 
-    ref = np.full([size], 7.0)
-    assert np.allclose(ref, out.array())
 
+def test_nested_scalar_arg(unstructured_case):  # noqa: F811 # fixtures
+    @field_operator
+    def testee_inner(a: int) -> cases.VField:
+        return broadcast(a + 1, (Vertex,))
 
-def test_scalar_arg_with_field(fieldview_backend):
-    if fieldview_backend in [gtfn_cpu.run_gtfn, gtfn_cpu.run_gtfn_imperative]:
-        pytest.xfail("IndexFields and ConstantFields are not supported yet.")
+    @field_operator
+    def testee(a: int) -> cases.VField:
+        return testee_inner(a + 1)
 
-    inp = index_field(Edge, dtype=float64)
-    factor = 3.0
-    out = np_as_located_field(Edge)(np.zeros((size,), dtype=np.float64))
+    cases.verify_with_default_data(
+        unstructured_case,
+        testee,
+        ref=lambda a: np.full([unstructured_case.default_sizes[Vertex]], a + 2, dtype=int),
+    )
 
-    @field_operator(backend=fieldview_backend)
-    def scalar_and_field_args(
-        inp: Field[[Edge], float64], factor: float64
-    ) -> Field[[Edge], float64]:
-        tmp = factor * inp
-        return tmp(EdgeOffset[1])
 
-    @program(backend=fieldview_backend)
-    def fencil(out: Field[[Edge], float64], inp: Field[[Edge], float64], factor: float64) -> None:
-        scalar_and_field_args(inp, factor, out=out)
+def test_scalar_arg_with_field(cartesian_case):  # noqa: F811 # fixtures
+    @field_operator
+    def testee(a: cases.IJKField, b: int) -> cases.IJKField:
+        tmp = b * a
+        return tmp(Ioff[1])
 
-    fencil(out, inp, factor, offset_provider={"EdgeOffset": Edge})
+    a = cases.allocate(cartesian_case, testee, "a").extend({IDim: (0, 1)})()
+    b = cases.allocate(cartesian_case, testee, "b")()
+    out = cases.allocate(cartesian_case, testee, cases.RETURN)()
+    ref = a.array()[1:] * b
 
-    ref = np.arange(1, size + 1) * factor
-    assert np.allclose(ref, out.array())
+    cases.verify(cartesian_case, testee, a, b, out=out, ref=ref)
 
 
-def test_scalar_in_domain_spec_and_fo_call(fieldview_backend):
-    if fieldview_backend in [gtfn_cpu.run_gtfn, gtfn_cpu.run_gtfn_imperative]:
+def test_scalar_in_domain_spec_and_fo_call(cartesian_case):  # noqa: F811 # fixtures
+    if cartesian_case.backend in [gtfn_cpu.run_gtfn, gtfn_cpu.run_gtfn_imperative]:
         pytest.xfail(
             "Scalar arguments not supported to be used in both domain specification "
             "and as an argument to a field operator."
         )
 
-    size = 10
-    out = np_as_located_field(Vertex)(np.zeros(10, dtype=int))
-
-    @field_operator(backend=fieldview_backend)
-    def foo(size: int) -> Field[[Vertex], int]:
-        return broadcast(size, (Vertex,))
-
-    @program(backend=fieldview_backend)
-    def bar(size: int, out: Field[[Vertex], int]):
-        foo(size, out=out, domain={Vertex: (0, size)})
+    @field_operator
+    def testee_op(size: int) -> cases.IField:
+        return broadcast(size, (IDim,))
 
-    bar(size, out, offset_provider={})
+    @program
+    def testee(size: int, out: cases.IField):
+        testee_op(size, out=out, domain={IDim: (0, size)})
 
-    assert (out.array() == size).all()
+    size = cartesian_case.default_sizes[IDim]
+    out = cases.allocate(cartesian_case, testee, "out").zeros()()
 
+    cases.verify(
+        cartesian_case, testee, size, out=out, ref=np.full_like(out.array(), size, dtype=int)
+    )
 
-def test_scalar_scan(fieldview_backend):
-    size = 10
-    KDim = Dimension("K", kind=DimensionKind.VERTICAL)
-    qc = np_as_located_field(IDim, KDim)(np.zeros((size, size)))
-    scalar = 1.0
-    expected = np.full((size, size), np.arange(start=1, stop=11, step=1).astype(float64))
 
+def test_scalar_scan(cartesian_case):  # noqa: F811 # fixtures
     @scan_operator(axis=KDim, forward=True, init=(0.0))
-    def _scan_scalar(carry: float, qc_in: float, scalar: float) -> float:
-        qc = qc_in + carry + scalar
+    def testee_scan(state: float, qc_in: float, scalar: float) -> float:
+        qc = qc_in + state + scalar
         return qc
 
-    @program(backend=fieldview_backend)
-    def scan_scalar(qc: Field[[IDim, KDim], float], scalar: float):
-        _scan_scalar(qc, scalar, out=qc)
+    @program
+    def testee(qc: Field[[IDim, KDim], float], scalar: float):
+        testee_scan(qc, scalar, out=qc)
 
-    scan_scalar(qc, scalar, offset_provider={})
-    assert np.allclose(np.asarray(qc), expected)
+    qc = cases.allocate(cartesian_case, testee, "qc").zeros()()
+    scalar = 1.0
+    ksize = cartesian_case.default_sizes[KDim]
+    expected = np.full((ksize, ksize), np.arange(start=1, stop=11, step=1).astype(float64))
 
+    cases.verify(cartesian_case, testee, qc, scalar, inout=qc, ref=expected)
 
-def test_tuple_scalar_scan(fieldview_backend):
-    if fieldview_backend in [gtfn_cpu.run_gtfn, gtfn_cpu.run_gtfn_imperative]:
-        pytest.xfail("Scalar tuple arguments are not supported in gtfn yet.")
 
-    size = 10
-    KDim = Dimension("K", kind=DimensionKind.VERTICAL)
-    qc = np_as_located_field(IDim, KDim)(np.zeros((size, size)))
-    tuple_scalar = (1.0, (1.0, 0.0))
-    expected = np.full((size, size), np.arange(start=1, stop=11, step=1).astype(float64))
+def test_tuple_scalar_scan(cartesian_case):  # noqa: F811 # fixtures
+    if cartesian_case.backend in [gtfn_cpu.run_gtfn, gtfn_cpu.run_gtfn_imperative]:
+        pytest.xfail("Scalar tuple arguments are not supported in gtfn yet.")
 
     @scan_operator(axis=KDim, forward=True, init=0.0)
-    def _scan_tuple_scalar(
+    def testee_scan(
         state: float, qc_in: float, tuple_scalar: tuple[float, tuple[float, float]]
     ) -> float:
         return (qc_in + state + tuple_scalar[1][0] + tuple_scalar[1][1]) / tuple_scalar[0]
 
-    @field_operator(backend=fieldview_backend)
-    def scan_tuple_scalar(
+    @field_operator
+    def testee_op(
         qc: Field[[IDim, KDim], float], tuple_scalar: tuple[float, tuple[float, float]]
     ) -> Field[[IDim, KDim], float]:
-        return _scan_tuple_scalar(qc, tuple_scalar)
-
-    scan_tuple_scalar(qc, tuple_scalar, out=qc, offset_provider={})
-    assert np.allclose(np.asarray(qc), expected)
-
-
-def test_astype_int(fieldview_backend):
-    size = 10
-    b_float_64 = np_as_located_field(IDim)(np.full((size,), fill_value=1.5, dtype=np.float64))
-    c_int64 = np_as_located_field(IDim)(np.ones((size,), dtype=np.int64))
-    out_int_64 = np_as_located_field(IDim)(np.zeros((size,), dtype=np.int64))
-
-    @field_operator(backend=fieldview_backend)
-    def astype_fieldop_int(b: Field[[IDim], float64]) -> Field[[IDim], int64]:
-        d = astype(b, int64)
-        return d
-
-    astype_fieldop_int(b_float_64, out=out_int_64, offset_provider={})
-    assert np.allclose(c_int64.array(), out_int_64)
+        return testee_scan(qc, tuple_scalar)
 
+    qc = cases.allocate(cartesian_case, testee_op, "qc").zeros()()
+    tuple_scalar = (1.0, (1.0, 0.0))
+    ksize = cartesian_case.default_sizes[KDim]
+    expected = np.full((ksize, ksize), np.arange(start=1, stop=11, step=1).astype(float64))
 
-def test_astype_bool(fieldview_backend):
-    b_float_64 = np_as_located_field(IDim)(np.full((size,), fill_value=0.5, dtype=np.float64))
-    c_bool = np_as_located_field(IDim)(np.ones((size,), dtype=bool))
-    out_bool = np_as_located_field(IDim)(np.zeros((size,), dtype=bool))
+    cases.verify(cartesian_case, testee_op, qc, tuple_scalar, out=qc, ref=expected)
 
-    @field_operator(backend=fieldview_backend)
-    def astype_fieldop_bool(b: Field[[IDim], float64]) -> Field[[IDim], bool]:
-        d = astype(b, bool)
-        return d
 
-    astype_fieldop_bool(b_float_64, out=out_bool, offset_provider={})
-    assert np.allclose(c_bool, out_bool)
+def test_astype_int(cartesian_case):  # noqa: F811 # fixtures
+    @field_operator
+    def testee(a: Field[[IDim], float64]) -> Field[[IDim], int64]:
+        b = astype(a, int64)
+        return b
+
+    cases.verify_with_default_data(
+        cartesian_case,
+        testee,
+        ref=lambda a: a.astype(int),
+        comparison=lambda a, b: np.all(a == b),
+    )
 
 
-def test_astype_float(fieldview_backend):
-    c_float64 = np_as_located_field(IDim)(
-        np.full((size,), fill_value=np.float64("5e300"), dtype=np.float64)
+def test_astype_bool(cartesian_case):  # noqa: F811 # fixtures
+    @field_operator
+    def testee(a: Field[[IDim], float64]) -> Field[[IDim], bool]:
+        b = astype(a, bool)
+        return b
+
+    cases.verify_with_default_data(
+        cartesian_case,
+        testee,
+        ref=lambda a: a.astype(bool),
+        comparison=lambda a, b: np.all(a == b),
     )
-    out_int_32 = np_as_located_field(IDim)(np.zeros((size,), dtype=np.float32))
 
-    @field_operator(backend=fieldview_backend)
-    def astype_fieldop_float(b: Field[[IDim], float64]) -> Field[[IDim], np.float32]:
-        d = astype(b, float32)
-        return d
 
-    astype_fieldop_float(c_float64, out=out_int_32, offset_provider={})
-    assert np.all(out_int_32.array() == np.inf)
+def test_astype_float(cartesian_case):  # noqa: F811 # fixtures
+    @field_operator
+    def testee(a: Field[[IDim], float64]) -> Field[[IDim], np.float32]:
+        b = astype(a, float32)
+        return b
+
+    cases.verify_with_default_data(
+        cartesian_case,
+        testee,
+        ref=lambda a: a.astype(np.float32),
+        comparison=lambda a, b: np.all(a == b),
+    )
 
 
 def test_offset_field(fieldview_backend):
diff --git a/tests/next_tests/integration_tests/feature_tests/test_util_cases.py b/tests/next_tests/integration_tests/feature_tests/test_util_cases.py
new file mode 100644
index 00000000..3b5cf1a5
--- /dev/null
+++ b/tests/next_tests/integration_tests/feature_tests/test_util_cases.py
@@ -0,0 +1,103 @@
+# GT4Py - GridTools Framework
+#
+# Copyright (c) 2014-2023, ETH Zurich
+# All rights reserved.
+#
+# This file is part of the GT4Py project and the GridTools framework.
+# GT4Py is free software: you can redistribute it and/or modify it under
+# the terms of the GNU General Public License as published by the
+# Free Software Foundation, either version 3 of the License, or any later
+# version. See the LICENSE.txt file at the top-level directory of this
+# distribution for a copy of the license or check <https://www.gnu.org/licenses/>.
+#
+# SPDX-License-Identifier: GPL-3.0-or-later
+
+import numpy as np
+import pytest
+
+from gt4py.next import common
+from gt4py.next.ffront.decorator import field_operator
+from gt4py.next.program_processors.runners import roundtrip
+
+from next_tests.integration_tests.feature_tests import cases
+from next_tests.integration_tests.feature_tests.cases import (  # noqa: F401 # fixtures
+    cartesian_case,
+    fieldview_backend,
+    no_default_backend,
+)
+
+
+@field_operator
+def addition(a: cases.IJKField, b: cases.IJKField) -> cases.IJKField:
+    return a + b
+
+
+@field_operator
+def mixed_args(
+    a: cases.IJKField, b: np.float64, c: cases.IJKField
+) -> tuple[cases.IJKField, tuple[cases.IJKField, cases.IJKField]]:
+    return (a, (c, a))
+
+
+def test_allocate_default_unique(cartesian_case):  # noqa: F811 # fixtures
+    a = cases.allocate(cartesian_case, mixed_args, "a")()
+
+    assert np.min(a.array()) == 0
+    assert np.max(a.array()) == np.prod(tuple(cartesian_case.default_sizes.values())) - 1
+
+    b = cases.allocate(cartesian_case, mixed_args, "b")()
+
+    assert b == np.max(a.array()) + 1
+
+    c = cases.allocate(cartesian_case, mixed_args, "c")()
+
+    assert np.min(c.array()) == b + 1
+    assert np.max(c.array()) == np.prod(tuple(cartesian_case.default_sizes.values())) * 2
+
+
+def test_allocate_return_default_zeros(cartesian_case):  # noqa: F811 # fixtures
+    a, (b, c) = cases.allocate(cartesian_case, mixed_args, cases.RETURN)()
+
+    assert np.all(a.array() == 0)
+    assert np.all(a.array() == b.array())
+    assert np.all(b.array() == c.array())
+
+
+def test_allocate_const(cartesian_case):  # noqa: F811 # fixtures
+    a = cases.allocate(cartesian_case, mixed_args, "a").strategy(cases.ConstInitializer(42))()
+    assert np.all(a.array() == 42)
+
+    b = cases.allocate(cartesian_case, mixed_args, "b").strategy(cases.ConstInitializer(42))()
+    assert b == 42.0
+
+
+@pytest.mark.parametrize("fieldview_backend", [roundtrip.executor])
+def test_verify_fails_with_wrong_reference(cartesian_case):  # noqa: F811 # fixtures
+    a = cases.allocate(cartesian_case, addition, "a")()
+    b = cases.allocate(cartesian_case, addition, "b")()
+    out = cases.allocate(cartesian_case, addition, cases.RETURN)()
+    wrong_ref = a
+
+    with pytest.raises(AssertionError):
+        cases.verify(cartesian_case, addition, a, b, out=out, ref=wrong_ref)
+
+
+@pytest.mark.parametrize("fieldview_backend", [roundtrip.executor])
+def test_verify_fails_with_wrong_type(cartesian_case):  # noqa: F811 # fixtures
+    a = cases.allocate(cartesian_case, addition, "a").dtype(np.float32)()
+    b = cases.allocate(cartesian_case, addition, "b")()
+    out = cases.allocate(cartesian_case, addition, cases.RETURN)()
+
+    with pytest.raises(common.GTTypeError):
+        cases.verify(cartesian_case, addition, a, b, out=out, ref=a.array() + b.array())
+
+
+@pytest.mark.parametrize("fieldview_backend", [roundtrip.executor])
+def test_verify_with_default_data_fails_with_wrong_reference(
+    cartesian_case,  # noqa: F811 # fixtures
+):
+    def wrong_ref(a, b):
+        return a - b
+
+    with pytest.raises(AssertionError):
+        cases.verify_with_default_data(cartesian_case, addition, ref=wrong_ref)
diff --git a/tests/next_tests/unit_tests/iterator_tests/test_embedded_internals.py b/tests/next_tests/unit_tests/iterator_tests/test_embedded_internals.py
index 15b3d177..3a35570c 100644
--- a/tests/next_tests/unit_tests/iterator_tests/test_embedded_internals.py
+++ b/tests/next_tests/unit_tests/iterator_tests/test_embedded_internals.py
@@ -12,85 +12,138 @@
 #
 # SPDX-License-Identifier: GPL-3.0-or-later
 
+import contextvars as cvars
+import threading
+from typing import Any, Callable, Optional
+
 import numpy as np
 import pytest
 
-from gt4py.next.iterator.embedded import Column
+from gt4py.next.iterator import embedded
+
+
+def _run_within_context(
+    func: Callable[[], Any],
+    *,
+    column_range: Optional[range] = None,
+    offset_provider: Optional[embedded.OffsetProvider] = None,
+) -> Any:
+    def wrapped_func():
+        embedded.column_range_cvar.set(column_range)
+        embedded.offset_provider_cvar.set(offset_provider)
+        func()
+
+    cvars.copy_context().run(wrapped_func)
 
 
 def test_column_ufunc():
-    a = Column(1, np.asarray(range(0, 3)))
-    b = Column(1, np.asarray(range(3, 6)))
+    def test_func():
+        a = embedded.Column(1, np.asarray(range(0, 3)))
+        b = embedded.Column(1, np.asarray(range(3, 6)))
+        res = a + b
+
+        assert isinstance(res, embedded.Column)
+        assert np.array_equal(res.data, a.data + b.data)
+        assert res.kstart == 1
+
+    _run_within_context(test_func)
+
+    def test_func(data_a: int, data_b: int):
+        a = embedded.Column(1, data_a)
+        b = embedded.Column(1, data_b)
+        res = a + b
 
-    res = a + b
-    assert isinstance(res, Column)
-    assert np.array_equal(res.data, a.data + b.data)
-    assert res.kstart == 1
+        assert isinstance(res, embedded.Column)
+        assert np.array_equal(res.data, a.data + b.data)
+        assert res.kstart == 1
+
+    # Setting an invalid column_range here shouldn't affect other contexts
+    embedded.column_range_cvar.set(range(2, 999))
+    _run_within_context(lambda: test_func(2, 3), column_range=range(0, 3))
 
 
 def test_column_ufunc_with_scalar():
-    a = Column(1, np.asarray(range(0, 3)))
-    res = 1.0 + a
-    assert isinstance(res, Column)
-    assert np.array_equal(res.data, a.data + 1.0)
-    assert res.kstart == 1
+    def test_func():
+        a = embedded.Column(1, np.asarray(range(0, 3)))
+        res = 1.0 + a
+        assert isinstance(res, embedded.Column)
+        assert np.array_equal(res.data, a.data + 1.0)
+        assert res.kstart == 1
 
 
 def test_column_ufunc_wrong_kstart():
-    a = Column(1, np.asarray(range(0, 3)))
-    wrong_kstart = Column(2, np.asarray(range(3, 6)))
+    def test_func():
+        a = embedded.Column(1, np.asarray(range(0, 3)))
+        wrong_kstart = embedded.Column(2, np.asarray(range(3, 6)))
+
+        with pytest.raises(ValueError):
+            a + wrong_kstart
 
-    with pytest.raises(ValueError):
-        a + wrong_kstart
+    _run_within_context(test_func)
 
 
 def test_column_ufunc_wrong_shape():
-    a = Column(1, np.asarray(range(0, 3)))
-    wrong_shape = Column(1, np.asarray([1, 2]))
+    def test_func():
+        a = embedded.Column(1, np.asarray(range(0, 3)))
+        wrong_shape = embedded.Column(1, np.asarray([1, 2]))
 
-    with pytest.raises(ValueError):
-        a + wrong_shape
+        with pytest.raises(ValueError):
+            a + wrong_shape
+
+    _run_within_context(test_func)
 
 
 def test_column_array_function():
-    cond = Column(1, np.asarray([True, False]))
-    a = Column(1, np.asarray([1, 1]))
-    b = Column(1, np.asarray([2, 2]))
+    def test_func():
+        cond = embedded.Column(1, np.asarray([True, False]))
+        a = embedded.Column(1, np.asarray([1, 1]))
+        b = embedded.Column(1, np.asarray([2, 2]))
+
+        res = np.where(cond, a, b)
+        ref = np.asarray([1, 2])
 
-    res = np.where(cond, a, b)
-    ref = np.asarray([1, 2])
+        assert isinstance(res, embedded.Column)
+        assert np.array_equal(res.data, ref)
+        assert res.kstart == 1
 
-    assert isinstance(res, Column)
-    assert np.array_equal(res.data, ref)
-    assert res.kstart == 1
+    _run_within_context(test_func)
 
 
 def test_column_array_function_with_scalar():
-    cond = Column(1, np.asarray([True, False]))
-    a = 1
-    b = Column(1, np.asarray([2, 2]))
+    def test_func():
+        cond = embedded.Column(1, np.asarray([True, False]))
+        a = 1
+        b = embedded.Column(1, np.asarray([2, 2]))
 
-    res = np.where(cond, a, b)
-    ref = np.asarray([1, 2])
+        res = np.where(cond, a, b)
+        ref = np.asarray([1, 2])
 
-    assert isinstance(res, Column)
-    assert np.array_equal(res.data, ref)
-    assert res.kstart == 1
+        assert isinstance(res, embedded.Column)
+        assert np.array_equal(res.data, ref)
+        assert res.kstart == 1
+
+    _run_within_context(test_func)
 
 
 def test_column_array_function_wrong_kstart():
-    cond = Column(1, np.asarray([True, False]))
-    wrong_kstart = Column(2, np.asarray([1, 1]))
-    b = Column(1, np.asarray([2, 2]))
+    def test_func():
+        cond = embedded.Column(1, np.asarray([True, False]))
+        wrong_kstart = embedded.Column(2, np.asarray([1, 1]))
+        b = embedded.Column(1, np.asarray([2, 2]))
+
+        with pytest.raises(ValueError):
+            np.where(cond, wrong_kstart, b)
 
-    with pytest.raises(ValueError):
-        np.where(cond, wrong_kstart, b)
+    _run_within_context(test_func)
 
 
 def test_column_array_function_wrong_shape():
-    cond = Column(1, np.asarray([True, False]))
-    wrong_shape = Column(2, np.asarray([1, 1, 1]))
-    b = Column(1, np.asarray([2, 2]))
+    def test_func():
+        cond = embedded.Column(1, np.asarray([True, False]))
+        wrong_shape = embedded.Column(2, np.asarray([1, 1, 1]))
+        b = embedded.Column(1, np.asarray([2, 2]))
+
+        with pytest.raises(ValueError):
+            np.where(cond, wrong_shape, b)
 
-    with pytest.raises(ValueError):
-        np.where(cond, wrong_shape, b)
+    _run_within_context(test_func)