Repository: GridTools/gt4py
Pull request #1441 (2024-02-01 14:06:39)

feat[next]: Embedded support for skip value connectivities
- Introduces `common._DEFAULT_SKIP_VALUE = -1`, currently we don't support execution with skip_values which are different than that default.
- inverse_image works by ignoring the skip_value for the inverse image calculation. current result is the smallest hypercube that contains all non skip values that are requested.

Additional changes:
- Adds a ffront test for the fvm nabla example
- Fixes a lowering bug from past to itir, where domain didn't respect order of horizontal and vertical dimension



Patch:
diff --git a/src/gt4py/next/common.py b/src/gt4py/next/common.py
index fdf515d2..f4e35b55 100644
--- a/src/gt4py/next/common.py
+++ b/src/gt4py/next/common.py
@@ -732,7 +732,7 @@ class ConnectivityKind(enum.Flag):
 
 
 @extended_runtime_checkable
-# type: ignore[misc] # DimT should be covariant, but break in another place
+# type: ignore[misc] # DimT should be covariant, but breaks in another place
 class ConnectivityField(Field[DimsT, core_defs.IntegralScalar], Protocol[DimsT, DimT]):
     @property
     @abc.abstractmethod
@@ -749,6 +749,10 @@ class ConnectivityField(Field[DimsT, core_defs.IntegralScalar], Protocol[DimsT,
     @abc.abstractmethod
     def inverse_image(self, image_range: UnitRange | NamedRange) -> Sequence[NamedRange]: ...
 
+    @property
+    @abc.abstractmethod
+    def skip_value(self) -> Optional[core_defs.IntegralScalar]: ...
+
     # Operators
     def __abs__(self) -> Never:
         raise TypeError("'ConnectivityField' does not support this operation.")
@@ -840,6 +844,7 @@ def _connectivity(
     *,
     domain: Optional[DomainLike] = None,
     dtype: Optional[core_defs.DType] = None,
+    skip_value: Optional[core_defs.IntegralScalar] = None,
 ) -> ConnectivityField:
     raise NotImplementedError
 
@@ -918,6 +923,10 @@ class CartesianConnectivity(ConnectivityField[DimsT, DimT]):
     def codomain(self) -> DimT:
         return self.dimension
 
+    @property
+    def skip_value(self) -> None:
+        return None
+
     @functools.cached_property
     def kind(self) -> ConnectivityKind:
         return ConnectivityKind(0)
@@ -1083,4 +1092,4 @@ class FieldBuiltinFuncRegistry:
 #: Numeric value used to represent missing values in connectivities.
 #: Equivalent to the `_FillValue` attribute in the UGRID Conventions
 #: (see: http://ugrid-conventions.github.io/ugrid-conventions/).
-SKIP_VALUE: Final[int] = -1
+_DEFAULT_SKIP_VALUE: Final[int] = -1
diff --git a/src/gt4py/next/constructors.py b/src/gt4py/next/constructors.py
index 8b41bf7c..f8e7b9bf 100644
--- a/src/gt4py/next/constructors.py
+++ b/src/gt4py/next/constructors.py
@@ -308,6 +308,7 @@ def as_connectivity(
     *,
     allocator: Optional[next_allocators.FieldBufferAllocatorProtocol] = None,
     device: Optional[core_defs.Device] = None,
+    skip_value: Optional[core_defs.IntegralScalar] = None,
     # copy=False, TODO
 ) -> common.ConnectivityField:
     """
@@ -330,6 +331,9 @@ def as_connectivity(
     Raises:
         ValueError: If the domain or codomain is invalid, or if the shape of the data does not match the domain shape.
     """
+    assert (
+        skip_value is None or skip_value == common._DEFAULT_SKIP_VALUE
+    )  # TODO(havogt): not yet configurable
     if isinstance(domain, Sequence) and all(isinstance(dim, common.Dimension) for dim in domain):
         domain = cast(Sequence[common.Dimension], domain)
         if len(domain) != data.ndim:
@@ -359,7 +363,7 @@ def as_connectivity(
     # TODO(havogt): consider adding MutableNDArrayObject
     buffer.ndarray[...] = storage_utils.asarray(data)  # type: ignore[index]
     connectivity_field = common._connectivity(
-        buffer.ndarray, codomain=codomain, domain=actual_domain
+        buffer.ndarray, codomain=codomain, domain=actual_domain, skip_value=skip_value
     )
     assert isinstance(connectivity_field, nd_array_field.NdArrayConnectivityField)
 
diff --git a/src/gt4py/next/embedded/nd_array_field.py b/src/gt4py/next/embedded/nd_array_field.py
index c39408ba..1bdb7161 100644
--- a/src/gt4py/next/embedded/nd_array_field.py
+++ b/src/gt4py/next/embedded/nd_array_field.py
@@ -26,8 +26,9 @@ from numpy import typing as npt
 from gt4py._core import definitions as core_defs
 from gt4py.eve.extended_typing import Any, Never, Optional, ParamSpec, TypeAlias, TypeVar
 from gt4py.next import common
-from gt4py.next.embedded import common as embedded_common
+from gt4py.next.embedded import common as embedded_common, context as embedded_context
 from gt4py.next.ffront import fbuiltins
+from gt4py.next.iterator import embedded as itir_embedded
 
 
 try:
@@ -170,9 +171,12 @@ class NdArrayField(
         if not common.is_connectivity_field(connectivity):
             assert isinstance(connectivity, fbuiltins.FieldOffset)
             connectivity = connectivity.as_connectivity_field()
-
         assert common.is_connectivity_field(connectivity)
 
+        # Current implementation relies on skip_value == -1:
+        # if we assume the indexed array has at least one element, we wrap around without out of bounds
+        assert connectivity.skip_value is None or connectivity.skip_value == -1
+
         # Compute the new domain
         dim = connectivity.codomain
         dim_idx = self.domain.dim_index(dim)
@@ -315,6 +319,7 @@ class NdArrayConnectivityField(  # type: ignore[misc] # for __ne__, __eq__
     NdArrayField[common.DimsT, core_defs.IntegralScalar],
 ):
     _codomain: common.DimT
+    _skip_value: Optional[core_defs.IntegralScalar]
 
     @functools.cached_property
     def _cache(self) -> dict:
@@ -329,6 +334,10 @@ class NdArrayConnectivityField(  # type: ignore[misc] # for __ne__, __eq__
     def codomain(self) -> common.DimT:
         return self._codomain
 
+    @property
+    def skip_value(self) -> Optional[core_defs.IntegralScalar]:
+        return self._skip_value
+
     @functools.cached_property
     def kind(self) -> common.ConnectivityKind:
         kind = common.ConnectivityKind.MODIFY_STRUCTURE
@@ -349,6 +358,7 @@ class NdArrayConnectivityField(  # type: ignore[misc] # for __ne__, __eq__
         *,
         domain: common.DomainLike,
         dtype: Optional[core_defs.DTypeLike] = None,
+        skip_value: Optional[core_defs.IntegralScalar] = None,
     ) -> NdArrayConnectivityField:
         domain = common.domain(domain)
         xp = cls.array_ns
@@ -367,7 +377,12 @@ class NdArrayConnectivityField(  # type: ignore[misc] # for __ne__, __eq__
 
         assert isinstance(codomain, common.Dimension)
 
-        return cls(domain, array, codomain)
+        return cls(
+            domain,
+            array,
+            codomain,
+            _skip_value=skip_value,
+        )
 
     def inverse_image(
         self, image_range: common.UnitRange | common.NamedRange
@@ -390,47 +405,16 @@ class NdArrayConnectivityField(  # type: ignore[misc] # for __ne__, __eq__
             assert isinstance(image_range, common.UnitRange)
 
             assert common.UnitRange.is_finite(image_range)
-            restricted_mask = (self._ndarray >= image_range.start) & (
-                self._ndarray < image_range.stop
-            )
-            # indices of non-zero elements in each dimension
-            nnz: tuple[core_defs.NDArrayObject, ...] = xp.nonzero(restricted_mask)
-
-            new_dims = []
-            non_contiguous_dims = []
-
-            for i, dim_nnz_indices in enumerate(nnz):
-                # Check if the indices are contiguous
-                first_data_index = dim_nnz_indices[0]
-                assert isinstance(first_data_index, core_defs.INTEGRAL_TYPES)
-                last_data_index = dim_nnz_indices[-1]
-                assert isinstance(last_data_index, core_defs.INTEGRAL_TYPES)
-                indices, counts = xp.unique(dim_nnz_indices, return_counts=True)
-                dim_range = self._domain[i]
-
-                if len(xp.unique(counts)) == 1 and (
-                    len(indices) == last_data_index - first_data_index + 1
-                ):
-                    idx_offset = dim_range[1].start
-                    start = idx_offset + first_data_index
-                    assert common.is_int_index(start)
-                    stop = idx_offset + last_data_index + 1
-                    assert common.is_int_index(stop)
-                    new_dims.append(
-                        common.named_range(
-                            (
-                                dim_range[0],
-                                (start, stop),
-                            )
-                        )
-                    )
-                else:
-                    non_contiguous_dims.append(dim_range[0])
 
-            if non_contiguous_dims:
-                raise ValueError(
-                    f"Restriction generates non-contiguous dimensions '{non_contiguous_dims}'."
-                )
+            relative_ranges = _hypercube(self._ndarray, image_range, xp, self.skip_value)
+
+            if relative_ranges is None:
+                raise ValueError("Restriction generates non-contiguous dimensions.")
+
+            new_dims = [
+                common.named_range((d, rr + ar.start))
+                for d, ar, rr in zip(self.domain.dims, self.domain.ranges, relative_ranges)
+            ]
 
             self._cache[cache_key] = new_dims
 
@@ -444,7 +428,7 @@ class NdArrayConnectivityField(  # type: ignore[misc] # for __ne__, __eq__
             xp = cls.array_ns
             new_domain, buffer_slice = self._slice(index)
             new_buffer = xp.asarray(self.ndarray[buffer_slice])
-            restricted_connectivity = cls(new_domain, new_buffer, self.codomain)
+            restricted_connectivity = cls(new_domain, new_buffer, self.codomain, self.skip_value)
             self._cache[cache_key] = restricted_connectivity
 
         return restricted_connectivity
@@ -452,6 +436,41 @@ class NdArrayConnectivityField(  # type: ignore[misc] # for __ne__, __eq__
     __getitem__ = restrict
 
 
+def _hypercube(
+    index_array: core_defs.NDArrayObject,
+    image_range: common.UnitRange,
+    xp: ModuleType,
+    skip_value: Optional[core_defs.IntegralScalar] = None,
+) -> Optional[list[common.UnitRange]]:
+    """
+    Return the hypercube that contains all indices in `index_array` that are within `image_range`, or `None` if no such hypercube exists.
+
+    If `skip_value` is given, the selected values are ignored. It returns the smallest hypercube.
+    A bigger hypercube could be constructed by adding lines that contain only `skip_value`s.
+    Example:
+    index_array =  0  1 -1
+                   3  4 -1
+                  -1 -1 -1
+    skip_value = -1
+    would currently select the 2x2 range [0,2], [0,2], but could also select the 3x3 range [0,3], [0,3].
+    """
+    select_mask = (index_array >= image_range.start) & (index_array < image_range.stop)
+
+    nnz: tuple[core_defs.NDArrayObject, ...] = xp.nonzero(select_mask)
+
+    slices = tuple(
+        slice(xp.min(dim_nnz_indices), xp.max(dim_nnz_indices) + 1) for dim_nnz_indices in nnz
+    )
+    hcube = select_mask[tuple(slices)]
+    if skip_value is not None:
+        ignore_mask = index_array == skip_value
+        hcube |= ignore_mask[tuple(slices)]
+    if not xp.all(hcube):
+        return None
+
+    return [common.UnitRange(s.start, s.stop) for s in slices]
+
+
 # -- Specialized implementations for builtin operations on array fields --
 
 NdArrayField.register_builtin_func(
@@ -483,21 +502,49 @@ NdArrayField.register_builtin_func(
 NdArrayField.register_builtin_func(fbuiltins.where, _make_builtin("where", "where"))
 
 
-def _make_reduction(builtin_name: str, array_builtin_name: str) -> Callable[
+def _make_reduction(
+    builtin_name: str, array_builtin_name: str, initial_value_op: Callable
+) -> Callable[
     ...,
     NdArrayField[common.DimsT, core_defs.ScalarT],
 ]:
     def _builtin_op(
         field: NdArrayField[common.DimsT, core_defs.ScalarT], axis: common.Dimension
     ) -> NdArrayField[common.DimsT, core_defs.ScalarT]:
+        xp = field.array_ns
+
         if not axis.kind == common.DimensionKind.LOCAL:
             raise ValueError("Can only reduce local dimensions.")
         if axis not in field.domain.dims:
             raise ValueError(f"Field can not be reduced as it doesn't have dimension '{axis}'.")
+        if len([d for d in field.domain.dims if d.kind is common.DimensionKind.LOCAL]) > 1:
+            raise NotImplementedError(
+                "Reducing a field with more than one local dimension is not supported."
+            )
         reduce_dim_index = field.domain.dims.index(axis)
+        current_offset_provider = embedded_context.offset_provider.get(None)
+        assert current_offset_provider is not None
+        offset_definition = current_offset_provider[
+            axis.value
+        ]  # assumes offset and local dimension have same name
+        assert isinstance(offset_definition, itir_embedded.NeighborTableOffsetProvider)
         new_domain = common.Domain(*[nr for nr in field.domain if nr[0] != axis])
+
+        broadcast_slice = tuple(
+            slice(None) if d in [axis, offset_definition.origin_axis] else xp.newaxis
+            for d in field.domain.dims
+        )
+        masked_array = xp.where(
+            xp.asarray(offset_definition.table[broadcast_slice]) != common._DEFAULT_SKIP_VALUE,
+            field.ndarray,
+            initial_value_op(field),
+        )
+
         return field.__class__.from_array(
-            getattr(field.array_ns, array_builtin_name)(field.ndarray, axis=reduce_dim_index),
+            getattr(xp, array_builtin_name)(
+                masked_array,
+                axis=reduce_dim_index,
+            ),
             domain=new_domain,
         )
 
@@ -505,9 +552,15 @@ def _make_reduction(builtin_name: str, array_builtin_name: str) -> Callable[
     return _builtin_op
 
 
-NdArrayField.register_builtin_func(fbuiltins.neighbor_sum, _make_reduction("neighbor_sum", "sum"))
-NdArrayField.register_builtin_func(fbuiltins.max_over, _make_reduction("max_over", "max"))
-NdArrayField.register_builtin_func(fbuiltins.min_over, _make_reduction("min_over", "min"))
+NdArrayField.register_builtin_func(
+    fbuiltins.neighbor_sum, _make_reduction("neighbor_sum", "sum", lambda x: x.dtype.scalar_type(0))
+)
+NdArrayField.register_builtin_func(
+    fbuiltins.max_over, _make_reduction("max_over", "max", lambda x: x.array_ns.min(x._ndarray))
+)
+NdArrayField.register_builtin_func(
+    fbuiltins.min_over, _make_reduction("min_over", "min", lambda x: x.array_ns.max(x._ndarray))
+)
 
 
 # -- Concrete array implementations --
diff --git a/src/gt4py/next/ffront/decorator.py b/src/gt4py/next/ffront/decorator.py
index a556d0ea..0cf1611b 100644
--- a/src/gt4py/next/ffront/decorator.py
+++ b/src/gt4py/next/ffront/decorator.py
@@ -221,6 +221,10 @@ class Program:
                 f"The following closure variables are undefined: {', '.join(undefined_symbols)}."
             )
 
+    @property
+    def __name__(self) -> str:
+        return self.definition.__name__
+
     @functools.cached_property
     def __gt_allocator__(
         self,
@@ -603,6 +607,10 @@ class FieldOperator(GTCallable, Generic[OperatorNodeT]):
             operator_attributes=operator_attributes,
         )
 
+    @property
+    def __name__(self) -> str:
+        return self.definition.__name__
+
     def __gt_type__(self) -> ts.CallableType:
         type_ = self.foast_node.type
         assert isinstance(type_, ts.CallableType)
diff --git a/src/gt4py/next/ffront/fbuiltins.py b/src/gt4py/next/ffront/fbuiltins.py
index 493493f6..34251a36 100644
--- a/src/gt4py/next/ffront/fbuiltins.py
+++ b/src/gt4py/next/ffront/fbuiltins.py
@@ -360,12 +360,14 @@ class FieldOffset(runtime.Offset):
             if common.is_connectivity_field(offset_definition):
                 connectivity = offset_definition
             elif isinstance(offset_definition, gtx.NeighborTableOffsetProvider):
-                assert not offset_definition.has_skip_values
                 connectivity = gtx.as_connectivity(
                     domain=self.target,
                     codomain=self.source,
                     data=offset_definition.table,
                     dtype=offset_definition.index_type,
+                    skip_value=(
+                        common._DEFAULT_SKIP_VALUE if offset_definition.has_skip_values else None
+                    ),
                 )
             else:
                 raise NotImplementedError()
diff --git a/src/gt4py/next/ffront/past_to_itir.py b/src/gt4py/next/ffront/past_to_itir.py
index 8be93096..e092538c 100644
--- a/src/gt4py/next/ffront/past_to_itir.py
+++ b/src/gt4py/next/ffront/past_to_itir.py
@@ -231,7 +231,6 @@ class ProgramLowering(
         node_domain: Optional[past.Expr],
         slices: Optional[list[past.Slice]] = None,
     ) -> itir.FunCall:
-        domain_args = []
 
         assert isinstance(out_field.type, ts.TypeSpec)
         out_field_types = type_info.primitive_constituents(out_field.type).to_list()
@@ -246,6 +245,8 @@ class ProgramLowering(
                 " caught in type deduction already."
             )
 
+        domain_args = []
+        domain_args_kind = []
         for dim_i, dim in enumerate(out_dims):
             # an expression for the size of a dimension
             dim_size = itir.SymRef(id=_size_arg_from_field(out_field.id, dim_i))
@@ -271,11 +272,17 @@ class ProgramLowering(
                     args=[itir.AxisLiteral(value=dim.value), lower, upper],
                 )
             )
+            domain_args_kind.append(dim.kind)
 
         if self.grid_type == GridType.CARTESIAN:
             domain_builtin = "cartesian_domain"
         elif self.grid_type == GridType.UNSTRUCTURED:
             domain_builtin = "unstructured_domain"
+            # for no good reason, the domain arguments for unstructured need to be in order (horizontal, vertical)
+            if domain_args_kind[0] == DimensionKind.VERTICAL:
+                assert len(domain_args) == 2
+                assert domain_args_kind[1] == DimensionKind.HORIZONTAL
+                domain_args[0], domain_args[1] = domain_args[1], domain_args[0]
         else:
             raise AssertionError()
 
diff --git a/src/gt4py/next/iterator/atlas_utils.py b/src/gt4py/next/iterator/atlas_utils.py
index 500c9253..a23a4b51 100644
--- a/src/gt4py/next/iterator/atlas_utils.py
+++ b/src/gt4py/next/iterator/atlas_utils.py
@@ -17,6 +17,8 @@ try:
 except ImportError:
     IrregularConnectivity = None
 
+from gt4py.next import common
+
 
 # TODO(tehrengruber): make this a proper Connectivity instead of faking a numpy array
 class AtlasTable:
@@ -29,7 +31,7 @@ class AtlasTable:
             if neigh_index < self.atlas_connectivity.cols(primary_index):
                 return self.atlas_connectivity[primary_index, neigh_index]
             else:
-                return None
+                return common._DEFAULT_SKIP_VALUE
         else:
             if neigh_index < 2:
                 return self.atlas_connectivity[primary_index, neigh_index]
@@ -53,3 +55,12 @@ class AtlasTable:
                 if v is not None:
                     maximum = max(maximum, v)
         return maximum
+
+    def asnumpy(self):
+        import numpy as np
+
+        res = np.empty(self.shape, dtype=self.dtype)
+        for i in range(self.shape[0]):
+            for j in range(self.shape[1]):
+                res[i, j] = self[i, j]
+        return res
diff --git a/src/gt4py/next/iterator/embedded.py b/src/gt4py/next/iterator/embedded.py
index a45b81a7..80fc5392 100644
--- a/src/gt4py/next/iterator/embedded.py
+++ b/src/gt4py/next/iterator/embedded.py
@@ -532,7 +532,7 @@ def execute_shift(
                 assert common.is_int_index(cur_index)
                 if offset_implementation.mapped_index(cur_index, index) in [
                     None,
-                    common.SKIP_VALUE,
+                    common._DEFAULT_SKIP_VALUE,
                 ]:
                     return None
 
@@ -559,7 +559,7 @@ def execute_shift(
         assert common.is_int_index(cur_index)
         if offset_implementation.mapped_index(cur_index, index) in [
             None,
-            common.SKIP_VALUE,
+            common._DEFAULT_SKIP_VALUE,
         ]:
             return None
         else:
diff --git a/src/gt4py/next/program_processors/runners/dace_iterator/__init__.py b/src/gt4py/next/program_processors/runners/dace_iterator/__init__.py
index aba56561..1263cff5 100644
--- a/src/gt4py/next/program_processors/runners/dace_iterator/__init__.py
+++ b/src/gt4py/next/program_processors/runners/dace_iterator/__init__.py
@@ -22,7 +22,6 @@ import numpy as np
 from dace.codegen.compiled_sdfg import CompiledSDFG
 from dace.sdfg import utils as sdutils
 from dace.transformation.auto import auto_optimize as autoopt
-from dace.transformation.interstate import RefineNestedAccess
 
 import gt4py.next.allocators as next_allocators
 import gt4py.next.iterator.ir as itir
@@ -53,8 +52,19 @@ def get_sorted_dim_ranges(domain: common.Domain) -> Sequence[common.FiniteUnitRa
 _build_type = "Release"
 
 
-def convert_arg(arg: Any):
+def convert_arg(arg: Any, sdfg_param: str):
     if common.is_field(arg):
+        # field domain offsets are not supported
+        non_zero_offsets = [
+            (dim, dim_range)
+            for dim, dim_range in zip(arg.domain.dims, arg.domain.ranges)
+            if dim_range.start != 0
+        ]
+        if non_zero_offsets:
+            dim, dim_range = non_zero_offsets[0]
+            raise RuntimeError(
+                f"Field '{sdfg_param}' passed as array slice with offset {dim_range.start} on dimension {dim.value}."
+            )
         sorted_dims = get_sorted_dims(arg.domain.dims)
         ndim = len(sorted_dims)
         dim_indices = [dim_index for dim_index, _ in sorted_dims]
@@ -99,7 +109,7 @@ def preprocess_program(
 
 def get_args(sdfg: dace.SDFG, args: Sequence[Any]) -> dict[str, Any]:
     sdfg_params: Sequence[str] = sdfg.arg_names
-    return {sdfg_param: convert_arg(arg) for sdfg_param, arg in zip(sdfg_params, args)}
+    return {sdfg_param: convert_arg(arg, sdfg_param) for sdfg_param, arg in zip(sdfg_params, args)}
 
 
 def _ensure_is_on_device(
@@ -140,27 +150,6 @@ def get_shape_args(
     return shape_args
 
 
-def get_offset_args(sdfg: dace.SDFG, args: Sequence[Any]) -> Mapping[str, int]:
-    sdfg_arrays: Mapping[str, dace.data.Array] = sdfg.arrays
-    sdfg_params: Sequence[str] = sdfg.arg_names
-    field_args = {param: arg for param, arg in zip(sdfg_params, args) if common.is_field(arg)}
-
-    # assume that arrays for connectivity tables do not use offset
-    assert all(
-        drange.start == 0
-        for sdfg_param, arg in field_args.items()
-        if sdfg_param.startswith("__connectivity")
-        for drange in arg.domain.ranges
-    )
-
-    return {
-        str(sym): -drange.start
-        for sdfg_param, arg in field_args.items()
-        if not sdfg_param.startswith("__connectivity")
-        for sym, drange in zip(sdfg_arrays[sdfg_param].offset, get_sorted_dim_ranges(arg.domain))
-    }
-
-
 def get_stride_args(
     arrays: Mapping[str, dace.data.Array], args: Mapping[str, Any]
 ) -> Mapping[str, int]:
@@ -243,7 +232,6 @@ def get_sdfg_args(sdfg: dace.SDFG, *args, check_args: bool = False, **kwargs) ->
     dace_conn_shapes = get_shape_args(sdfg.arrays, dace_conn_args)
     dace_strides = get_stride_args(sdfg.arrays, dace_field_args)
     dace_conn_strides = get_stride_args(sdfg.arrays, dace_conn_args)
-    dace_offsets = get_offset_args(sdfg, args)
     all_args = {
         **dace_args,
         **dace_conn_args,
@@ -251,7 +239,6 @@ def get_sdfg_args(sdfg: dace.SDFG, *args, check_args: bool = False, **kwargs) ->
         **dace_conn_shapes,
         **dace_strides,
         **dace_conn_strides,
-        **dace_offsets,
     }
 
     if check_args:
@@ -308,9 +295,6 @@ def build_sdfg_from_itir(
     sdfg = sdfg_genenerator.visit(program)
     if sdfg is None:
         raise RuntimeError(f"Visit failed for program {program.id}.")
-    elif tmps:
-        # This pass is needed to avoid transformation errors in SDFG inlining, because temporaries are using offsets
-        sdfg.apply_transformations_repeated(RefineNestedAccess)
 
     for nested_sdfg in sdfg.all_sdfgs_recursive():
         if not nested_sdfg.debuginfo:
diff --git a/src/gt4py/next/program_processors/runners/dace_iterator/itir_to_sdfg.py b/src/gt4py/next/program_processors/runners/dace_iterator/itir_to_sdfg.py
index 05a10444..ed71522d 100644
--- a/src/gt4py/next/program_processors/runners/dace_iterator/itir_to_sdfg.py
+++ b/src/gt4py/next/program_processors/runners/dace_iterator/itir_to_sdfg.py
@@ -11,6 +11,7 @@
 # distribution for a copy of the license or check <https://www.gnu.org/licenses/>.
 #
 # SPDX-License-Identifier: GPL-3.0-or-later
+import warnings
 from typing import Any, Mapping, Optional, Sequence, cast
 
 import dace
@@ -104,7 +105,7 @@ def _make_array_shape_and_strides(
     dims: Sequence[Dimension],
     neighbor_tables: Mapping[str, NeighborTable],
     sort_dims: bool,
-) -> tuple[list[dace.symbol], list[dace.symbol], list[dace.symbol]]:
+) -> tuple[list[dace.symbol], list[dace.symbol]]:
     """
     Parse field dimensions and allocate symbols for array shape and strides.
 
@@ -127,9 +128,8 @@ def _make_array_shape_and_strides(
         )
         for i, dim in sorted_dims
     ]
-    offset = [dace.symbol(unique_name(f"{name}_offset{i}"), dtype) for i, _ in sorted_dims]
     strides = [dace.symbol(unique_name(f"{name}_stride{i}"), dtype) for i, _ in sorted_dims]
-    return shape, offset, strides
+    return shape, strides
 
 
 def _check_no_lifts(node: itir.StencilClosure):
@@ -177,11 +177,10 @@ class ItirToSDFG(eve.NodeVisitor):
         name: str,
         type_: ts.TypeSpec,
         neighbor_tables: Mapping[str, NeighborTable],
-        has_offset: bool = True,
         sort_dimensions: bool = True,
     ):
         if isinstance(type_, ts.FieldType):
-            shape, offset, strides = _make_array_shape_and_strides(
+            shape, strides = _make_array_shape_and_strides(
                 name, type_.dims, neighbor_tables, sort_dimensions
             )
             dtype = as_dace_type(type_.dtype)
@@ -189,7 +188,6 @@ class ItirToSDFG(eve.NodeVisitor):
                 name,
                 shape=shape,
                 strides=strides,
-                offset=(offset if has_offset else None),
                 dtype=dtype,
             )
 
@@ -249,22 +247,27 @@ class ItirToSDFG(eve.NodeVisitor):
             # Another option, in the future, is to use symbolic strides and apply auto-tuning or some heuristics
             # to assign optimal stride values.
             tmp_shape, _ = new_array_symbols(tmp_name, len(dims))
-            tmp_offset = [
-                dace.symbol(unique_name(f"{tmp_name}_offset{i}")) for i in range(len(dims))
-            ]
             _, tmp_array = program_sdfg.add_array(
-                tmp_name, tmp_shape, as_dace_type(type_.dtype), offset=tmp_offset, transient=True
+                tmp_name, tmp_shape, as_dace_type(type_.dtype), transient=True
             )
 
             # Loop through all dimensions to visit the symbolic expressions for array shape and offset.
             # These expressions are later mapped to interstate symbols.
-            for (_, (begin, end)), offset_sym, shape_sym in zip(
+            for (_, (begin, end)), shape_sym in zip(
                 tmp_domain,
-                tmp_array.offset,
                 tmp_array.shape,
             ):
-                tmp_symbols[str(offset_sym)] = f"0 - {begin.value}"
-                tmp_symbols[str(shape_sym)] = f"{end.value} - {begin.value}"
+                """
+                The temporary field has a dimension range defined by `begin` and `end` values.
+                Therefore, the actual size is given by the difference `end.value - begin.value`.
+                Instead of allocating the actual size, we allocate space to enable indexing from 0
+                because we want to avoid using dace array offsets (which will be deprecated soon).
+                The result should still be valid, but the stencil will be using only a subset
+                of the array.
+                """
+                if not (isinstance(begin, SymbolExpr) and begin.value == "0"):
+                    warnings.warn(f"Domain start offset for temporary {tmp_name} is ignored.")
+                tmp_symbols[str(shape_sym)] = end.value
 
         return tmp_symbols
 
@@ -319,7 +322,6 @@ class ItirToSDFG(eve.NodeVisitor):
                 connectivity_identifier(offset),
                 type_,
                 neighbor_tables,
-                has_offset=False,
                 sort_dimensions=False,
             )
 
@@ -663,7 +665,6 @@ class ItirToSDFG(eve.NodeVisitor):
             output_name,
             shape=(array_table[node.output.id].shape[scan_dim_index],),
             strides=(array_table[node.output.id].strides[scan_dim_index],),
-            offset=(array_table[node.output.id].offset[scan_dim_index],),
             dtype=array_table[node.output.id].dtype,
         )
 
@@ -770,7 +771,7 @@ class ItirToSDFG(eve.NodeVisitor):
 
     def _visit_domain(
         self, node: itir.FunCall, context: Context
-    ) -> tuple[tuple[str, tuple[ValueExpr, ValueExpr]], ...]:
+    ) -> tuple[tuple[str, tuple[SymbolExpr | ValueExpr, SymbolExpr | ValueExpr]], ...]:
         assert isinstance(node.fun, itir.SymRef)
         assert node.fun.id == "cartesian_domain" or node.fun.id == "unstructured_domain"
 
diff --git a/src/gt4py/next/program_processors/runners/dace_iterator/itir_to_tasklet.py b/src/gt4py/next/program_processors/runners/dace_iterator/itir_to_tasklet.py
index 3a33ee1e..cf6d7ab0 100644
--- a/src/gt4py/next/program_processors/runners/dace_iterator/itir_to_tasklet.py
+++ b/src/gt4py/next/program_processors/runners/dace_iterator/itir_to_tasklet.py
@@ -22,7 +22,7 @@ import numpy as np
 import gt4py.eve.codegen
 from gt4py import eve
 from gt4py.next import Dimension, StridedNeighborOffsetProvider, type_inference as next_typing
-from gt4py.next.common import SKIP_VALUE as neighbor_skip_value
+from gt4py.next.common import _DEFAULT_SKIP_VALUE as neighbor_skip_value
 from gt4py.next.iterator import ir as itir, type_inference as itir_typing
 from gt4py.next.iterator.embedded import NeighborTableOffsetProvider
 from gt4py.next.iterator.ir import FunCall, Lambda


Test patch:
diff --git a/tests/next_tests/definitions.py b/tests/next_tests/definitions.py
index 56b220e0..a8ba72f3 100644
--- a/tests/next_tests/definitions.py
+++ b/tests/next_tests/definitions.py
@@ -170,7 +170,6 @@ EMBEDDED_SKIP_LIST = [
         XFAIL,
         UNSUPPORTED_MESSAGE,
     ),  # we can't extract the field type from scan args
-    (USES_MESH_WITH_SKIP_VALUES, XFAIL, UNSUPPORTED_MESSAGE),
 ]
 GTFN_SKIP_TEST_LIST = COMMON_SKIP_TEST_LIST + [
     # floordiv not yet supported, see https://github.com/GridTools/gt4py/issues/1136
diff --git a/tests/next_tests/integration_tests/feature_tests/ffront_tests/ffront_test_utils.py b/tests/next_tests/integration_tests/feature_tests/ffront_tests/ffront_test_utils.py
index d8c46960..aca601d7 100644
--- a/tests/next_tests/integration_tests/feature_tests/ffront_tests/ffront_test_utils.py
+++ b/tests/next_tests/integration_tests/feature_tests/ffront_tests/ffront_test_utils.py
@@ -239,12 +239,12 @@ def skip_value_mesh() -> MeshDescriptor:
 
     v2e_arr = np.array(
         [
-            [1, 8, 7, 0, -1],
-            [2, 8, 1, -1, -1],
-            [3, 9, 8, 2, -1],
-            [4, 10, 3, -1, -1],
-            [5, 11, 4, -1, -1],
-            [0, 6, 4, -1, -1],
+            [1, 8, 7, 0, common._DEFAULT_SKIP_VALUE],
+            [2, 8, 1, common._DEFAULT_SKIP_VALUE, common._DEFAULT_SKIP_VALUE],
+            [3, 9, 8, 2, common._DEFAULT_SKIP_VALUE],
+            [4, 10, 3, common._DEFAULT_SKIP_VALUE, common._DEFAULT_SKIP_VALUE],
+            [5, 11, 4, common._DEFAULT_SKIP_VALUE, common._DEFAULT_SKIP_VALUE],
+            [0, 6, 4, common._DEFAULT_SKIP_VALUE, common._DEFAULT_SKIP_VALUE],
             [6, 7, 9, 10, 11],
         ],
         dtype=gtx.IndexType,
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 e499f83f..5ded38ab 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
@@ -530,7 +530,7 @@ def test_nested_reduction(unstructured_case):
                 initial=0,
             )[unstructured_case.offset_provider["V2E"].table],
             axis=1,
-            where=unstructured_case.offset_provider["V2E"].table != common.SKIP_VALUE,
+            where=unstructured_case.offset_provider["V2E"].table != common._DEFAULT_SKIP_VALUE,
         ),
         comparison=lambda a, tmp_2: np.all(a == tmp_2),
     )
@@ -595,7 +595,7 @@ def test_tuple_with_local_field_in_reduction_shifted(unstructured_case):
             e[v2e.table] + np.tile(v, (v2e.max_neighbors, 1)).T,
             axis=1,
             initial=0,
-            where=v2e.table != common.SKIP_VALUE,
+            where=v2e.table != common._DEFAULT_SKIP_VALUE,
         )[unstructured_case.offset_provider["E2V"].table[:, 0]],
     )
 
@@ -735,7 +735,7 @@ def test_ternary_builtin_neighbor_sum(unstructured_case):
                 b[v2e_table],
                 axis=1,
                 initial=0,
-                where=v2e_table != common.SKIP_VALUE,
+                where=v2e_table != common._DEFAULT_SKIP_VALUE,
             )
         ),
     )
@@ -1079,9 +1079,9 @@ def test_where_k_offset(cartesian_case):
     )()
     out = cases.allocate(cartesian_case, fieldop_where_k_offset, "inp")()
 
-    ref = np.where(k_index.asnumpy() > 0, np.roll(inp.asnumpy(), 1, axis=1), 2)
+    ref = np.where(k_index.asnumpy() > 0, np.roll(inp.asnumpy(), 1, axis=1), out.asnumpy())
 
-    cases.verify(cartesian_case, prog, inp, k_index, out=out[:, 1:], ref=ref[:, 1:])
+    cases.verify(cartesian_case, prog, inp, k_index, out=out, ref=ref)
 
 
 def test_undefined_symbols(cartesian_case):
diff --git a/tests/next_tests/integration_tests/feature_tests/ffront_tests/test_external_local_field.py b/tests/next_tests/integration_tests/feature_tests/ffront_tests/test_external_local_field.py
index 569d7b56..16d40384 100644
--- a/tests/next_tests/integration_tests/feature_tests/ffront_tests/test_external_local_field.py
+++ b/tests/next_tests/integration_tests/feature_tests/ffront_tests/test_external_local_field.py
@@ -54,7 +54,7 @@ def test_external_local_field(unstructured_case):
             v2e_table,
             axis=1,
             initial=0,
-            where=v2e_table != common.SKIP_VALUE,
+            where=v2e_table != common._DEFAULT_SKIP_VALUE,
         ),
     )
 
diff --git a/tests/next_tests/integration_tests/feature_tests/ffront_tests/test_gt4py_builtins.py b/tests/next_tests/integration_tests/feature_tests/ffront_tests/test_gt4py_builtins.py
index 05824fa7..583246f5 100644
--- a/tests/next_tests/integration_tests/feature_tests/ffront_tests/test_gt4py_builtins.py
+++ b/tests/next_tests/integration_tests/feature_tests/ffront_tests/test_gt4py_builtins.py
@@ -63,7 +63,7 @@ def test_maxover_execution_(unstructured_case, strategy):
         inp.asnumpy()[v2e_table],
         axis=1,
         initial=np.min(inp.asnumpy()),
-        where=v2e_table != common.SKIP_VALUE,
+        where=v2e_table != common._DEFAULT_SKIP_VALUE,
     )
     cases.verify(unstructured_case, testee, inp, ref=ref, out=out)
 
@@ -83,31 +83,65 @@ def test_minover_execution(unstructured_case):
             edge_f[v2e_table],
             axis=1,
             initial=np.max(edge_f),
-            where=v2e_table != common.SKIP_VALUE,
+            where=v2e_table != common._DEFAULT_SKIP_VALUE,
         ),
     )
 
 
-@pytest.mark.uses_unstructured_shift
-def test_reduction_execution(unstructured_case):
-    @gtx.field_operator
-    def reduction(edge_f: cases.EField) -> cases.VField:
-        return neighbor_sum(edge_f(V2E), axis=V2EDim)
+@gtx.field_operator
+def reduction_e_field(edge_f: cases.EField) -> cases.VField:
+    return neighbor_sum(edge_f(V2E), axis=V2EDim)
 
-    @gtx.program
-    def fencil(edge_f: cases.EField, out: cases.VField):
-        reduction(edge_f, out=out)
 
+@gtx.field_operator
+def reduction_ek_field(
+    edge_f: common.Field[[Edge, KDim], np.int32]
+) -> common.Field[[Vertex, KDim], np.int32]:
+    return neighbor_sum(edge_f(V2E), axis=V2EDim)
+
+
+@gtx.field_operator
+def reduction_ke_field(
+    edge_f: common.Field[[KDim, Edge], np.int32]
+) -> common.Field[[KDim, Vertex], np.int32]:
+    return neighbor_sum(edge_f(V2E), axis=V2EDim)
+
+
+@pytest.mark.uses_unstructured_shift
+@pytest.mark.parametrize(
+    "fop", [reduction_e_field, reduction_ek_field, reduction_ke_field], ids=lambda fop: fop.__name__
+)
+def test_neighbor_sum(unstructured_case, fop):
     v2e_table = unstructured_case.offset_provider["V2E"].table
-    cases.verify_with_default_data(
+
+    edge_f = cases.allocate(unstructured_case, fop, "edge_f")()
+
+    local_dim_idx = edge_f.domain.dims.index(Edge) + 1
+    adv_indexing = tuple(
+        slice(None) if dim is not Edge else v2e_table for dim in edge_f.domain.dims
+    )
+
+    broadcast_slice = []
+    for dim in edge_f.domain.dims:
+        if dim is Edge:
+            broadcast_slice.append(slice(None))
+            broadcast_slice.append(slice(None))
+        else:
+            broadcast_slice.append(None)
+
+    broadcasted_table = v2e_table[tuple(broadcast_slice)]
+    ref = np.sum(
+        edge_f.asnumpy()[adv_indexing],
+        axis=local_dim_idx,
+        initial=0,
+        where=broadcasted_table != common._DEFAULT_SKIP_VALUE,
+    )
+    cases.verify(
         unstructured_case,
-        fencil,
-        ref=lambda edge_f: np.sum(
-            edge_f[v2e_table],
-            axis=1,
-            initial=0,
-            where=v2e_table != common.SKIP_VALUE,
-        ),
+        fop,
+        edge_f,
+        out=cases.allocate(unstructured_case, fop, cases.RETURN)(),
+        ref=ref,
     )
 
 
@@ -143,7 +177,7 @@ def test_reduction_execution_with_offset(unstructured_case):
             field.asnumpy()[:, 1][v2e_table],
             axis=1,
             initial=0,
-            where=v2e_table != common.SKIP_VALUE,
+            where=v2e_table != common._DEFAULT_SKIP_VALUE,
         ).reshape(out.shape),
         offset_provider=unstructured_case.offset_provider | {"Koff": KDim},
     )
@@ -171,7 +205,7 @@ def test_reduction_expression_in_call(unstructured_case):
             -edge_f[v2e_table] ** 2 * 2,
             axis=1,
             initial=0,
-            where=v2e_table != common.SKIP_VALUE,
+            where=v2e_table != common._DEFAULT_SKIP_VALUE,
         ),
     )
 
@@ -190,7 +224,7 @@ def test_reduction_with_common_expression(unstructured_case):
             flux[v2e_table] * 2,
             axis=1,
             initial=0,
-            where=v2e_table != common.SKIP_VALUE,
+            where=v2e_table != common._DEFAULT_SKIP_VALUE,
         ),
     )
 
diff --git a/tests/next_tests/integration_tests/multi_feature_tests/ffront_tests/test_ffront_fvm_nabla.py b/tests/next_tests/integration_tests/multi_feature_tests/ffront_tests/test_ffront_fvm_nabla.py
new file mode 100644
index 00000000..aeed607b
--- /dev/null
+++ b/tests/next_tests/integration_tests/multi_feature_tests/ffront_tests/test_ffront_fvm_nabla.py
@@ -0,0 +1,124 @@
+# 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 typing import Tuple
+
+import numpy as np
+import pytest
+
+
+pytest.importorskip("atlas4py")  # isort: skip
+
+from gt4py import next as gtx
+from gt4py.next import allocators, neighbor_sum
+from gt4py.next.iterator import atlas_utils
+from gt4py.next.program_processors import processor_interface as ppi
+
+from next_tests.integration_tests.feature_tests.ffront_tests.ffront_test_utils import (
+    exec_alloc_descriptor,
+)
+from next_tests.integration_tests.multi_feature_tests.fvm_nabla_setup import (
+    assert_close,
+    nabla_setup,
+)
+
+
+Vertex = gtx.Dimension("Vertex")
+Edge = gtx.Dimension("Edge")
+V2EDim = gtx.Dimension("V2E", kind=gtx.DimensionKind.LOCAL)
+E2VDim = gtx.Dimension("E2V", kind=gtx.DimensionKind.LOCAL)
+
+V2E = gtx.FieldOffset("V2E", source=Edge, target=(Vertex, V2EDim))
+E2V = gtx.FieldOffset("E2V", source=Vertex, target=(Edge, E2VDim))
+
+
+@gtx.field_operator
+def compute_zavgS(
+    pp: gtx.Field[[Vertex], float], S_M: gtx.Field[[Edge], float]
+) -> gtx.Field[[Edge], float]:
+    zavg = 0.5 * (pp(E2V[0]) + pp(E2V[1]))
+    return S_M * zavg
+
+
+@gtx.field_operator
+def compute_pnabla(
+    pp: gtx.Field[[Vertex], float],
+    S_M: gtx.Field[[Edge], float],
+    sign: gtx.Field[[Vertex, V2EDim], float],
+    vol: gtx.Field[[Vertex], float],
+) -> gtx.Field[[Vertex], float]:
+    zavgS = compute_zavgS(pp, S_M)
+    pnabla_M = neighbor_sum(zavgS(V2E) * sign, axis=V2EDim)
+    return pnabla_M / vol
+
+
+@gtx.field_operator
+def pnabla(
+    pp: gtx.Field[[Vertex], float],
+    S_M: Tuple[gtx.Field[[Edge], float], gtx.Field[[Edge], float]],
+    sign: gtx.Field[[Vertex, V2EDim], float],
+    vol: gtx.Field[[Vertex], float],
+) -> Tuple[gtx.Field[[Vertex], float], gtx.Field[[Vertex], float]]:
+    return compute_pnabla(pp, S_M[0], sign, vol), compute_pnabla(pp, S_M[1], sign, vol)
+
+
+def test_ffront_compute_zavgS(exec_alloc_descriptor):
+    executor, allocator = exec_alloc_descriptor.executor, exec_alloc_descriptor.allocator
+
+    setup = nabla_setup()
+
+    pp = gtx.as_field([Vertex], setup.input_field, allocator=allocator)
+    S_M = tuple(map(gtx.as_field.partial([Edge], allocator=allocator), setup.S_fields))
+
+    zavgS = gtx.zeros({Edge: setup.edges_size}, allocator=allocator)
+
+    e2v = gtx.NeighborTableOffsetProvider(
+        atlas_utils.AtlasTable(setup.edges2node_connectivity).asnumpy(), Edge, Vertex, 2, False
+    )
+
+    compute_zavgS.with_backend(executor)(pp, S_M[0], out=zavgS, offset_provider={"E2V": e2v})
+
+    assert_close(-199755464.25741270, np.min(zavgS.asnumpy()))
+    assert_close(388241977.58389181, np.max(zavgS.asnumpy()))
+
+
+def test_ffront_nabla(exec_alloc_descriptor):
+    executor, allocator = exec_alloc_descriptor.executor, exec_alloc_descriptor.allocator
+
+    setup = nabla_setup()
+
+    sign = gtx.as_field([Vertex, V2EDim], setup.sign_field, allocator=allocator)
+    pp = gtx.as_field([Vertex], setup.input_field, allocator=allocator)
+    S_M = tuple(map(gtx.as_field.partial([Edge], allocator=allocator), setup.S_fields))
+    vol = gtx.as_field([Vertex], setup.vol_field, allocator=allocator)
+
+    pnabla_MXX = gtx.zeros({Vertex: setup.nodes_size}, allocator=allocator)
+    pnabla_MYY = gtx.zeros({Vertex: setup.nodes_size}, allocator=allocator)
+
+    e2v = gtx.NeighborTableOffsetProvider(
+        atlas_utils.AtlasTable(setup.edges2node_connectivity).asnumpy(), Edge, Vertex, 2, False
+    )
+    v2e = gtx.NeighborTableOffsetProvider(
+        atlas_utils.AtlasTable(setup.nodes2edge_connectivity).asnumpy(), Vertex, Edge, 7
+    )
+
+    pnabla.with_backend(executor)(
+        pp, S_M, sign, vol, out=(pnabla_MXX, pnabla_MYY), offset_provider={"E2V": e2v, "V2E": v2e}
+    )
+
+    # TODO this check is not sensitive enough, need to implement a proper numpy reference!
+    assert_close(-3.5455427772566003e-003, np.min(pnabla_MXX.asnumpy()))
+    assert_close(3.5455427772565435e-003, np.max(pnabla_MXX.asnumpy()))
+    assert_close(-3.3540113705465301e-003, np.min(pnabla_MYY.asnumpy()))
+    assert_close(3.3540113705465301e-003, np.max(pnabla_MYY.asnumpy()))
diff --git a/tests/next_tests/integration_tests/multi_feature_tests/iterator_tests/fvm_nabla_setup.py b/tests/next_tests/integration_tests/multi_feature_tests/fvm_nabla_setup.py
similarity index 100%
rename from tests/next_tests/integration_tests/multi_feature_tests/iterator_tests/fvm_nabla_setup.py
rename to tests/next_tests/integration_tests/multi_feature_tests/fvm_nabla_setup.py
diff --git a/tests/next_tests/integration_tests/multi_feature_tests/iterator_tests/test_fvm_nabla.py b/tests/next_tests/integration_tests/multi_feature_tests/iterator_tests/test_fvm_nabla.py
index e1d959ab..9ee364c0 100644
--- a/tests/next_tests/integration_tests/multi_feature_tests/iterator_tests/test_fvm_nabla.py
+++ b/tests/next_tests/integration_tests/multi_feature_tests/iterator_tests/test_fvm_nabla.py
@@ -38,7 +38,7 @@ from gt4py.next.iterator.builtins import (
 from gt4py.next.iterator.runtime import closure, fendef, fundef, offset
 from gt4py.next.iterator.transforms.pass_manager import LiftMode
 
-from next_tests.integration_tests.multi_feature_tests.iterator_tests.fvm_nabla_setup import (
+from next_tests.integration_tests.multi_feature_tests.fvm_nabla_setup import (
     assert_close,
     nabla_setup,
 )
diff --git a/tests/next_tests/unit_tests/embedded_tests/test_nd_array_field.py b/tests/next_tests/unit_tests/embedded_tests/test_nd_array_field.py
index 79830a75..66189cf9 100644
--- a/tests/next_tests/unit_tests/embedded_tests/test_nd_array_field.py
+++ b/tests/next_tests/unit_tests/embedded_tests/test_nd_array_field.py
@@ -669,60 +669,60 @@ def test_connectivity_field_inverse_image():
 
 def test_connectivity_field_inverse_image_2d_domain():
     V = Dimension("V")
-    E = Dimension("E")
-    E2V = Dimension("E2V")
+    C = Dimension("C")
+    C2V = Dimension("C2V")
 
     V_START, V_STOP = 0, 3
-    E_START, E_STOP = 0, 3
-    E2V_START, E2V_STOP = 0, 3
+    C_START, C_STOP = 0, 3
+    C2V_START, C2V_STOP = 0, 3
 
-    e2v_conn = common._connectivity(
+    c2v_conn = common._connectivity(
         np.asarray([[0, 0, 2], [1, 1, 2], [2, 2, 2]]),
         domain=common.domain(
             [
-                common.named_range((E, (E_START, E_STOP))),
-                common.named_range((E2V, (E2V_START, E2V_STOP))),
+                common.named_range((C, (C_START, C_STOP))),
+                common.named_range((C2V, (C2V_START, C2V_STOP))),
             ]
         ),
         codomain=V,
     )
 
-    # e2c_conn:
-    #  ---E2V----
+    # c2v_conn:
+    #  ---C2V----
     #  |[[0 0 2]
-    #  E [1 1 2]
+    #  C [1 1 2]
     #  | [2 2 2]]
 
     # Test contiguous and non-contiguous ranges.
-    # For the 'e2c_conn' defined above, the only valid range including 2
+    # For the 'c2v_conn' defined above, the only valid range including 2
     # is [0, 3). Otherwise, the inverse image would be non-contiguous.
     image_range = UnitRange(V_START, V_STOP)
-    result = e2v_conn.inverse_image(image_range)
+    result = c2v_conn.inverse_image(image_range)
 
     assert len(result) == 2
-    assert result[0] == (E, UnitRange(E_START, E_STOP))
-    assert result[1] == (E2V, UnitRange(E2V_START, E2V_STOP))
+    assert result[0] == (C, UnitRange(C_START, C_STOP))
+    assert result[1] == (C2V, UnitRange(C2V_START, C2V_STOP))
 
-    result = e2v_conn.inverse_image(UnitRange(0, 2))
+    result = c2v_conn.inverse_image(UnitRange(0, 2))
     assert len(result) == 2
-    assert result[0] == (E, UnitRange(0, 2))
-    assert result[1] == (E2V, UnitRange(0, 2))
+    assert result[0] == (C, UnitRange(0, 2))
+    assert result[1] == (C2V, UnitRange(0, 2))
 
-    result = e2v_conn.inverse_image(UnitRange(0, 1))
+    result = c2v_conn.inverse_image(UnitRange(0, 1))
     assert len(result) == 2
-    assert result[0] == (E, UnitRange(0, 1))
-    assert result[1] == (E2V, UnitRange(0, 2))
+    assert result[0] == (C, UnitRange(0, 1))
+    assert result[1] == (C2V, UnitRange(0, 2))
 
-    result = e2v_conn.inverse_image(UnitRange(1, 2))
+    result = c2v_conn.inverse_image(UnitRange(1, 2))
     assert len(result) == 2
-    assert result[0] == (E, UnitRange(1, 2))
-    assert result[1] == (E2V, UnitRange(0, 2))
+    assert result[0] == (C, UnitRange(1, 2))
+    assert result[1] == (C2V, UnitRange(0, 2))
 
     with pytest.raises(ValueError, match="generates non-contiguous dimensions"):
-        result = e2v_conn.inverse_image(UnitRange(1, 3))
+        result = c2v_conn.inverse_image(UnitRange(1, 3))
 
     with pytest.raises(ValueError, match="generates non-contiguous dimensions"):
-        result = e2v_conn.inverse_image(UnitRange(2, 3))
+        result = c2v_conn.inverse_image(UnitRange(2, 3))
 
 
 def test_connectivity_field_inverse_image_non_contiguous():
@@ -746,3 +746,85 @@ def test_connectivity_field_inverse_image_non_contiguous():
 
     with pytest.raises(ValueError, match="generates non-contiguous dimensions"):
         e2v_conn.inverse_image(UnitRange(V_START, V_STOP))
+
+
+def test_connectivity_field_inverse_image_2d_domain_skip_values():
+    V = Dimension("V")
+    C = Dimension("C")
+    C2V = Dimension("C2V")
+
+    V_START, V_STOP = 0, 3
+    C_START, C_STOP = 0, 4
+    C2V_START, C2V_STOP = 0, 4
+
+    c2v_conn = common._connectivity(
+        np.asarray([[-1, 0, 2, -1], [1, 1, 2, 2], [2, 2, -1, -1], [-1, 2, -1, -1]]),
+        domain=common.domain(
+            [
+                common.named_range((C, (C_START, C_STOP))),
+                common.named_range((C2V, (C2V_START, C2V_STOP))),
+            ]
+        ),
+        codomain=V,
+        skip_value=-1,
+    )
+
+    # c2v_conn:
+    #  ---C2V---------
+    #  |[[-1  0  2 -1]
+    #  C [ 1  1  2  2]
+    #  | [ 2  2 -1 -1]
+    #  | [-1  2 -1 -1]]
+
+    image_range = UnitRange(V_START, V_STOP)
+    result = c2v_conn.inverse_image(image_range)
+
+    assert len(result) == 2
+    assert result[0] == (C, UnitRange(C_START, C_STOP))
+    assert result[1] == (C2V, UnitRange(C2V_START, C2V_STOP))
+
+    result = c2v_conn.inverse_image(UnitRange(0, 2))
+    assert len(result) == 2
+    assert result[0] == (C, UnitRange(0, 2))
+    assert result[1] == (C2V, UnitRange(0, 2))
+
+    result = c2v_conn.inverse_image(UnitRange(0, 1))
+    assert len(result) == 2
+    assert result[0] == (C, UnitRange(0, 1))
+    assert result[1] == (C2V, UnitRange(1, 2))
+
+    result = c2v_conn.inverse_image(UnitRange(1, 2))
+    assert len(result) == 2
+    assert result[0] == (C, UnitRange(1, 2))
+    assert result[1] == (C2V, UnitRange(0, 2))
+
+    with pytest.raises(ValueError, match="generates non-contiguous dimensions"):
+        result = c2v_conn.inverse_image(UnitRange(1, 3))
+
+    with pytest.raises(ValueError, match="generates non-contiguous dimensions"):
+        result = c2v_conn.inverse_image(UnitRange(2, 3))
+
+
+@pytest.mark.parametrize(
+    "index_array, expected",
+    [
+        ([0, 0, 1], [(0, 2)]),
+        ([0, 1, 0], None),
+        ([0, -1, 0], [(0, 3)]),
+        ([[1, 1, 1], [1, 0, 0]], [(1, 2), (1, 3)]),
+        (
+            [[1, 0, -1], [1, 0, 0]],
+            [(0, 2), (1, 3)],
+        ),
+    ],
+)
+def test_hypercube(index_array, expected):
+    index_array = np.asarray(index_array)
+    image_range = common.UnitRange(0, 1)
+    skip_value = -1
+
+    expected = [common.unit_range(e) for e in expected] if expected is not None else None
+
+    result = nd_array_field._hypercube(index_array, image_range, np, skip_value)
+
+    assert result == expected