Repository: GridTools/gt4py
Pull request #1296 (2023-07-24 13:20:24)

feat[next]: Basic support for field operations (same domain, no shifting)
Adds basic support for embedded field view:
- `_BaseNdArrayField` implements Fields with NumPy, CuPy and JAX buffers
- `fbuiltins.BuiltInFunction` gets support for dispatching to the implementation of a Field implementation
- `decorator.Program` and `decorator.FieldOperator` get support for executing the wrapped function directly

Support is limited to operations that don't change the domain, e.g. binary operations are only supported between fields with the same domain, shift/remap is not supported.

havogt: > One of the things I dislike about my CartesianSet implementation is that it does not give you a straightforward way to name the dimensions.

@tehrengruber One of the things we are still discussing is to make the `tuple[UnitRange]` a CartesianSet, because then we can move some generic functionality from the very specific Domain to CartesianSet (e.g. intersect). But we decided to go for the minimal implementation here and then see what other applications we have.namely 

Patch:
diff --git a/src/gt4py/_core/__init__.py b/src/gt4py/_core/__init__.py
index 6c43e2f1..373401e5 100644
--- a/src/gt4py/_core/__init__.py
+++ b/src/gt4py/_core/__init__.py
@@ -11,3 +11,5 @@
 # distribution for a copy of the license or check <https://www.gnu.org/licenses/>.
 #
 # SPDX-License-Identifier: GPL-3.0-or-later
+
+"""GridTools core definitions and utilities."""
diff --git a/src/gt4py/_core/definitions.py b/src/gt4py/_core/definitions.py
index d6b34c42..2546ae3e 100644
--- a/src/gt4py/_core/definitions.py
+++ b/src/gt4py/_core/definitions.py
@@ -101,6 +101,10 @@ ScalarT = TypeVar("ScalarT", bound=Union[BoolScalar, IntegralScalar, FloatingSca
 SCALAR_TYPES: Final[tuple[type, ...]] = (*BOOL_TYPES, *INTEGRAL_TYPES, *FLOAT_TYPES)
 
 
+def is_scalar_type(t: Any) -> TypeGuard[Scalar]:
+    return isinstance(t, SCALAR_TYPES)
+
+
 class BooleanIntegral(numbers.Integral):
     """Abstract base class for boolean integral types."""
 
@@ -390,46 +394,46 @@ class NDArrayObjectProto(Protocol):
     def dtype(self) -> Any:
         ...
 
-    def __getitem__(self, item: SliceLike) -> NDArrayObjectProto:
+    def __getitem__(self, item: SliceLike) -> NDArrayObject:
         ...
 
-    def __abs__(self) -> NDArrayObjectProto:
+    def __abs__(self) -> NDArrayObject:
         ...
 
-    def __neg__(self) -> NDArrayObjectProto:
+    def __neg__(self) -> NDArrayObject:
         ...
 
-    def __add__(self, other: NDArrayObject | Scalar) -> NDArrayObjectProto:
+    def __add__(self, other: NDArrayObject | Scalar) -> NDArrayObject:
         ...
 
-    def __radd__(self, other: NDArrayObject | Scalar) -> NDArrayObjectProto:
+    def __radd__(self, other: NDArrayObject | Scalar) -> NDArrayObject:
         ...
 
-    def __sub__(self, other: NDArrayObject | Scalar) -> NDArrayObjectProto:
+    def __sub__(self, other: NDArrayObject | Scalar) -> NDArrayObject:
         ...
 
-    def __rsub__(self, other: NDArrayObject | Scalar) -> NDArrayObjectProto:
+    def __rsub__(self, other: NDArrayObject | Scalar) -> NDArrayObject:
         ...
 
-    def __mul__(self, other: NDArrayObject | Scalar) -> NDArrayObjectProto:
+    def __mul__(self, other: NDArrayObject | Scalar) -> NDArrayObject:
         ...
 
-    def __rmul__(self, other: NDArrayObject | Scalar) -> NDArrayObjectProto:
+    def __rmul__(self, other: NDArrayObject | Scalar) -> NDArrayObject:
         ...
 
-    def __floordiv__(self, other: NDArrayObject | Scalar) -> NDArrayObjectProto:
+    def __floordiv__(self, other: NDArrayObject | Scalar) -> NDArrayObject:
         ...
 
-    def __rfloordiv__(self, other: NDArrayObject | Scalar) -> NDArrayObjectProto:
+    def __rfloordiv__(self, other: NDArrayObject | Scalar) -> NDArrayObject:
         ...
 
-    def __truediv__(self, other: NDArrayObject | Scalar) -> NDArrayObjectProto:
+    def __truediv__(self, other: NDArrayObject | Scalar) -> NDArrayObject:
         ...
 
-    def __rtruediv__(self, other: NDArrayObject | Scalar) -> NDArrayObjectProto:
+    def __rtruediv__(self, other: NDArrayObject | Scalar) -> NDArrayObject:
         ...
 
-    def __pow__(self, other: NDArrayObject | Scalar) -> NDArrayObjectProto:
+    def __pow__(self, other: NDArrayObject | Scalar) -> NDArrayObject:
         ...
 
 
diff --git a/src/gt4py/eve/extended_typing.py b/src/gt4py/eve/extended_typing.py
index 5daa7e1f..34829317 100644
--- a/src/gt4py/eve/extended_typing.py
+++ b/src/gt4py/eve/extended_typing.py
@@ -301,6 +301,10 @@ class NonProtocolABCMeta(_typing._ProtocolMeta):
     __subclasshook__ = None  # type: ignore[assignment]
 
 
+class NonProtocolABC(metaclass=NonProtocolABCMeta):
+    pass
+
+
 _ProtoT = TypeVar("_ProtoT", bound=_abc.ABCMeta)
 
 
diff --git a/src/gt4py/next/common.py b/src/gt4py/next/common.py
index b0f0b8ac..0e7f7aef 100644
--- a/src/gt4py/next/common.py
+++ b/src/gt4py/next/common.py
@@ -14,32 +14,48 @@
 
 from __future__ import annotations
 
+import abc
+import dataclasses
 import enum
-from collections.abc import Sequence
-from dataclasses import dataclass
-from typing import (
+import functools
+import sys
+from collections.abc import Sequence, Set
+from typing import overload
+
+import numpy as np
+import numpy.typing as npt
+
+from gt4py._core import definitions as core_defs
+from gt4py.eve.extended_typing import (
+    TYPE_CHECKING,
     Any,
-    Generic,
+    Callable,
+    ClassVar,
     Optional,
+    ParamSpec,
     Protocol,
-    SupportsFloat,
-    SupportsInt,
     TypeAlias,
     TypeVar,
+    Union,
+    extended_runtime_checkable,
+    final,
     runtime_checkable,
 )
-
-import numpy as np
-import numpy.typing as npt
-
 from gt4py.eve.type_definitions import StrEnum
 
 
 DimT = TypeVar("DimT", bound="Dimension")
-DimsT = TypeVar("DimsT", bound=Sequence["Dimension"])
-DT = TypeVar("DT", bound="Scalar")
+DimsT = TypeVar("DimsT", bound=Sequence["Dimension"], covariant=True)
+
+
+class Infinity(int):
+    @classmethod
+    def positive(cls) -> "Infinity":
+        return cls(sys.maxsize)
 
-Scalar: TypeAlias = SupportsInt | SupportsFloat | np.int32 | np.int64 | np.float32 | np.float64
+    @classmethod
+    def negative(cls) -> "Infinity":
+        return cls(-sys.maxsize)
 
 
 @enum.unique
@@ -52,7 +68,7 @@ class DimensionKind(StrEnum):
         return f"{type(self).__name__}.{self.name}"
 
 
-@dataclass(frozen=True)
+@dataclasses.dataclass(frozen=True)
 class Dimension:
     value: str
     kind: DimensionKind = DimensionKind.HORIZONTAL
@@ -61,21 +77,266 @@ class Dimension:
         return f'Dimension(value="{self.value}", kind={self.kind})'
 
 
-class DType:
-    ...
+DomainLike: TypeAlias = Union[
+    Sequence[Dimension], Dimension, str
+]  # TODO(havogt): revisit once embedded implementation is concluded
+
+
+@dataclasses.dataclass(frozen=True)
+class UnitRange(Sequence[int], Set[int]):
+    """Range from `start` to `stop` with step size one."""
+
+    start: int
+    stop: int
+
+    def __post_init__(self):
+        if self.stop <= self.start:
+            # make UnitRange(0,0) the single empty UnitRange
+            object.__setattr__(self, "start", 0)
+            object.__setattr__(self, "stop", 0)
+
+    def __len__(self) -> int:
+        if Infinity.positive() in (abs(self.start), abs(self.stop)):
+            return Infinity.positive()
+        return max(0, self.stop - self.start)
+
+    def __repr__(self) -> str:
+        return f"UnitRange({self.start}, {self.stop})"
+
+    @overload
+    def __getitem__(self, index: int) -> int:
+        ...
+
+    @overload
+    def __getitem__(self, index: slice) -> UnitRange:
+        ...
+
+    def __getitem__(self, index: int | slice) -> int | UnitRange:
+        if isinstance(index, slice):
+            start, stop, step = index.indices(len(self))
+            if step != 1:
+                raise ValueError("UnitRange: step required to be `1`.")
+            new_start = self.start + (start or 0)
+            new_stop = (self.start if stop > 0 else self.stop) + stop
+            return UnitRange(new_start, new_stop)
+        else:
+            if index < 0:
+                index += len(self)
+
+            if 0 <= index < len(self):
+                return self.start + index
+            else:
+                raise IndexError("UnitRange index out of range")
+
+    def __and__(self, other: Set[Any]) -> UnitRange:
+        if isinstance(other, UnitRange):
+            start = max(self.start, other.start)
+            stop = min(self.stop, other.stop)
+            return UnitRange(start, stop)
+        else:
+            raise NotImplementedError("Can only find the intersection between UnitRange instances.")
+
+
+NamedRange: TypeAlias = tuple[Dimension, UnitRange]
+
+
+@dataclasses.dataclass(frozen=True)
+class Domain(Sequence[NamedRange]):
+    dims: tuple[Dimension, ...]
+    ranges: tuple[UnitRange, ...]
+
+    def __post_init__(self):
+        if len(set(self.dims)) != len(self.dims):
+            raise NotImplementedError(f"Domain dimensions must be unique, not {self.dims}.")
+
+    def __len__(self) -> int:
+        return len(self.ranges)
+
+    @overload
+    def __getitem__(self, index: int) -> NamedRange:
+        ...
+
+    @overload
+    def __getitem__(self, index: slice) -> "Domain":
+        ...
+
+    @overload
+    def __getitem__(self, index: Dimension) -> NamedRange:
+        ...
+
+    def __getitem__(self, index: int | slice | Dimension) -> NamedRange | Domain:
+        if isinstance(index, int):
+            return self.dims[index], self.ranges[index]
+        elif isinstance(index, slice):
+            dims_slice = self.dims[index]
+            ranges_slice = self.ranges[index]
+            return Domain(dims_slice, ranges_slice)
+        elif isinstance(index, Dimension):
+            try:
+                index_pos = self.dims.index(index)
+                return self.dims[index_pos], self.ranges[index_pos]
+            except ValueError:
+                raise KeyError(f"No Dimension of type {index} is present in the Domain.")
+        else:
+            raise KeyError("Invalid index type, must be either int, slice, or Dimension.")
+
+    def __and__(self, other: "Domain") -> "Domain":
+        broadcast_dims = tuple(promote_dims(self.dims, other.dims))
+        intersected_ranges = tuple(
+            rng1 & rng2
+            for rng1, rng2 in zip(
+                _broadcast_ranges(broadcast_dims, self.dims, self.ranges),
+                _broadcast_ranges(broadcast_dims, other.dims, other.ranges),
+            )
+        )
+        return Domain(broadcast_dims, intersected_ranges)
+
+
+def _broadcast_ranges(
+    broadcast_dims: Sequence[Dimension], dims: Sequence[Dimension], ranges: Sequence[UnitRange]
+) -> tuple[UnitRange, ...]:
+    return tuple(
+        ranges[dims.index(d)] if d in dims else UnitRange(Infinity.negative(), Infinity.positive())
+        for d in broadcast_dims
+    )
+
+
+if TYPE_CHECKING:
+    import gt4py.next.ffront.fbuiltins as fbuiltins
+
+    _Value: TypeAlias = "Field" | core_defs.ScalarT
+    _P = ParamSpec("_P")
+    _R = TypeVar("_R", _Value, tuple[_Value, ...])
+
+    class GTBuiltInFuncDispatcher(Protocol):
+        def __call__(self, func: fbuiltins.BuiltInFunction[_R, _P], /) -> Callable[_P, _R]:
+            ...
 
 
-class Field(Generic[DimsT, DT]):
-    ...
+@extended_runtime_checkable
+class Field(Protocol[DimsT, core_defs.ScalarT]):
+    __gt_builtin_func__: ClassVar[GTBuiltInFuncDispatcher]
 
+    @property
+    def domain(self) -> Domain:
+        ...
 
-@dataclass(frozen=True)
+    @property
+    def dtype(self) -> core_defs.DType[core_defs.ScalarT]:
+        ...
+
+    @property
+    def value_type(self) -> type[core_defs.ScalarT]:
+        ...
+
+    @property
+    def ndarray(self) -> core_defs.NDArrayObject:
+        ...
+
+    def __str__(self) -> str:
+        codomain = self.value_type.__name__
+        return f"⟨{self.domain!s} → {codomain}⟩"
+
+    @abc.abstractmethod
+    def remap(self, index_field: Field) -> Field:
+        ...
+
+    @abc.abstractmethod
+    def restrict(self, item: "DomainLike") -> Field:
+        ...
+
+    # Operators
+    @abc.abstractmethod
+    def __call__(self, index_field: Field) -> Field:
+        ...
+
+    @abc.abstractmethod
+    def __getitem__(self, item: "DomainLike") -> Field:
+        ...
+
+    @abc.abstractmethod
+    def __abs__(self) -> Field:
+        ...
+
+    @abc.abstractmethod
+    def __neg__(self) -> Field:
+        ...
+
+    @abc.abstractmethod
+    def __add__(self, other: Field | core_defs.ScalarT) -> Field:
+        ...
+
+    @abc.abstractmethod
+    def __radd__(self, other: Field | core_defs.ScalarT) -> Field:
+        ...
+
+    @abc.abstractmethod
+    def __sub__(self, other: Field | core_defs.ScalarT) -> Field:
+        ...
+
+    @abc.abstractmethod
+    def __rsub__(self, other: Field | core_defs.ScalarT) -> Field:
+        ...
+
+    @abc.abstractmethod
+    def __mul__(self, other: Field | core_defs.ScalarT) -> Field:
+        ...
+
+    @abc.abstractmethod
+    def __rmul__(self, other: Field | core_defs.ScalarT) -> Field:
+        ...
+
+    @abc.abstractmethod
+    def __floordiv__(self, other: Field | core_defs.ScalarT) -> Field:
+        ...
+
+    @abc.abstractmethod
+    def __rfloordiv__(self, other: Field | core_defs.ScalarT) -> Field:
+        ...
+
+    @abc.abstractmethod
+    def __truediv__(self, other: Field | core_defs.ScalarT) -> Field:
+        ...
+
+    @abc.abstractmethod
+    def __rtruediv__(self, other: Field | core_defs.ScalarT) -> Field:
+        ...
+
+    @abc.abstractmethod
+    def __pow__(self, other: Field | core_defs.ScalarT) -> Field:
+        ...
+
+
+class FieldABC(Field[DimsT, core_defs.ScalarT]):
+    """Abstract base class for implementations of the :class:`Field` protocol."""
+
+    @final
+    def __setattr__(self, key, value) -> None:
+        raise TypeError("Immutable type")
+
+    @final
+    def __setitem__(self, key, value) -> None:
+        raise TypeError("Immutable type")
+
+
+@functools.singledispatch
+def field(
+    definition: Any,
+    /,
+    *,
+    domain: Optional[Any] = None,  # TODO(havogt): provide domain_like to Domain conversion
+    value_type: Optional[type] = None,
+) -> Field:
+    raise NotImplementedError
+
+
+@dataclasses.dataclass(frozen=True)
 class GTInfo:
     definition: Any
     ir: Any
 
 
-@dataclass(frozen=True)
+@dataclasses.dataclass(frozen=True)
 class Backend:
     def __init__(self, **kwargs: Any) -> None:
         super().__init__(**kwargs)
@@ -108,3 +369,79 @@ class NeighborTable(Connectivity, Protocol):
 class GridType(StrEnum):
     CARTESIAN = "cartesian"
     UNSTRUCTURED = "unstructured"
+
+
+def promote_dims(*dims_list: Sequence[Dimension]) -> list[Dimension]:
+    """
+    Find a unique ordering of multiple (individually ordered) lists of dimensions.
+
+    The resulting list of dimensions contains all dimensions of the arguments
+    in the order they originally appear. If no unique order exists or a
+    contradicting order is found an exception is raised.
+
+    A modified version (ensuring uniqueness of the order) of
+    `Kahn's algorithm <https://en.wikipedia.org/wiki/Topological_sorting#Kahn's_algorithm>`_
+    is used to topologically sort the arguments.
+    >>> from gt4py.next.common import Dimension
+    >>> I, J, K = (Dimension(value=dim) for dim in ["I", "J", "K"])
+    >>> promote_dims([I, J], [I, J, K]) == [I, J, K]
+    True
+    >>> promote_dims([I, J], [K]) # doctest: +ELLIPSIS
+    Traceback (most recent call last):
+     ...
+    ValueError: Dimensions can not be promoted. Could not determine order of the following dimensions: J, K.
+    >>> promote_dims([I, J], [J, I]) # doctest: +ELLIPSIS
+    Traceback (most recent call last):
+     ...
+    ValueError: Dimensions can not be promoted. The following dimensions appear in contradicting order: I, J.
+    """
+    # build a graph with the vertices being dimensions and edges representing
+    #  the order between two dimensions. The graph is encoded as a dictionary
+    #  mapping dimensions to their predecessors, i.e. a dictionary containing
+    #  adjacency lists. Since graphlib.TopologicalSorter uses predecessors
+    #  (contrary to successors) we also use this directionality here.
+    graph: dict[Dimension, set[Dimension]] = {}
+    for dims in dims_list:
+        if len(dims) == 0:
+            continue
+        # create a vertex for each dimension
+        for dim in dims:
+            graph.setdefault(dim, set())
+        # add edges
+        predecessor = dims[0]
+        for dim in dims[1:]:
+            graph[dim].add(predecessor)
+            predecessor = dim
+
+    # modified version of Kahn's algorithm
+    topologically_sorted_list: list[Dimension] = []
+
+    # compute in-degree for each vertex
+    in_degree = {v: 0 for v in graph.keys()}
+    for v1 in graph:
+        for v2 in graph[v1]:
+            in_degree[v2] += 1
+
+    # process vertices with in-degree == 0
+    # TODO(tehrengruber): avoid recomputation of zero_in_degree_vertex_list
+    while zero_in_degree_vertex_list := [v for v, d in in_degree.items() if d == 0]:
+        if len(zero_in_degree_vertex_list) != 1:
+            raise ValueError(
+                f"Dimensions can not be promoted. Could not determine "
+                f"order of the following dimensions: "
+                f"{', '.join((dim.value for dim in zero_in_degree_vertex_list))}."
+            )
+        v = zero_in_degree_vertex_list[0]
+        del in_degree[v]
+        topologically_sorted_list.insert(0, v)
+        # update in-degree
+        for predecessor in graph[v]:
+            in_degree[predecessor] -= 1
+
+    if len(in_degree.items()) > 0:
+        raise ValueError(
+            f"Dimensions can not be promoted. The following dimensions "
+            f"appear in contradicting order: {', '.join((dim.value for dim in in_degree.keys()))}."
+        )
+
+    return topologically_sorted_list
diff --git a/src/gt4py/next/embedded/__init__.py b/src/gt4py/next/embedded/__init__.py
new file mode 100644
index 00000000..6c43e2f1
--- /dev/null
+++ b/src/gt4py/next/embedded/__init__.py
@@ -0,0 +1,13 @@
+# 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
diff --git a/src/gt4py/next/embedded/nd_array_field.py b/src/gt4py/next/embedded/nd_array_field.py
new file mode 100644
index 00000000..66c09795
--- /dev/null
+++ b/src/gt4py/next/embedded/nd_array_field.py
@@ -0,0 +1,256 @@
+# 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
+from collections.abc import Callable
+from types import ModuleType
+from typing import ClassVar, Optional, ParamSpec, TypeAlias, TypeVar, overload
+
+import numpy as np
+from numpy import typing as npt
+
+from gt4py.next import common
+
+
+try:
+    import cupy as cp
+except ImportError:
+    cp: Optional[ModuleType] = None  # type:ignore[no-redef]
+
+try:
+    from jax import numpy as jnp
+except ImportError:
+    jnp: Optional[ModuleType] = None  # type:ignore[no-redef]
+
+
+from gt4py._core import definitions
+from gt4py._core.definitions import ScalarT
+from gt4py.next.common import DimsT, Domain
+from gt4py.next.ffront import fbuiltins
+
+
+def _make_unary_array_field_intrinsic_func(builtin_name: str, array_builtin_name: str) -> Callable:
+    def _builtin_unary_op(a: _BaseNdArrayField) -> common.Field:
+        xp = a.__class__.array_ns
+        op = getattr(xp, array_builtin_name)
+        new_data = op(a.ndarray)
+
+        return a.__class__.from_array(new_data, domain=a.domain)
+
+    _builtin_unary_op.__name__ = builtin_name
+    return _builtin_unary_op
+
+
+def _make_binary_array_field_intrinsic_func(builtin_name: str, array_builtin_name: str) -> Callable:
+    def _builtin_binary_op(a: _BaseNdArrayField, b: common.Field) -> common.Field:
+        xp = a.__class__.array_ns
+        op = getattr(xp, array_builtin_name)
+        if hasattr(b, "__gt_builtin_func__"):  # isinstance(b, common.Field):
+            if not a.domain == b.domain:
+                raise NotImplementedError(
+                    f"support for different domain not implemented: {a.domain}, {b.domain}"
+                )
+            new_data = op(a.ndarray, xp.asarray(b.ndarray))
+        else:
+            assert isinstance(b, definitions.SCALAR_TYPES)
+            new_data = op(a.ndarray, b)
+
+        return a.__class__.from_array(new_data, domain=a.domain)
+
+    _builtin_binary_op.__name__ = builtin_name
+    return _builtin_binary_op
+
+
+_Value: TypeAlias = common.Field | ScalarT
+_P = ParamSpec("_P")
+_R = TypeVar("_R", _Value, tuple[_Value, ...])
+
+
+@dataclasses.dataclass(frozen=True)
+class _BaseNdArrayField(common.FieldABC[DimsT, ScalarT]):
+    """
+    Shared field implementation for NumPy-like fields.
+
+    Builtin function implementations are registered in a dictionary.
+    Note: Currently, all concrete NdArray-implementations share
+    the same implementation, dispatching is handled inside of the registered
+    function via its namespace.
+    """
+
+    _domain: Domain
+    _ndarray: definitions.NDArrayObject
+    _value_type: type[ScalarT]
+
+    array_ns: ClassVar[
+        ModuleType
+    ]  # TODO(havogt) after storage PR is merged, update to the NDArrayNamespace protocol
+
+    _builtin_func_map: ClassVar[dict[fbuiltins.BuiltInFunction, Callable]] = {}
+
+    @classmethod
+    def __gt_builtin_func__(cls, func: fbuiltins.BuiltInFunction[_R, _P], /) -> Callable[_P, _R]:
+        return cls._builtin_func_map.get(func, NotImplemented)
+
+    @overload
+    @classmethod
+    def register_builtin_func(
+        cls, op: fbuiltins.BuiltInFunction[_R, _P], op_func: None
+    ) -> functools.partial[Callable[_P, _R]]:
+        ...
+
+    @overload
+    @classmethod
+    def register_builtin_func(
+        cls, op: fbuiltins.BuiltInFunction[_R, _P], op_func: Callable[_P, _R]
+    ) -> Callable[_P, _R]:
+        ...
+
+    @classmethod
+    def register_builtin_func(
+        cls, op: fbuiltins.BuiltInFunction[_R, _P], op_func: Optional[Callable[_P, _R]] = None
+    ) -> Callable[_P, _R] | functools.partial[Callable[_P, _R]]:
+        assert op not in cls._builtin_func_map
+        if op_func is None:  # when used as a decorator
+            return functools.partial(cls.register_builtin_func, op)  # type: ignore[arg-type]
+        return cls._builtin_func_map.setdefault(op, op_func)
+
+    @property
+    def domain(self) -> Domain:
+        return self._domain
+
+    @property
+    def ndarray(self) -> definitions.NDArrayObject:
+        return self._ndarray
+
+    @property
+    def value_type(self) -> type[definitions.ScalarT]:
+        return self._value_type
+
+    @classmethod
+    def from_array(
+        cls,
+        data: npt.ArrayLike,
+        /,
+        *,
+        domain: Domain,
+        value_type: Optional[type] = None,
+    ) -> _BaseNdArrayField:
+        xp = cls.array_ns
+        dtype = None
+        if value_type is not None:
+            dtype = xp.dtype(value_type)
+        array = xp.asarray(data, dtype=dtype)
+
+        value_type = array.dtype.type  # TODO add support for Dimensions as value_type
+
+        assert issubclass(array.dtype.type, definitions.SCALAR_TYPES)
+
+        assert all(isinstance(d, common.Dimension) for d, r in domain), domain
+        assert len(domain) == array.ndim
+        assert all(len(nr[1]) == s for nr, s in zip(domain, array.shape))
+
+        assert value_type is not None  # for mypy
+        return cls(domain, array, value_type)
+
+    def remap(self: _BaseNdArrayField, connectivity) -> _BaseNdArrayField:
+        raise NotImplementedError()
+
+    def restrict(self: _BaseNdArrayField, domain) -> _BaseNdArrayField:
+        raise NotImplementedError()
+
+    __call__ = None  # type: ignore[assignment]  # TODO: remap
+
+    __getitem__ = None  # type: ignore[assignment]  # TODO: restrict
+
+    __abs__ = _make_unary_array_field_intrinsic_func("abs", "abs")
+
+    __neg__ = _make_unary_array_field_intrinsic_func("neg", "negative")
+
+    __add__ = __radd__ = _make_binary_array_field_intrinsic_func("add", "add")
+
+    __sub__ = __rsub__ = _make_binary_array_field_intrinsic_func("sub", "subtract")
+
+    __mul__ = __rmul__ = _make_binary_array_field_intrinsic_func("mul", "multiply")
+
+    __truediv__ = __rtruediv__ = _make_binary_array_field_intrinsic_func("div", "divide")
+
+    __floordiv__ = __rfloordiv__ = _make_binary_array_field_intrinsic_func(
+        "floordiv", "floor_divide"
+    )
+
+    __pow__ = _make_binary_array_field_intrinsic_func("pow", "power")
+
+
+# -- Specialized implementations for intrinsic operations on array fields --
+
+_BaseNdArrayField.register_builtin_func(fbuiltins.abs, _BaseNdArrayField.__abs__)  # type: ignore[attr-defined]
+_BaseNdArrayField.register_builtin_func(fbuiltins.power, _BaseNdArrayField.__pow__)  # type: ignore[attr-defined]
+# TODO gamma
+
+for name in (
+    fbuiltins.UNARY_MATH_FP_BUILTIN_NAMES
+    + fbuiltins.UNARY_MATH_FP_PREDICATE_BUILTIN_NAMES
+    + fbuiltins.UNARY_MATH_NUMBER_BUILTIN_NAMES
+):
+    if name in ["abs", "power", "gamma"]:
+        continue
+    _BaseNdArrayField.register_builtin_func(
+        getattr(fbuiltins, name), _make_unary_array_field_intrinsic_func(name, name)
+    )
+
+_BaseNdArrayField.register_builtin_func(
+    fbuiltins.minimum, _make_binary_array_field_intrinsic_func("minimum", "minimum")  # type: ignore[attr-defined]
+)
+_BaseNdArrayField.register_builtin_func(
+    fbuiltins.maximum, _make_binary_array_field_intrinsic_func("maximum", "maximum")  # type: ignore[attr-defined]
+)
+_BaseNdArrayField.register_builtin_func(
+    fbuiltins.fmod, _make_binary_array_field_intrinsic_func("fmod", "fmod")  # type: ignore[attr-defined]
+)
+
+# -- Concrete array implementations --
+# NumPy
+_nd_array_implementations = [np]
+
+
+@dataclasses.dataclass(frozen=True)
+class NumPyArrayField(_BaseNdArrayField):
+    array_ns: ClassVar[ModuleType] = np
+
+
+common.field.register(np.ndarray, NumPyArrayField.from_array)
+
+
+# CuPy
+if cp:
+    _nd_array_implementations.append(cp)
+
+    @dataclasses.dataclass(frozen=True)
+    class CuPyArrayField(_BaseNdArrayField):
+        array_ns: ClassVar[ModuleType] = cp
+
+    common.field.register(cp.ndarray, CuPyArrayField.from_array)
+
+# JAX
+if jnp:
+    _nd_array_implementations.append(jnp)
+
+    @dataclasses.dataclass(frozen=True)
+    class JaxArrayField(_BaseNdArrayField):
+        array_ns: ClassVar[ModuleType] = jnp
+
+    common.field.register(jnp.ndarray, JaxArrayField.from_array)
diff --git a/src/gt4py/next/ffront/decorator.py b/src/gt4py/next/ffront/decorator.py
index 2a343454..5b4d32b5 100644
--- a/src/gt4py/next/ffront/decorator.py
+++ b/src/gt4py/next/ffront/decorator.py
@@ -29,9 +29,10 @@ from typing import Generator, Generic, TypeVar
 
 from devtools import debug
 
+from gt4py._core import definitions as core_defs
 from gt4py.eve.extended_typing import Any, Optional
 from gt4py.eve.utils import UIDGenerator
-from gt4py.next.common import Dimension, DimensionKind, GridType, Scalar
+from gt4py.next.common import Dimension, DimensionKind, GridType
 from gt4py.next.ffront import (
     dialect_ast_enums,
     field_operator_ast as foast,
@@ -236,6 +237,12 @@ class Program:
         )
 
     def __call__(self, *args, offset_provider: dict[str, Dimension], **kwargs) -> None:
+        if (
+            self.backend is None and DEFAULT_BACKEND is None
+        ):  # TODO(havogt): for now enable embedded execution by setting DEFAULT_BACKEND to None
+            self.definition(*args, **kwargs)
+            return
+
         rewritten_args, size_args, kwargs = self._process_args(args, kwargs)
 
         if not self.backend:
@@ -552,23 +559,36 @@ class FieldOperator(GTCallable, Generic[OperatorNodeT]):
     def __call__(
         self,
         *args,
-        out,
-        offset_provider: dict[str, Dimension],
         **kwargs,
     ) -> None:
-        args, kwargs = type_info.canonicalize_arguments(self.foast_node.type, args, kwargs)
-        # TODO(tehrengruber): check all offset providers are given
-        # deduce argument types
-        arg_types = []
-        for arg in args:
-            arg_types.append(type_translation.from_value(arg))
-        kwarg_types = {}
-        for name, arg in kwargs.items():
-            kwarg_types[name] = type_translation.from_value(arg)
-
-        return self.as_program(arg_types, kwarg_types)(
-            *args, out, offset_provider=offset_provider, **kwargs
-        )
+        # TODO(havogt): Don't select mode based on existence of kwargs,
+        # because now we cannot provide nice error messages. E.g. set context var
+        # if we are reaching this from a program call.
+        if "out" in kwargs:
+            out = kwargs.pop("out")
+            if "offset_provider" in kwargs:
+                # "out" and "offset_provider" -> field_operator as program
+                offset_provider = kwargs.pop("offset_provider")
+                args, kwargs = type_info.canonicalize_arguments(self.foast_node.type, args, kwargs)
+                # TODO(tehrengruber): check all offset providers are given
+                # deduce argument types
+                arg_types = []
+                for arg in args:
+                    arg_types.append(type_translation.from_value(arg))
+                kwarg_types = {}
+                for name, arg in kwargs.items():
+                    kwarg_types[name] = type_translation.from_value(arg)
+
+                return self.as_program(arg_types, kwarg_types)(
+                    *args, out, offset_provider=offset_provider, **kwargs
+                )
+            else:
+                # "out" -> field_operator called from program in embedded execution
+                out.ndarray[:] = self.definition(*args, **kwargs).ndarray[:]
+                return
+        else:
+            # field_operator called from other field_operator in embedded execution
+            return self.definition(*args, **kwargs)
 
 
 @typing.overload
@@ -614,7 +634,7 @@ def scan_operator(
     *,
     axis: Dimension,
     forward: bool,
-    init: Scalar,
+    init: core_defs.Scalar,
     backend: Optional[str],
 ) -> FieldOperator[foast.ScanOperator]:
     ...
@@ -625,7 +645,7 @@ def scan_operator(
     *,
     axis: Dimension,
     forward: bool,
-    init: Scalar,
+    init: core_defs.Scalar,
     backend: Optional[str],
 ) -> Callable[[types.FunctionType], FieldOperator[foast.ScanOperator]]:
     ...
@@ -636,7 +656,7 @@ def scan_operator(
     *,
     axis: Dimension,
     forward: bool = True,
-    init: Scalar = 0.0,
+    init: core_defs.Scalar = 0.0,
     backend=None,
 ) -> (
     FieldOperator[foast.ScanOperator]
diff --git a/src/gt4py/next/ffront/fbuiltins.py b/src/gt4py/next/ffront/fbuiltins.py
index 52842336..2129bd6d 100644
--- a/src/gt4py/next/ffront/fbuiltins.py
+++ b/src/gt4py/next/ffront/fbuiltins.py
@@ -12,12 +12,25 @@
 #
 # SPDX-License-Identifier: GPL-3.0-or-later
 
+import dataclasses
+import inspect
 from builtins import bool, float, int, tuple
-from dataclasses import dataclass
-from typing import TypeAlias
+from typing import (
+    Any,
+    Callable,
+    Generic,
+    Optional,
+    ParamSpec,
+    Tuple,
+    TypeAlias,
+    TypeVar,
+    Union,
+    cast,
+)
 
 from numpy import float32, float64, int32, int64
 
+from gt4py._core import definitions as gt4py_defs
 from gt4py.next.common import Dimension, DimensionKind, Field
 from gt4py.next.ffront.experimental import as_offset  # noqa F401
 from gt4py.next.iterator import runtime
@@ -37,103 +50,138 @@ IndexType: TypeAlias = int32
 TYPE_ALIAS_NAMES = ["IndexType"]
 
 
-@dataclass
-class BuiltInFunction:
-    __gt_type: ts.FunctionType
+_P = ParamSpec("_P")
+_R = TypeVar("_R")
 
-    def __call__(self, *args, **kwargs):
-        """Act as an empty place holder for the built in function."""
 
-    def __gt_type__(self):
-        return self.__gt_type
+def _type_conversion_helper(t: type) -> type[ts.TypeSpec] | tuple[type[ts.TypeSpec], ...]:
+    if t is Field:
+        return ts.FieldType
+    elif t is Dimension:
+        return ts.DimensionType
+    elif t is gt4py_defs.ScalarT:
+        return ts.ScalarType
+    elif t is type:
+        return (
+            ts.FunctionType
+        )  # our type of type is currently represented by the type constructor function
+    elif t is Tuple:
+        return ts.TupleType
+    elif hasattr(t, "__origin__") and t.__origin__ is Union:
+        types = [_type_conversion_helper(e) for e in t.__args__]  # type: ignore[attr-defined]
+        assert all(type(t) is type and issubclass(t, ts.TypeSpec) for t in types)
+        return cast(tuple[type[ts.TypeSpec], ...], tuple(types))  # `cast` to break the recursion
+    else:
+        raise AssertionError("Illegal type encountered.")
 
 
-_reduction_like = BuiltInFunction(
-    ts.FunctionType(
-        pos_only_args=[ts.DeferredType(constraint=ts.FieldType)],
-        pos_or_kw_args={"axis": ts.DeferredType(constraint=ts.DimensionType)},
-        kw_only_args={},
-        returns=ts.DeferredType(constraint=ts.FieldType),
-    )
-)
+def _type_conversion(t: type) -> ts.DeferredType:
+    return ts.DeferredType(constraint=_type_conversion_helper(t))
 
-neighbor_sum = _reduction_like
-max_over = _reduction_like
-min_over = _reduction_like
-
-broadcast = BuiltInFunction(
-    ts.FunctionType(
-        pos_only_args=[
-            ts.DeferredType(constraint=(ts.FieldType, ts.ScalarType)),
-            ts.DeferredType(constraint=ts.TupleType),
-        ],
-        pos_or_kw_args={},
-        kw_only_args={},
-        returns=ts.DeferredType(constraint=ts.FieldType),
-    )
-)
 
-where = BuiltInFunction(
-    ts.FunctionType(
-        pos_only_args=[
-            ts.DeferredType(constraint=ts.FieldType),
-            ts.DeferredType(constraint=(ts.FieldType, ts.ScalarType, ts.TupleType)),
-            ts.DeferredType(constraint=(ts.FieldType, ts.ScalarType, ts.TupleType)),
-        ],
-        pos_or_kw_args={},
-        kw_only_args={},
-        returns=ts.DeferredType(constraint=(ts.FieldType, ts.TupleType)),
-    )
-)
-
-astype = BuiltInFunction(
-    ts.FunctionType(
-        pos_only_args=[
-            ts.DeferredType(constraint=ts.FieldType),
-            ts.DeferredType(constraint=ts.FunctionType),
-        ],
-        pos_or_kw_args={},
-        kw_only_args={},
-        returns=ts.DeferredType(constraint=ts.FieldType),
-    )
-)
+@dataclasses.dataclass(frozen=True)
+class BuiltInFunction(Generic[_R, _P]):
+    name: str = dataclasses.field(init=False)
+    # `function` can be used to provide a default implementation for all `Field` implementations,
+    # e.g. a fused multiply add could have a default implementation as a*b+c, but an optimized implementation for a specific `Field`
+    function: Callable[_P, _R]
 
-_unary_math_builtin = BuiltInFunction(
-    ts.FunctionType(
-        pos_only_args=[ts.DeferredType(constraint=(ts.ScalarType, ts.FieldType))],
-        pos_or_kw_args={},
-        kw_only_args={},
-        returns=ts.DeferredType(constraint=(ts.ScalarType, ts.FieldType)),
-    )
-)
+    def __post_init__(self):
+        object.__setattr__(self, "name", f"{self.function.__module__}.{self.function.__name__}")
+
+    def __call__(self, *args: _P.args, **kwargs: _P.kwargs) -> _R:
+        impl = self.dispatch(*args)
+        return impl(*args, **kwargs)
+
+    def dispatch(self, *args: Any) -> Callable[_P, _R]:
+        arg_types = tuple(type(arg) for arg in args)
+        for atype in arg_types:
+            # current strategy is to select the implementation of the first arg that supports the operation
+            # TODO: define a strategy that converts or prevents conversion
+            if (dispatcher := getattr(atype, "__gt_builtin_func__", None)) is not None and (
+                op_func := dispatcher(self)
+            ) is not NotImplemented:
+                return op_func
+        else:
+            return self.function
+
+    def __gt_type__(self) -> ts.FunctionType:
+        signature = inspect.signature(self.function)
+        params = signature.parameters
+
+        return ts.FunctionType(
+            pos_only_args=[
+                _type_conversion(param.annotation)
+                for param in params.values()
+                if param.kind == inspect.Parameter.POSITIONAL_ONLY
+            ],
+            pos_or_kw_args={
+                k: _type_conversion(v.annotation)
+                for k, v in params.items()
+                if v.kind == inspect.Parameter.POSITIONAL_OR_KEYWORD
+            },
+            kw_only_args={
+                k: _type_conversion(v.annotation)
+                for k, v in params.items()
+                if v.kind == inspect.Parameter.KEYWORD_ONLY
+            },
+            returns=_type_conversion(signature.return_annotation),
+        )
+
+
+def builtin_function(fun: Callable[_P, _R]) -> BuiltInFunction[_R, _P]:
+    return BuiltInFunction(fun)
+
+
+@builtin_function
+def neighbor_sum(
+    field: Field,
+    /,
+    axis: Dimension,
+) -> Field:
+    raise NotImplementedError()
+
+
+@builtin_function
+def max_over(
+    field: Field,
+    /,
+    axis: Dimension,
+) -> Field:
+    raise NotImplementedError()
+
+
+@builtin_function
+def min_over(
+    field: Field,
+    /,
+    axis: Dimension,
+) -> Field:
+    raise NotImplementedError()
+
+
+@builtin_function
+def broadcast(field: Field | gt4py_defs.ScalarT, dims: Tuple, /) -> Field:
+    raise NotImplementedError()
+
+
+@builtin_function
+def where(
+    mask: Field,
+    true_field: Field | gt4py_defs.ScalarT | Tuple,
+    false_field: Field | gt4py_defs.ScalarT | Tuple,
+    /,
+) -> Field | Tuple:
+    raise NotImplementedError()
+
+
+@builtin_function
+def astype(field: Field, type_: type, /) -> Field:
+    raise NotImplementedError()
 
-# unary math builtins (number) -> number
-abs = _unary_math_builtin  # noqa: A001
 
 UNARY_MATH_NUMBER_BUILTIN_NAMES = ["abs"]
 
-# unary math builtins (float) -> float
-sin = _unary_math_builtin
-cos = _unary_math_builtin
-tan = _unary_math_builtin
-arcsin = _unary_math_builtin
-arccos = _unary_math_builtin
-arctan = _unary_math_builtin
-sinh = _unary_math_builtin
-cosh = _unary_math_builtin
-tanh = _unary_math_builtin
-arcsinh = _unary_math_builtin
-arccosh = _unary_math_builtin
-arctanh = _unary_math_builtin
-sqrt = _unary_math_builtin
-exp = _unary_math_builtin
-log = _unary_math_builtin
-gamma = _unary_math_builtin
-cbrt = _unary_math_builtin
-floor = _unary_math_builtin
-ceil = _unary_math_builtin
-trunc = _unary_math_builtin
-
 UNARY_MATH_FP_BUILTIN_NAMES = [
     "sin",
     "cos",
@@ -157,42 +205,44 @@ UNARY_MATH_FP_BUILTIN_NAMES = [
     "trunc",
 ]
 
-# unary math predicates (float) -> bool
-_unary_math_predicate_builtin = BuiltInFunction(
-    ts.FunctionType(
-        pos_only_args=[ts.DeferredType(constraint=(ts.ScalarType, ts.FieldType))],
-        pos_or_kw_args={},
-        kw_only_args={},
-        returns=ts.DeferredType(constraint=(ts.ScalarType, ts.FieldType)),
-    )
-)
-
-isfinite = _unary_math_predicate_builtin
-isinf = _unary_math_predicate_builtin
-isnan = _unary_math_predicate_builtin
 
 UNARY_MATH_FP_PREDICATE_BUILTIN_NAMES = ["isfinite", "isinf", "isnan"]
 
-# binary math builtins (number, number) -> number
-_binary_math_builtin = BuiltInFunction(
-    ts.FunctionType(
-        pos_only_args=[
-            ts.DeferredType(constraint=(ts.ScalarType, ts.FieldType)),
-            ts.DeferredType(constraint=(ts.ScalarType, ts.FieldType)),
-        ],
-        pos_or_kw_args={},
-        kw_only_args={},
-        returns=ts.DeferredType(constraint=(ts.ScalarType, ts.FieldType)),
-    )
-)
 
-minimum = _binary_math_builtin
-maximum = _binary_math_builtin
-fmod = _binary_math_builtin
-power = _binary_math_builtin
+def _make_unary_math_builtin(name):
+    def impl(value: Field | gt4py_defs.ScalarT, /) -> Field | gt4py_defs.ScalarT:
+        raise NotImplementedError()
+
+    impl.__name__ = name
+    globals()[name] = builtin_function(impl)
+
+
+for f in (
+    UNARY_MATH_NUMBER_BUILTIN_NAMES
+    + UNARY_MATH_FP_BUILTIN_NAMES
+    + UNARY_MATH_FP_PREDICATE_BUILTIN_NAMES
+):
+    _make_unary_math_builtin(f)
+
 
 BINARY_MATH_NUMBER_BUILTIN_NAMES = ["minimum", "maximum", "fmod", "power"]
 
+
+def _make_binary_math_builtin(name):
+    def impl(
+        lhs: Field | gt4py_defs.ScalarT,
+        rhs: Field | gt4py_defs.ScalarT,
+        /,
+    ) -> Field | gt4py_defs.ScalarT:
+        raise NotImplementedError()
+
+    impl.__name__ = name
+    globals()[name] = builtin_function(impl)
+
+
+for f in BINARY_MATH_NUMBER_BUILTIN_NAMES:
+    _make_binary_math_builtin(f)
+
 MATH_BUILTIN_NAMES = (
     UNARY_MATH_NUMBER_BUILTIN_NAMES
     + UNARY_MATH_FP_BUILTIN_NAMES
@@ -221,10 +271,11 @@ __all__ = [*((set(BUILTIN_NAMES) | set(TYPE_ALIAS_NAMES)) - {"Dimension", "Field
 #  match. Revisit if we want to continue subclassing here. If we split
 #  them also check whether Dimension should continue to be the shared or define
 #  guidelines for decision.
-@dataclass(frozen=True)
+@dataclasses.dataclass(frozen=True)
 class FieldOffset(runtime.Offset):
     source: Dimension
     target: tuple[Dimension] | tuple[Dimension, Dimension]
+    connectivity: Optional[Any] = None  # TODO
 
     def __post_init__(self):
         if len(self.target) == 2 and self.target[1].kind != DimensionKind.LOCAL:
diff --git a/src/gt4py/next/ffront/type_info.py b/src/gt4py/next/ffront/type_info.py
index c81390ca..7f56f5d9 100644
--- a/src/gt4py/next/ffront/type_info.py
+++ b/src/gt4py/next/ffront/type_info.py
@@ -17,7 +17,7 @@ from typing import Iterator, cast
 
 import gt4py.next.ffront.type_specifications as ts_ffront
 import gt4py.next.type_system.type_specifications as ts
-from gt4py.next.common import Dimension
+from gt4py.next import common
 from gt4py.next.type_system import type_info
 
 
@@ -156,7 +156,7 @@ def _scan_param_promotion(param: ts.TypeSpec, arg: ts.TypeSpec) -> ts.FieldType
     --------
     >>> _scan_param_promotion(
     ...     ts.ScalarType(kind=ts.ScalarKind.INT64),
-    ...     ts.FieldType(dims=[Dimension("I")], dtype=ts.ScalarKind.FLOAT64)
+    ...     ts.FieldType(dims=[common.Dimension("I")], dtype=ts.ScalarKind.FLOAT64)
     ... )
     FieldType(dims=[Dimension(value='I', kind=<DimensionKind.HORIZONTAL: 'horizontal'>)], dtype=ScalarType(kind=<ScalarKind.INT64: 64>, shape=None))
     """
@@ -173,7 +173,7 @@ def _scan_param_promotion(param: ts.TypeSpec, arg: ts.TypeSpec) -> ts.FieldType
             # argument type differ. As such we can not extract the dimensions
             # and just return a generic field shown in the error later on.
             # TODO: we want some generic field type here, but our type system does not support it yet.
-            return ts.FieldType(dims=[Dimension("...")], dtype=dtype)
+            return ts.FieldType(dims=[common.Dimension("...")], dtype=dtype)
 
     return type_info.apply_to_primitive_constituents(param, _as_field, with_path_arg=True)
 
@@ -216,7 +216,7 @@ def function_signature_incompatibilities_scanop(
         for el in type_info.primitive_constituents(arg)
     ]
     try:
-        type_info.promote_dims(*arg_dims)
+        common.promote_dims(*arg_dims)
     except ValueError as e:
         yield e.args[0]
 
@@ -279,7 +279,7 @@ def return_type_scanop(
     with_kwargs: dict[str, ts.TypeSpec],
 ):
     carry_dtype = callable_type.definition.returns
-    promoted_dims = type_info.promote_dims(
+    promoted_dims = common.promote_dims(
         *(
             type_info.extract_dims(el)
             for arg in with_args + list(with_kwargs.values())
diff --git a/src/gt4py/next/iterator/ir_makers.py b/src/gt4py/next/iterator/ir_makers.py
index e74ce232..7a7ee993 100644
--- a/src/gt4py/next/iterator/ir_makers.py
+++ b/src/gt4py/next/iterator/ir_makers.py
@@ -14,7 +14,7 @@
 
 from typing import Callable, Union
 
-from gt4py.next import common
+from gt4py._core import definitions as core_defs
 from gt4py.next.iterator import ir as itir
 from gt4py.next.type_system import type_specifications as ts, type_translation
 
@@ -53,7 +53,7 @@ def ref(ref_or_name: Union[str, itir.SymRef]) -> itir.SymRef:
     return itir.SymRef(id=ref_or_name)
 
 
-def ensure_expr(literal_or_expr: Union[str, common.Scalar, itir.Expr]) -> itir.Expr:
+def ensure_expr(literal_or_expr: Union[str, core_defs.Scalar, itir.Expr]) -> itir.Expr:
     """
     Convert literals into a SymRef or Literal and let expressions pass unchanged.
 
@@ -70,9 +70,10 @@ def ensure_expr(literal_or_expr: Union[str, common.Scalar, itir.Expr]) -> itir.E
     """
     if isinstance(literal_or_expr, str):
         return ref(literal_or_expr)
-    elif isinstance(literal_or_expr, common.Scalar):  # type: ignore[arg-type] # mypy bug #11673
-        return literal_from_value(literal_or_expr)  # type: ignore[arg-type] # mypy bug #11673
-    return literal_or_expr  # type: ignore[return-value] # mypy bug #11673
+    elif core_defs.is_scalar_type(literal_or_expr):
+        return literal_from_value(literal_or_expr)
+    assert isinstance(literal_or_expr, itir.Expr)
+    return literal_or_expr
 
 
 def ensure_offset(str_or_offset: Union[str, int, itir.OffsetLiteral]) -> itir.OffsetLiteral:
@@ -275,7 +276,7 @@ def literal(value: str, typename: str):
     return itir.Literal(value=value, type=typename)
 
 
-def literal_from_value(val: common.Scalar) -> itir.Literal:
+def literal_from_value(val: core_defs.Scalar) -> itir.Literal:
     """
     Make a literal node from a value.
 
@@ -288,7 +289,7 @@ def literal_from_value(val: common.Scalar) -> itir.Literal:
     >>> literal_from_value(True)
     Literal(value='True', type='bool')
     """
-    if not isinstance(val, common.Scalar):  # type: ignore[arg-type] # mypy bug #11673
+    if not isinstance(val, core_defs.Scalar):  # type: ignore[arg-type] # mypy bug #11673
         raise ValueError(f"Value must be a scalar, but got {type(val).__name__}")
 
     # At the time this has been written the iterator module has its own type system that is
diff --git a/src/gt4py/next/iterator/tracing.py b/src/gt4py/next/iterator/tracing.py
index 4a96c8f4..81c97d02 100644
--- a/src/gt4py/next/iterator/tracing.py
+++ b/src/gt4py/next/iterator/tracing.py
@@ -16,8 +16,9 @@ import dataclasses
 import inspect
 from typing import List
 
+from gt4py._core import definitions as core_defs
 from gt4py.eve import Node
-from gt4py.next import common, iterator
+from gt4py.next import iterator
 from gt4py.next.iterator import builtins, ir_makers as im
 from gt4py.next.iterator.ir import (
     AxisLiteral,
@@ -151,7 +152,7 @@ def make_node(o):
             return lambdadef(o)
     if isinstance(o, iterator.runtime.Offset):
         return OffsetLiteral(value=o.value)
-    if isinstance(o, common.Scalar):
+    if isinstance(o, core_defs.Scalar):
         return im.literal_from_value(o)
     if isinstance(o, CartesianAxis):
         return AxisLiteral(value=o.value)
diff --git a/src/gt4py/next/type_system/type_info.py b/src/gt4py/next/type_system/type_info.py
index e4ce2e91..5a228247 100644
--- a/src/gt4py/next/type_system/type_info.py
+++ b/src/gt4py/next/type_system/type_info.py
@@ -19,7 +19,7 @@ from typing import Any, Callable, Iterator, Type, TypeGuard, cast
 import numpy as np
 
 from gt4py.eve.utils import XIterable, xiter
-from gt4py.next.common import Dimension, DimensionKind
+from gt4py.next import common
 from gt4py.next.type_system import type_specifications as ts
 
 
@@ -65,8 +65,8 @@ def primitive_constituents(
     """
     Return the primitive types contained in a composite type.
 
-    >>> from gt4py.next import Dimension
-    >>> I = Dimension(value="I")
+    >>> from gt4py.next import common
+    >>> I = common.Dimension(value="I")
     >>> int_type = ts.ScalarType(kind=ts.ScalarKind.INT64)
     >>> field_type = ts.FieldType(dims=[I], dtype=int_type)
 
@@ -275,7 +275,7 @@ def is_tuple_of_type(type_: ts.TypeSpec, expected_type: type | tuple) -> TypeGua
     return isinstance(type_, ts.TupleType) and is_type_or_tuple_of_type(type_, expected_type)
 
 
-def extract_dims(symbol_type: ts.TypeSpec) -> list[Dimension]:
+def extract_dims(symbol_type: ts.TypeSpec) -> list[common.Dimension]:
     """
     Try to extract field dimensions if possible.
 
@@ -285,8 +285,8 @@ def extract_dims(symbol_type: ts.TypeSpec) -> list[Dimension]:
     ---------
     >>> extract_dims(ts.ScalarType(kind=ts.ScalarKind.INT64, shape=[3, 4]))
     []
-    >>> I = Dimension(value="I")
-    >>> J = Dimension(value="J")
+    >>> I = common.Dimension(value="I")
+    >>> J = common.Dimension(value="J")
     >>> extract_dims(ts.FieldType(dims=[I, J], dtype=ts.ScalarType(kind=ts.ScalarKind.INT64)))
     [Dimension(value='I', kind=<DimensionKind.HORIZONTAL: 'horizontal'>), Dimension(value='J', kind=<DimensionKind.HORIZONTAL: 'horizontal'>)]
     """
@@ -304,14 +304,14 @@ def is_local_field(type_: ts.FieldType) -> bool:
 
     Examples:
     ---------
-    >>> V = Dimension(value="V")
-    >>> V2E = Dimension(value="V2E", kind=DimensionKind.LOCAL)
+    >>> V = common.Dimension(value="V")
+    >>> V2E = common.Dimension(value="V2E", kind=common.DimensionKind.LOCAL)
     >>> is_local_field(ts.FieldType(dims=[V, V2E], dtype=ts.ScalarType(kind=ts.ScalarKind.INT64)))
     True
     >>> is_local_field(ts.FieldType(dims=[V], dtype=ts.ScalarType(kind=ts.ScalarKind.INT64)))
     False
     """
-    return any(dim.kind == DimensionKind.LOCAL for dim in type_.dims)
+    return any(dim.kind == common.DimensionKind.LOCAL for dim in type_.dims)
 
 
 def contains_local_field(type_: ts.TypeSpec) -> bool:
@@ -379,10 +379,10 @@ def promote(*types: ts.FieldType | ts.ScalarType) -> ts.FieldType | ts.ScalarTyp
 
     The resulting type is defined on all dimensions of the arguments, respecting
     the individual order of the dimensions of each argument (see
-    :func:`promote_dims` for more details).
+    :func:`common.promote_dims` for more details).
 
     >>> dtype = ts.ScalarType(kind=ts.ScalarKind.INT64)
-    >>> I, J, K = (Dimension(value=dim) for dim in ["I", "J", "K"])
+    >>> I, J, K = (common.Dimension(value=dim) for dim in ["I", "J", "K"])
     >>> promoted: ts.FieldType = promote(
     ...     ts.FieldType(dims=[I, J], dtype=dtype),
     ...     ts.FieldType(dims=[I, J, K], dtype=dtype),
@@ -406,89 +406,13 @@ def promote(*types: ts.FieldType | ts.ScalarType) -> ts.FieldType | ts.ScalarTyp
             raise NotImplementedError("Shape promotion not implemented.")
         return types[0]
     elif all(isinstance(type_, (ts.ScalarType, ts.FieldType)) for type_ in types):
-        dims = promote_dims(*(extract_dims(type_) for type_ in types))
+        dims = common.promote_dims(*(extract_dims(type_) for type_ in types))
         dtype = cast(ts.ScalarType, promote(*(extract_dtype(type_) for type_ in types)))
 
         return ts.FieldType(dims=dims, dtype=dtype)
     raise TypeError("Expected a FieldType or ScalarType.")
 
 
-def promote_dims(*dims_list: list[Dimension]) -> list[Dimension]:
-    """
-    Find a unique ordering of multiple (individually ordered) lists of dimensions.
-
-    The resulting list of dimensions contains all dimensions of the arguments
-    in the order they originally appear. If no unique order exists or a
-    contradicting order is found an exception is raised.
-
-    A modified version (ensuring uniqueness of the order) of
-    `Kahn's algorithm <https://en.wikipedia.org/wiki/Topological_sorting#Kahn's_algorithm>`_
-    is used to topologically sort the arguments.
-
-    >>> I, J, K = (Dimension(value=dim) for dim in ["I", "J", "K"])
-    >>> promote_dims([I, J], [I, J, K]) == [I, J, K]
-    True
-    >>> promote_dims([I, J], [K]) # doctest: +ELLIPSIS
-    Traceback (most recent call last):
-     ...
-    ValueError: Dimensions can not be promoted. Could not determine order of the following dimensions: J, K.
-    >>> promote_dims([I, J], [J, I]) # doctest: +ELLIPSIS
-    Traceback (most recent call last):
-     ...
-    ValueError: Dimensions can not be promoted. The following dimensions appear in contradicting order: I, J.
-    """
-    # build a graph with the vertices being dimensions and edges representing
-    #  the order between two dimensions. The graph is encoded as a dictionary
-    #  mapping dimensions to their predecessors, i.e. a dictionary containing
-    #  adjacency lists. Since graphlib.TopologicalSorter uses predecessors
-    #  (contrary to successors) we also use this directionality here.
-    graph: dict[Dimension, set[Dimension]] = {}
-    for dims in dims_list:
-        if len(dims) == 0:
-            continue
-        # create a vertex for each dimension
-        for dim in dims:
-            graph.setdefault(dim, set())
-        # add edges
-        predecessor = dims[0]
-        for dim in dims[1:]:
-            graph[dim].add(predecessor)
-            predecessor = dim
-
-    # modified version of Kahn's algorithm
-    topologically_sorted_list: list[Dimension] = []
-
-    # compute in-degree for each vertex
-    in_degree = {v: 0 for v in graph.keys()}
-    for v1 in graph:
-        for v2 in graph[v1]:
-            in_degree[v2] += 1
-
-    # process vertices with in-degree == 0
-    # TODO(tehrengruber): avoid recomputation of zero_in_degree_vertex_list
-    while zero_in_degree_vertex_list := [v for v, d in in_degree.items() if d == 0]:
-        if len(zero_in_degree_vertex_list) != 1:
-            raise ValueError(
-                f"Dimensions can not be promoted. Could not determine "
-                f"order of the following dimensions: "
-                f"{', '.join((dim.value for dim in zero_in_degree_vertex_list))}."
-            )
-        v = zero_in_degree_vertex_list[0]
-        del in_degree[v]
-        topologically_sorted_list.insert(0, v)
-        # update in-degree
-        for predecessor in graph[v]:
-            in_degree[predecessor] -= 1
-
-    if len(in_degree.items()) > 0:
-        raise ValueError(
-            f"Dimensions can not be promoted. The following dimensions "
-            f"appear in contradicting order: {', '.join((dim.value for dim in in_degree.keys()))}."
-        )
-
-    return topologically_sorted_list
-
-
 @functools.singledispatch
 def return_type(
     callable_type: ts.CallableType,
@@ -690,7 +614,7 @@ def function_signature_incompatibilities_func(  # noqa: C901
             if i < len(func_type.pos_only_args):
                 arg_repr = f"{i}-th argument"
             else:
-                arg_repr = f"argument `{list(func_type.pos_or_kw_args.keys())[i-len(func_type.pos_only_args)]}`"
+                arg_repr = f"argument `{list(func_type.pos_or_kw_args.keys())[i - len(func_type.pos_only_args)]}`"
             yield f"Expected {arg_repr} to be of type `{a_arg}`, but got `{b_arg}`."
 
     for kwarg in set(func_type.kw_only_args.keys()) & set(kwargs.keys()):


Test patch:
diff --git a/tests/next_tests/integration_tests/feature_tests/ffront_tests/test_type_deduction.py b/tests/next_tests/integration_tests/feature_tests/ffront_tests/test_type_deduction.py
index c5013d1f..70b24ace 100644
--- a/tests/next_tests/integration_tests/feature_tests/ffront_tests/test_type_deduction.py
+++ b/tests/next_tests/integration_tests/feature_tests/ffront_tests/test_type_deduction.py
@@ -24,6 +24,7 @@ from gt4py.next import (
     FieldOffset,
     astype,
     broadcast,
+    common,
     errors,
     float32,
     float64,
@@ -461,51 +462,6 @@ def test_unpack_assign():
     )
 
 
-def dimension_promotion_cases() -> (
-    list[tuple[list[list[Dimension]], list[Dimension] | None, None | Pattern]]
-):
-    raw_list = [
-        # list of list of dimensions, expected result, expected error message
-        ([["I", "J"], ["I"]], ["I", "J"], None),
-        ([["I", "J"], ["J"]], ["I", "J"], None),
-        ([["I", "J"], ["J", "K"]], ["I", "J", "K"], None),
-        (
-            [["I", "J"], ["J", "I"]],
-            None,
-            r"The following dimensions appear in contradicting order: I, J.",
-        ),
-        (
-            [["I", "K"], ["J", "K"]],
-            None,
-            r"Could not determine order of the following dimensions: I, J",
-        ),
-    ]
-    # transform dimension names into Dimension objects
-    return [
-        (
-            [[Dimension(el) for el in arg] for arg in args],
-            [Dimension(el) for el in result] if result else result,
-            msg,
-        )
-        for args, result, msg in raw_list
-    ]
-
-
-@pytest.mark.parametrize("dim_list,expected_result,expected_error_msg", dimension_promotion_cases())
-def test_dimension_promotion(
-    dim_list: list[list[Dimension]],
-    expected_result: Optional[list[Dimension]],
-    expected_error_msg: Optional[str],
-):
-    if expected_result:
-        assert type_info.promote_dims(*dim_list) == expected_result
-    else:
-        with pytest.raises(Exception) as exc_info:
-            type_info.promote_dims(*dim_list)
-
-        assert exc_info.match(expected_error_msg)
-
-
 def test_assign_tuple():
     def temp_tuple(a: Field[[TDim], float64], b: Field[[TDim], int64]):
         tmp = a, b
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
new file mode 100644
index 00000000..6df39654
--- /dev/null
+++ b/tests/next_tests/unit_tests/embedded_tests/test_nd_array_field.py
@@ -0,0 +1,127 @@
+# 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 itertools
+import math
+import operator
+from typing import Callable, Iterable
+
+import numpy as np
+import pytest
+
+from gt4py.next import common
+from gt4py.next.embedded import nd_array_field
+from gt4py.next.ffront import fbuiltins
+
+from next_tests.integration_tests.feature_tests.math_builtin_test_data import math_builtin_test_data
+
+
+@pytest.fixture(params=nd_array_field._nd_array_implementations)
+def nd_array_implementation(request):
+    yield request.param
+
+
+@pytest.fixture(
+    params=[operator.add, operator.sub, operator.mul, operator.truediv, operator.floordiv],
+)
+def binary_op(request):
+    yield request.param
+
+
+def _make_field(lst: Iterable, nd_array_implementation):
+    return common.field(
+        nd_array_implementation.asarray(lst, dtype=nd_array_implementation.float32),
+        domain=((common.Dimension("foo"), common.UnitRange(0, len(lst))),),
+    )
+
+
+@pytest.mark.parametrize("builtin_name, inputs", math_builtin_test_data())
+def test_math_function_builtins(builtin_name: str, inputs, nd_array_implementation):
+    if builtin_name == "gamma":
+        # numpy has no gamma function
+        pytest.xfail("TODO: implement gamma")
+        ref_impl: Callable = np.vectorize(math.gamma)
+    else:
+        ref_impl: Callable = getattr(np, builtin_name)
+
+    expected = ref_impl(*[np.asarray(inp, dtype=np.float32) for inp in inputs])
+
+    field_inputs = [_make_field(inp, nd_array_implementation) for inp in inputs]
+
+    builtin = getattr(fbuiltins, builtin_name)
+    result = builtin(*field_inputs)
+
+    assert np.allclose(result.ndarray, expected)
+
+
+def test_binary_ops(binary_op, nd_array_implementation):
+    inp_a = [-1.0, 4.2, 42]
+    inp_b = [2.0, 3.0, -3.0]
+    inputs = [inp_a, inp_b]
+
+    expected = binary_op(*[np.asarray(inp, dtype=np.float32) for inp in inputs])
+
+    field_inputs = [_make_field(inp, nd_array_implementation) for inp in inputs]
+
+    result = binary_op(*field_inputs)
+
+    assert np.allclose(result.ndarray, expected)
+
+
+@pytest.fixture(
+    params=itertools.product(
+        nd_array_field._nd_array_implementations, nd_array_field._nd_array_implementations
+    ),
+    ids=lambda param: f"{param[0].__name__}-{param[1].__name__}",
+)
+def product_nd_array_implementation(request):
+    yield request.param
+
+
+def test_mixed_fields(product_nd_array_implementation):
+    first_impl, second_impl = product_nd_array_implementation
+    if (first_impl.__name__ == "cupy" and second_impl.__name__ == "numpy") or (
+        first_impl.__name__ == "numpy" and second_impl.__name__ == "cupy"
+    ):
+        pytest.skip("Binary operation between CuPy and NumPy requires explicit conversion.")
+
+    inp_a = [-1.0, 4.2, 42]
+    inp_b = [2.0, 3.0, -3.0]
+
+    expected = np.asarray(inp_a) + np.asarray(inp_b)
+
+    field_inp_a = _make_field(inp_a, first_impl)
+    field_inp_b = _make_field(inp_b, second_impl)
+
+    result = field_inp_a + field_inp_b
+    assert np.allclose(result.ndarray, expected)
+
+
+def test_non_dispatched_function():
+    @fbuiltins.builtin_function
+    def fma(a: common.Field, b: common.Field, c: common.Field, /) -> common.Field:
+        return a * b + c
+
+    inp_a = [-1.0, 4.2, 42]
+    inp_b = [2.0, 3.0, -3.0]
+    inp_c = [-2.0, -3.0, 3.0]
+
+    expected = np.asarray(inp_a) * np.asarray(inp_b) + np.asarray(inp_c)
+
+    field_inp_a = _make_field(inp_a, np)
+    field_inp_b = _make_field(inp_b, np)
+    field_inp_c = _make_field(inp_c, np)
+
+    result = fma(field_inp_a, field_inp_b, field_inp_c)
+    assert np.allclose(result.ndarray, expected)
diff --git a/tests/next_tests/unit_tests/test_common.py b/tests/next_tests/unit_tests/test_common.py
new file mode 100644
index 00000000..75c793f9
--- /dev/null
+++ b/tests/next_tests/unit_tests/test_common.py
@@ -0,0 +1,297 @@
+# 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 Optional, Pattern
+
+import pytest
+
+from gt4py.next.common import Dimension, DimensionKind, Domain, Infinity, UnitRange, promote_dims
+
+
+IDim = Dimension("IDim")
+ECDim = Dimension("ECDim")
+JDim = Dimension("JDim")
+KDim = Dimension("KDim", kind=DimensionKind.VERTICAL)
+
+
+@pytest.fixture
+def domain():
+    range1 = UnitRange(0, 10)
+    range2 = UnitRange(5, 15)
+    range3 = UnitRange(20, 30)
+
+    dimensions = (IDim, JDim, KDim)
+    ranges = (range1, range2, range3)
+
+    return Domain(dimensions, ranges)
+
+
+def test_empty_range():
+    expected = UnitRange(0, 0)
+    assert UnitRange(1, 1) == expected
+    assert UnitRange(1, -1) == expected
+
+
+@pytest.fixture
+def rng():
+    return UnitRange(-5, 5)
+
+
+def test_unit_range_length(rng):
+    assert rng.start == -5
+    assert rng.stop == 5
+    assert len(rng) == 10
+
+
+def test_unit_range_repr(rng):
+    assert repr(rng) == "UnitRange(-5, 5)"
+
+
+def test_unit_range_iter(rng):
+    actual = []
+    expected = list(range(-5, 5))
+
+    for elt in rng:
+        actual.append(elt)
+
+    assert actual == expected
+
+
+def test_unit_range_get_item(rng):
+    assert rng[-1] == 4
+    assert rng[0] == -5
+    assert rng[0:4] == UnitRange(-5, -1)
+    assert rng[-4:] == UnitRange(1, 5)
+
+
+def test_unit_range_index_error(rng):
+    with pytest.raises(IndexError):
+        rng[10]
+
+
+def test_unit_range_slice_error(rng):
+    with pytest.raises(ValueError):
+        rng[1:2:5]
+
+
+def test_unit_range_set_intersection(rng):
+    with pytest.raises(
+        NotImplementedError, match="Can only find the intersection between UnitRange instances."
+    ):
+        rng & {1, 5}
+
+
+@pytest.mark.parametrize(
+    "rng1, rng2, expected",
+    [
+        (UnitRange(0, 5), UnitRange(10, 15), UnitRange(0, 0)),
+        (UnitRange(0, 5), UnitRange(5, 10), UnitRange(5, 5)),
+        (UnitRange(0, 5), UnitRange(3, 7), UnitRange(3, 5)),
+        (UnitRange(0, 5), UnitRange(1, 6), UnitRange(1, 5)),
+        (UnitRange(0, 5), UnitRange(-5, 5), UnitRange(0, 5)),
+        (UnitRange(0, 0), UnitRange(0, 5), UnitRange(0, 0)),
+        (UnitRange(0, 0), UnitRange(0, 0), UnitRange(0, 0)),
+    ],
+)
+def test_unit_range_intersection(rng1, rng2, expected):
+    result = rng1 & rng2
+    assert result == expected
+
+
+@pytest.mark.parametrize(
+    "rng1, rng2, expected",
+    [
+        (UnitRange(20, Infinity.positive()), UnitRange(10, 15), UnitRange(0, 0)),
+        (UnitRange(Infinity.negative(), 0), UnitRange(5, 10), UnitRange(0, 0)),
+        (UnitRange(Infinity.negative(), 0), UnitRange(-10, 0), UnitRange(-10, 0)),
+        (UnitRange(0, Infinity.positive()), UnitRange(Infinity.negative(), 5), UnitRange(0, 5)),
+        (
+            UnitRange(Infinity.negative(), 0),
+            UnitRange(Infinity.negative(), 5),
+            UnitRange(Infinity.negative(), 0),
+        ),
+        (
+            UnitRange(Infinity.negative(), Infinity.positive()),
+            UnitRange(Infinity.negative(), Infinity.positive()),
+            UnitRange(Infinity.negative(), Infinity.positive()),
+        ),
+    ],
+)
+def test_unit_range_infinite_intersection(rng1, rng2, expected):
+    result = rng1 & rng2
+    assert result == expected
+
+
+def test_positive_infinity_range():
+    pos_inf_range = UnitRange(Infinity.positive(), Infinity.positive())
+    assert len(pos_inf_range) == 0
+
+
+def test_mixed_infinity_range():
+    mixed_inf_range = UnitRange(Infinity.negative(), Infinity.positive())
+    assert len(mixed_inf_range) == Infinity.positive()
+
+
+def test_domain_length(domain):
+    assert len(domain) == 3
+
+
+def test_domain_iteration(domain):
+    iterated_values = [val for val in domain]
+    assert iterated_values == list(zip(domain.dims, domain.ranges))
+
+
+def test_domain_contains_named_range(domain):
+    assert (IDim, UnitRange(0, 10)) in domain
+    assert (IDim, UnitRange(-5, 5)) not in domain
+
+
+@pytest.mark.parametrize(
+    "second_domain, expected",
+    [
+        (
+            Domain((IDim, JDim), (UnitRange(2, 12), UnitRange(7, 17))),
+            Domain((IDim, JDim, KDim), (UnitRange(2, 10), UnitRange(7, 15), UnitRange(20, 30))),
+        ),
+        (
+            Domain((IDim, KDim), (UnitRange(2, 12), UnitRange(7, 27))),
+            Domain((IDim, JDim, KDim), (UnitRange(2, 10), UnitRange(5, 15), UnitRange(20, 27))),
+        ),
+        (
+            Domain((JDim, KDim), (UnitRange(2, 12), UnitRange(4, 27))),
+            Domain((IDim, JDim, KDim), (UnitRange(0, 10), UnitRange(5, 12), UnitRange(20, 27))),
+        ),
+    ],
+)
+def test_domain_intersection_different_dimensions(domain, second_domain, expected):
+    result_domain = domain & second_domain
+    print(result_domain)
+
+    assert result_domain == expected
+
+
+def test_domain_intersection_reversed_dimensions(domain):
+    dimensions = (JDim, IDim)
+    ranges = (UnitRange(2, 12), UnitRange(7, 17))
+    domain2 = Domain(dimensions, ranges)
+
+    with pytest.raises(
+        ValueError,
+        match="Dimensions can not be promoted. The following dimensions appear in contradicting order: IDim, JDim.",
+    ):
+        domain & domain2
+
+
+@pytest.mark.parametrize(
+    "index, expected",
+    [
+        (0, (IDim, UnitRange(0, 10))),
+        (1, (JDim, UnitRange(5, 15))),
+        (2, (KDim, UnitRange(20, 30))),
+        (-1, (KDim, UnitRange(20, 30))),
+        (-2, (JDim, UnitRange(5, 15))),
+    ],
+)
+def test_domain_integer_indexing(domain, index, expected):
+    result = domain[index]
+    assert result == expected
+
+
+@pytest.mark.parametrize(
+    "slice_obj, expected",
+    [
+        (slice(0, 2), ((IDim, UnitRange(0, 10)), (JDim, UnitRange(5, 15)))),
+        (slice(1, None), ((JDim, UnitRange(5, 15)), (KDim, UnitRange(20, 30)))),
+    ],
+)
+def test_domain_slice_indexing(domain, slice_obj, expected):
+    result = domain[slice_obj]
+    assert isinstance(result, Domain)
+    assert len(result) == len(expected)
+    assert all(res == exp for res, exp in zip(result, expected))
+
+
+@pytest.mark.parametrize(
+    "index, expected_result",
+    [
+        (JDim, (JDim, UnitRange(5, 15))),
+        (KDim, (KDim, UnitRange(20, 30))),
+    ],
+)
+def test_domain_dimension_indexing(domain, index, expected_result):
+    result = domain[index]
+    assert result == expected_result
+
+
+def test_domain_indexing_dimension_missing(domain):
+    with pytest.raises(KeyError, match=r"No Dimension of type .* is present in the Domain."):
+        domain[ECDim]
+
+
+def test_domain_indexing_invalid_type(domain):
+    with pytest.raises(
+        KeyError, match="Invalid index type, must be either int, slice, or Dimension."
+    ):
+        domain["foo"]
+
+
+def test_domain_repeat_dims():
+    dims = (IDim, JDim, IDim)
+    ranges = (UnitRange(0, 5), UnitRange(0, 8), UnitRange(0, 3))
+    with pytest.raises(NotImplementedError, match=r"Domain dimensions must be unique, not .*"):
+        Domain(dims, ranges)
+
+
+def dimension_promotion_cases() -> (
+    list[tuple[list[list[Dimension]], list[Dimension] | None, None | Pattern]]
+):
+    raw_list = [
+        # list of list of dimensions, expected result, expected error message
+        ([["I", "J"], ["I"]], ["I", "J"], None),
+        ([["I", "J"], ["J"]], ["I", "J"], None),
+        ([["I", "J"], ["J", "K"]], ["I", "J", "K"], None),
+        (
+            [["I", "J"], ["J", "I"]],
+            None,
+            r"The following dimensions appear in contradicting order: I, J.",
+        ),
+        (
+            [["I", "K"], ["J", "K"]],
+            None,
+            r"Could not determine order of the following dimensions: I, J",
+        ),
+    ]
+    # transform dimension names into Dimension objects
+    return [
+        (
+            [[Dimension(el) for el in arg] for arg in args],
+            [Dimension(el) for el in result] if result else result,
+            msg,
+        )
+        for args, result, msg in raw_list
+    ]
+
+
+@pytest.mark.parametrize("dim_list,expected_result,expected_error_msg", dimension_promotion_cases())
+def test_dimension_promotion(
+    dim_list: list[list[Dimension]],
+    expected_result: Optional[list[Dimension]],
+    expected_error_msg: Optional[str],
+):
+    if expected_result:
+        assert promote_dims(*dim_list) == expected_result
+    else:
+        with pytest.raises(Exception) as exc_info:
+            promote_dims(*dim_list)
+
+        assert exc_info.match(expected_error_msg)