Repository: GridTools/gt4py
Pull request #1303 (2023-07-27 13:54:27)

feat: low-level buffer allocation utilities for new storages
## Description

Add new buffer allocation utilities and concepts definitions.
Main changes:
- New definition of core concepts dealing with data: DType, Device, ScalarTypes
- New BullerAllocator protocol to support different buffer allocation strategies depending on the device type.
- Refactoring of existing cartesian storages allocation tools using the new buffer allocators

## Pending
- [x] Make current unit tests pass
- [ ] Add unit tests for new buffer allocation and core definitions
- [x] Fix typing with python 3.8 compatibility and code style errors


havogt: cscs-ci run
havogt: cscs-ci run
havogt: cscs-ci run
havogt: cscs-ci run
havogt: cscs-ci run
havogt: cscs-ci run
havogt: cscs-ci run

Patch:
diff --git a/ci/base.Dockerfile b/ci/base.Dockerfile
index f32fe4cb..2de6aa7f 100644
--- a/ci/base.Dockerfile
+++ b/ci/base.Dockerfile
@@ -54,5 +54,4 @@ RUN pyenv update && \
 ENV PATH="/root/.pyenv/shims:${PATH}"
 
 
-RUN pip install --upgrade pip setuptools wheel tox cupy-cuda112
-
+RUN pip install --upgrade pip setuptools wheel tox cupy-cuda11x
diff --git a/ci/cscs-ci.yml b/ci/cscs-ci.yml
index 4474069a..3dc38bcd 100644
--- a/ci/cscs-ci.yml
+++ b/ci/cscs-ci.yml
@@ -88,7 +88,7 @@ test py38:
     matrix:
     - SUBPACKAGE: [cartesian, storage]
       VARIANT: [-internal, -dace]
-      SUBVARIANT: [-cuda112, -cpu]
+      SUBVARIANT: [-cuda11x, -cpu]
     - SUBPACKAGE: eve
   variables:
     CRAY_CUDA_MPS: 1
@@ -113,7 +113,7 @@ test py310:
     matrix:
     - SUBPACKAGE: [cartesian, storage]
       VARIANT: [-internal, -dace]
-      SUBVARIANT: [-cuda112, -cpu]
+      SUBVARIANT: [-cuda11x, -cpu]
     - SUBPACKAGE: eve
     - SUBPACKAGE: next
       VARIANT: [-nomesh, -atlas]
diff --git a/pyproject.toml b/pyproject.toml
index 38505ad5..acc04ee8 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -66,16 +66,9 @@ readme = 'README.md'
 requires-python = '>=3.8'
 
 [project.optional-dependencies]
-cuda = ['cupy']
-cuda110 = ['cupy-cuda110']
-cuda111 = ['cupy-cuda111']
-cuda112 = ['cupy-cuda112']
-cuda113 = ['cupy-cuda113']
-cuda114 = ['cupy-cuda114']
-cuda115 = ['cupy-cuda115']
-cuda116 = ['cupy-cuda116']
-cuda117 = ['cupy-cuda117']
-cuda11x = ['cupy-cuda11x']
+cuda = ['cupy>=12.0']
+cuda11x = ['cupy-cuda11x>=12.0']
+cuda12x = ['cupy-cuda12x>=12.0']
 dace = ['dace>=0.14.2,<0.15', 'sympy>=1.7']
 formatting = ['clang-format>=9.0']
 # Always add all extra packages to 'full' for a simple full gt4py installation
diff --git a/src/gt4py/_core/__init__.py b/src/gt4py/_core/__init__.py
new file mode 100644
index 00000000..6c43e2f1
--- /dev/null
+++ b/src/gt4py/_core/__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/_core/definitions.py b/src/gt4py/_core/definitions.py
new file mode 100644
index 00000000..d6b34c42
--- /dev/null
+++ b/src/gt4py/_core/definitions.py
@@ -0,0 +1,439 @@
+# 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 collections.abc
+import dataclasses
+import enum
+import functools
+import math
+import numbers
+from typing import overload
+
+import numpy as np
+import numpy.typing as npt
+
+from gt4py.eve.extended_typing import (
+    TYPE_CHECKING,
+    Any,
+    Final,
+    Generic,
+    Iterator,
+    Literal,
+    Protocol,
+    Sequence,
+    Tuple,
+    Type,
+    TypeAlias,
+    TypeGuard,
+    TypeVar,
+    Union,
+    cast,
+)
+
+
+if TYPE_CHECKING:
+    import cupy as cp
+
+    CuPyNDArray: TypeAlias = cp.ndarray
+
+    import jax.numpy as jnp
+
+    JaxNDArray: TypeAlias = jnp.ndarray
+
+
+# -- Scalar types supported by GT4Py --
+bool_ = np.bool_
+
+int8 = np.int8
+int16 = np.int16
+int32 = np.int32
+int64 = np.int64
+
+uint8 = np.uint8
+uint16 = np.uint16
+uint32 = np.uint32
+uint64 = np.uint64
+
+float32 = np.float32
+float64 = np.float64
+
+BoolScalar: TypeAlias = Union[bool_, bool]
+BoolT = TypeVar("BoolT", bound=Union[bool_, bool])
+BOOL_TYPES: Final[Tuple[type, ...]] = cast(Tuple[type, ...], BoolScalar.__args__)  # type: ignore[attr-defined]
+
+
+IntScalar: TypeAlias = Union[int8, int16, int32, int64, int]
+IntT = TypeVar("IntT", bound=Union[int8, int16, int32, int64, int])
+INT_TYPES: Final[Tuple[type, ...]] = cast(Tuple[type, ...], IntScalar.__args__)  # type: ignore[attr-defined]
+
+
+UnsignedIntScalar: TypeAlias = Union[uint8, uint16, uint32, uint64]
+UnsignedIntT = TypeVar("UnsignedIntT", bound=Union[uint8, uint16, uint32, uint64])
+UINT_TYPES: Final[Tuple[type, ...]] = cast(Tuple[type, ...], UnsignedIntScalar.__args__)  # type: ignore[attr-defined]
+
+
+IntegralScalar: TypeAlias = Union[IntScalar, UnsignedIntScalar]
+IntegralT = TypeVar("IntegralT", bound=Union[IntScalar, UnsignedIntScalar])
+INTEGRAL_TYPES: Final[Tuple[type, ...]] = (*INT_TYPES, *UINT_TYPES)
+
+
+FloatingScalar: TypeAlias = Union[float32, float64, float]
+FloatingT = TypeVar("FloatingT", bound=Union[float32, float64, float])
+FLOAT_TYPES: Final[Tuple[type, ...]] = cast(Tuple[type, ...], FloatingScalar.__args__)  # type: ignore[attr-defined]
+
+
+#: Type alias for all scalar types supported by GT4Py
+Scalar: TypeAlias = Union[BoolScalar, IntegralScalar, FloatingScalar]
+ScalarT = TypeVar("ScalarT", bound=Union[BoolScalar, IntegralScalar, FloatingScalar])
+SCALAR_TYPES: Final[tuple[type, ...]] = (*BOOL_TYPES, *INTEGRAL_TYPES, *FLOAT_TYPES)
+
+
+class BooleanIntegral(numbers.Integral):
+    """Abstract base class for boolean integral types."""
+
+    ...
+
+
+class PositiveIntegral(numbers.Integral):
+    """Abstract base class representing positive integral numbers.."""
+
+    ...
+
+
+def is_boolean_integral_type(integral_type: type) -> TypeGuard[Type[BooleanIntegral]]:
+    return issubclass(integral_type, BOOL_TYPES)
+
+
+def is_positive_integral_type(integral_type: type) -> TypeGuard[Type[PositiveIntegral]]:
+    return issubclass(integral_type, UINT_TYPES)
+
+
+TensorShape: TypeAlias = Sequence[
+    int
+]  # TODO(egparedes) figure out if PositiveIntegral can be made to work
+
+
+def is_valid_tensor_shape(
+    value: Sequence[IntegralScalar],
+) -> TypeGuard[TensorShape]:
+    return isinstance(value, collections.abc.Sequence) and all(
+        isinstance(v, numbers.Integral) and v > 0 for v in value
+    )
+
+
+# -- Data type descriptors --
+class DTypeKind(enum.Enum):
+    """
+    Kind of a specific data type.
+
+    Character codes match the type code for the corresponding kind in the
+    array interface protocol (https://numpy.org/doc/stable/reference/arrays.interface.html).
+    """
+
+    BOOL = "b"
+    INT = "i"
+    UINT = "u"
+    FLOAT = "f"
+    OPAQUE_POINTER = "V"
+    COMPLEX = "c"
+
+
+@overload
+def dtype_kind(sc_type: Type[BoolT]) -> Literal[DTypeKind.BOOL]:
+    ...
+
+
+@overload
+def dtype_kind(sc_type: Type[IntT]) -> Literal[DTypeKind.INT]:
+    ...
+
+
+@overload
+def dtype_kind(sc_type: Type[UnsignedIntT]) -> Literal[DTypeKind.UINT]:
+    ...
+
+
+@overload
+def dtype_kind(sc_type: Type[FloatingT]) -> Literal[DTypeKind.FLOAT]:
+    ...
+
+
+@overload
+def dtype_kind(sc_type: Type[ScalarT]) -> DTypeKind:
+    ...
+
+
+def dtype_kind(sc_type: Type[ScalarT]) -> DTypeKind:
+    """Return the data type kind of the given scalar type."""
+    if issubclass(sc_type, numbers.Integral):
+        if is_boolean_integral_type(sc_type):
+            return DTypeKind.BOOL
+        elif is_positive_integral_type(sc_type):
+            return DTypeKind.UINT
+        else:
+            return DTypeKind.INT
+    if issubclass(sc_type, numbers.Real):
+        return DTypeKind.FLOAT
+    if issubclass(sc_type, numbers.Complex):
+        return DTypeKind.COMPLEX
+
+    raise TypeError("Unknown scalar type kind")
+
+
+@dataclasses.dataclass(frozen=True)
+class DType(Generic[ScalarT]):
+    """
+    Descriptor of data type for Field elements.
+
+    This definition is based on DLPack and Array API standards. The data type
+    should be interpreted as packed `lanes` repetitions of elements from
+    `kind` data-category of `bit_width` width.
+
+    The Array API standard only requires DTypes to be comparable with `__eq__`.
+
+    Additionally, instances of this class can also be used as valid NumPy
+    `dtype`s definitions due to the `.dtype` attribute.
+    """
+
+    scalar_type: Type[ScalarT]
+    tensor_shape: TensorShape
+
+    def __init__(
+        self, scalar_type: Type[ScalarT], tensor_shape: Sequence[IntegralScalar] = ()
+    ) -> None:
+        if not isinstance(scalar_type, type):
+            raise TypeError(f"Invalid scalar type '{scalar_type}'")
+        if not is_valid_tensor_shape(tensor_shape):
+            raise TypeError(f"Invalid tensor shape '{tensor_shape}'")
+
+        object.__setattr__(self, "scalar_type", scalar_type)
+        object.__setattr__(self, "tensor_shape", tensor_shape)
+
+    @functools.cached_property
+    def kind(self) -> DTypeKind:
+        return dtype_kind(self.scalar_type)
+
+    @functools.cached_property
+    def dtype(self) -> np.dtype:
+        """The NumPy dtype corresponding to this DType."""
+        return np.dtype(f"{self.tensor_shape}{np.dtype(self.scalar_type).name}")
+
+    @functools.cached_property
+    def byte_size(self) -> int:
+        return np.dtype(self.scalar_type).itemsize * self.lanes
+
+    @property
+    def bit_width(self) -> int:
+        return 8 * np.dtype(self.scalar_type).itemsize
+
+    @property
+    def lanes(self) -> int:
+        return math.prod(self.tensor_shape or (1,))
+
+    @property
+    def subndim(self) -> int:
+        return len(self.tensor_shape)
+
+
+@dataclasses.dataclass(frozen=True)
+class IntegerDType(DType[IntegralT]):
+    pass
+
+
+@dataclasses.dataclass(frozen=True)
+class UnsignedIntDType(DType[UnsignedIntT]):
+    pass
+
+
+@dataclasses.dataclass(frozen=True)
+class UInt8DType(UnsignedIntDType[uint8]):
+    scalar_type: Final[Type[uint8]] = dataclasses.field(default=uint8, init=False)
+
+
+@dataclasses.dataclass(frozen=True)
+class UInt16DType(UnsignedIntDType[uint16]):
+    scalar_type: Final[Type[uint16]] = dataclasses.field(default=uint16, init=False)
+
+
+@dataclasses.dataclass(frozen=True)
+class UInt32DType(UnsignedIntDType[uint32]):
+    scalar_type: Final[Type[uint32]] = dataclasses.field(default=uint32, init=False)
+
+
+@dataclasses.dataclass(frozen=True)
+class UInt64DType(UnsignedIntDType[uint64]):
+    scalar_type: Final[Type[uint64]] = dataclasses.field(default=uint64, init=False)
+
+
+@dataclasses.dataclass(frozen=True)
+class SignedIntDType(DType[IntT]):
+    pass
+
+
+@dataclasses.dataclass(frozen=True)
+class Int8DType(SignedIntDType[int8]):
+    scalar_type: Final[Type[int8]] = dataclasses.field(default=int8, init=False)
+
+
+@dataclasses.dataclass(frozen=True)
+class Int16DType(SignedIntDType[int16]):
+    scalar_type: Final[Type[int16]] = dataclasses.field(default=int16, init=False)
+
+
+@dataclasses.dataclass(frozen=True)
+class Int32DType(SignedIntDType[int32]):
+    scalar_type: Final[Type[int32]] = dataclasses.field(default=int32, init=False)
+
+
+@dataclasses.dataclass(frozen=True)
+class Int64DType(SignedIntDType[int64]):
+    scalar_type: Final[Type[int64]] = dataclasses.field(default=int64, init=False)
+
+
+@dataclasses.dataclass(frozen=True)
+class FloatingDType(DType[FloatingT]):
+    pass
+
+
+@dataclasses.dataclass(frozen=True)
+class Float32DType(FloatingDType[float32]):
+    scalar_type: Final[Type[float32]] = dataclasses.field(default=float32, init=False)
+
+
+@dataclasses.dataclass(frozen=True)
+class Float64DType(FloatingDType[float64]):
+    scalar_type: Final[Type[float64]] = dataclasses.field(default=float64, init=False)
+
+
+DTypeLike = Union[DType, npt.DTypeLike]
+
+
+def dtype(dtype_like: DTypeLike) -> DType:
+    """Return the DType corresponding to the given dtype-like object."""
+    return dtype_like if isinstance(dtype_like, DType) else DType(np.dtype(dtype_like).type)
+
+
+# -- Custom protocols  --
+class GTDimsInterface(Protocol):
+    __gt_dims__: Tuple[str, ...]
+
+
+class GTOriginInterface(Protocol):
+    __gt_origin__: Tuple[int, ...]
+
+
+# -- Device representation --
+class DeviceType(enum.Enum):
+    """The type of the device where a memory buffer is allocated.
+
+    Enum values taken from DLPack reference implementation at:
+    https://github.com/dmlc/dlpack/blob/main/include/dlpack/dlpack.h
+    """
+
+    CPU = 1
+    CUDA = 2
+    CPU_PINNED = 3
+    OPENCL = 4
+    VULKAN = 7
+    METAL = 8
+    VPI = 9
+    ROCM = 10
+
+
+@dataclasses.dataclass(frozen=True)
+class Device:
+    """
+    Representation of a computing device.
+
+    This definition is based on the DLPack device definition. A device is
+    described by a pair of `DeviceType` and `device_id`. The `device_id`
+    is an integer that is interpreted differently depending on the
+    `DeviceType`. For example, for `DeviceType.CPU` it could be the CPU
+    core number, for `DeviceType.CUDA` it could be the CUDA device number, etc.
+    """
+
+    device_type: DeviceType
+    device_id: int
+
+    def __iter__(self) -> Iterator[DeviceType | int]:
+        yield self.device_type
+        yield self.device_id
+
+
+# -- NDArrays and slices --
+SliceLike = Union[int, Tuple[int, ...], None, slice, "NDArrayObject"]
+
+
+class NDArrayObjectProto(Protocol):
+    @property
+    def ndim(self) -> int:
+        ...
+
+    @property
+    def shape(self) -> tuple[int, ...]:
+        ...
+
+    @property
+    def dtype(self) -> Any:
+        ...
+
+    def __getitem__(self, item: SliceLike) -> NDArrayObjectProto:
+        ...
+
+    def __abs__(self) -> NDArrayObjectProto:
+        ...
+
+    def __neg__(self) -> NDArrayObjectProto:
+        ...
+
+    def __add__(self, other: NDArrayObject | Scalar) -> NDArrayObjectProto:
+        ...
+
+    def __radd__(self, other: NDArrayObject | Scalar) -> NDArrayObjectProto:
+        ...
+
+    def __sub__(self, other: NDArrayObject | Scalar) -> NDArrayObjectProto:
+        ...
+
+    def __rsub__(self, other: NDArrayObject | Scalar) -> NDArrayObjectProto:
+        ...
+
+    def __mul__(self, other: NDArrayObject | Scalar) -> NDArrayObjectProto:
+        ...
+
+    def __rmul__(self, other: NDArrayObject | Scalar) -> NDArrayObjectProto:
+        ...
+
+    def __floordiv__(self, other: NDArrayObject | Scalar) -> NDArrayObjectProto:
+        ...
+
+    def __rfloordiv__(self, other: NDArrayObject | Scalar) -> NDArrayObjectProto:
+        ...
+
+    def __truediv__(self, other: NDArrayObject | Scalar) -> NDArrayObjectProto:
+        ...
+
+    def __rtruediv__(self, other: NDArrayObject | Scalar) -> NDArrayObjectProto:
+        ...
+
+    def __pow__(self, other: NDArrayObject | Scalar) -> NDArrayObjectProto:
+        ...
+
+
+NDArrayObject = Union[npt.NDArray, "CuPyNDArray", "JaxNDArray", NDArrayObjectProto]
+NDArrayObjectT = TypeVar(
+    "NDArrayObjectT", npt.NDArray, "CuPyNDArray", "JaxNDArray", NDArrayObjectProto, covariant=True
+)
diff --git a/src/gt4py/cartesian/stencil_object.py b/src/gt4py/cartesian/stencil_object.py
index 59089693..69ce980b 100644
--- a/src/gt4py/cartesian/stencil_object.py
+++ b/src/gt4py/cartesian/stencil_object.py
@@ -12,6 +12,8 @@
 #
 # SPDX-License-Identifier: GPL-3.0-or-later
 
+from __future__ import annotations
+
 import abc
 import collections.abc
 import sys
@@ -19,12 +21,12 @@ import time
 from dataclasses import dataclass
 from numbers import Number
 from pickle import dumps
-from typing import Any, Callable, ClassVar, Dict, Optional, Tuple, Union, cast
+from typing import Any, Callable, ClassVar, Dict, Literal, Optional, Tuple, Union, cast
 
 import numpy as np
 
 import gt4py.cartesian.gtc.utils as gtc_utils
-import gt4py.storage.utils as storage_utils
+import gt4py.storage.cartesian.utils as storage_utils
 from gt4py import cartesian as gt4pyc
 from gt4py.cartesian.definitions import AccessKind, DomainInfo, FieldInfo, ParameterInfo
 from gt4py.cartesian.gtc.definitions import Index, Shape
@@ -39,17 +41,15 @@ FieldType = Union["cp.ndarray", np.ndarray]
 OriginType = Union[Tuple[int, int, int], Dict[str, Tuple[int, ...]]]
 
 
-def _compute_cache_key(
-    field_args: Dict[str, Optional[FieldType]],
+def _compute_domain_origin_cache_key(
+    field_args_info: Dict[str, Optional[ArgsInfo]],
     parameter_args: Dict[str, Optional[Number]],
     domain: Optional[Tuple[int, ...]],
     origin: Optional[OriginType],
-    device: str,
 ) -> int:
-    asarray = storage_utils.as_cupy if device == "gpu" else storage_utils.as_numpy
     field_data = tuple(
-        (name, asarray(arg).shape, getattr(arg, "__gt_origin__", (0, 0, 0)))
-        for name, arg in field_args.items()
+        (name, arg.array.shape, arg.origin or (0, 0, 0))
+        for name, arg in field_args_info.items()
         if arg is not None
     )
     return hash((field_data, *parameter_args.keys(), dumps(domain), dumps(origin)))
@@ -65,15 +65,14 @@ class ArgsInfo:
 
 
 def _extract_array_infos(
-    field_args: Dict[str, Optional[FieldType]], device: str
+    field_args: Dict[str, Optional[FieldType]], device: Literal["cpu", "gpu"]
 ) -> Dict[str, Optional[ArgsInfo]]:
-    asarray = storage_utils.as_cupy if device == "gpu" else storage_utils.as_numpy
     array_infos: Dict[str, Optional[ArgsInfo]] = {}
     for name, arg in field_args.items():
         if arg is None:
             array_infos[name] = None
         else:
-            array = asarray(arg)
+            array = storage_utils.asarray(arg, device=device)
             dimensions = storage_utils.get_dims(arg)
             if dimensions is not None:
                 sorted_dimensions = [d for d in "IJK" if d in dimensions]
@@ -571,7 +570,7 @@ class StencilObject(abc.ABC):
         device = backend_cls.storage_info["device"]
         array_infos = _extract_array_infos(field_args, device)
 
-        cache_key = _compute_cache_key(field_args, parameter_args, domain, origin, device)
+        cache_key = _compute_domain_origin_cache_key(array_infos, parameter_args, domain, origin)
         if cache_key not in self._domain_origin_cache:
             origin = self._normalize_origins(array_infos, self.field_info, origin)
 
diff --git a/src/gt4py/eve/extended_typing.py b/src/gt4py/eve/extended_typing.py
index 2002ae14..5daa7e1f 100644
--- a/src/gt4py/eve/extended_typing.py
+++ b/src/gt4py/eve/extended_typing.py
@@ -12,7 +12,8 @@
 #
 # SPDX-License-Identifier: GPL-3.0-or-later
 
-"""Typing definitions working across different Python versions (via `typing_extensions`).
+"""
+Typing definitions working across different Python versions (via `typing_extensions`).
 
 Definitions in 'typing_extensions' take priority over those in 'typing'.
 """
@@ -22,6 +23,7 @@ from __future__ import annotations
 import abc as _abc
 import array as _array
 import collections.abc as _collections_abc
+import ctypes as _ctypes
 import dataclasses as _dataclasses
 import enum as _enum
 import functools as _functools
@@ -234,6 +236,56 @@ class HashlibAlgorithm(Protocol):
 
 
 # -- Third party protocols --
+class ArrayInterfaceTypedDict(TypedDict):
+    shape: Tuple[int, ...]
+    typestr: str
+    descr: NotRequired[List[Tuple]]
+    data: NotRequired[Tuple[int, bool]]
+    strides: NotRequired[Optional[Tuple[int, ...]]]
+    mask: NotRequired[Optional["StrictArrayInterface"]]
+    offset: NotRequired[int]
+    version: int
+
+
+class StrictArrayInterface(Protocol):
+    __array_interface__: ArrayInterfaceTypedDict
+
+
+class ArrayInterface(Protocol):
+    __array_interface__: Dict[str, Any]
+
+
+class CUDAArrayInterfaceTypedDict(TypedDict):
+    shape: Tuple[int, ...]
+    typestr: str
+    data: Tuple[int, bool]
+    version: int
+    strides: NotRequired[Optional[Tuple[int, ...]]]
+    descr: NotRequired[List[Tuple]]
+    mask: NotRequired[Optional["StrictCUDAArrayInterface"]]
+    stream: NotRequired[Optional[int]]
+
+
+class StrictCUDAArrayInterface(Protocol):
+    __cuda_array_interface__: CUDAArrayInterfaceTypedDict
+
+
+class CUDAArrayInterface(Protocol):
+    __cuda_array_interface__: Dict[str, Any]
+
+
+PyCapsule = NewType("PyCapsule", object)
+DLPackDevice = Tuple[int, int]
+
+
+class DLPackBuffer(Protocol):
+    def __dlpack__(self, stream: Optional[int] = None) -> PyCapsule:
+        ...
+
+    def __dlpack_device__(self) -> DLPackDevice:
+        ...
+
+
 class DevToolsPrettyPrintable(Protocol):
     """Used by python-devtools (https://python-devtools.helpmanual.io/)."""
 
diff --git a/src/gt4py/next/iterator/ir.py b/src/gt4py/next/iterator/ir.py
index 5c1efca8..535648cc 100644
--- a/src/gt4py/next/iterator/ir.py
+++ b/src/gt4py/next/iterator/ir.py
@@ -147,7 +147,6 @@ BINARY_MATH_NUMBER_BUILTINS = {
     "minimum",
     "maximum",
     "fmod",
-    "power",
     "plus",
     "minus",
     "multiplies",
@@ -169,6 +168,17 @@ BINARY_LOGICAL_BUILTINS = {
     "xor_",
 }
 
+ARITHMETIC_BUILTINS = {
+    *UNARY_MATH_NUMBER_BUILTINS,
+    *UNARY_LOGICAL_BUILTINS,
+    *UNARY_MATH_FP_BUILTINS,
+    *UNARY_MATH_FP_PREDICATE_BUILTINS,
+    *BINARY_MATH_NUMBER_BUILTINS,
+    "power",
+    *BINARY_MATH_COMPARISON_BUILTINS,
+    *BINARY_LOGICAL_BUILTINS,
+}
+
 #: builtin / dtype used to construct integer indices, like domain bounds
 INTEGER_INDEX_BUILTIN = "int32"
 INTEGER_BUILTINS = {"int32", "int64"}
@@ -197,13 +207,7 @@ BUILTINS = {
     "can_deref",
     "scan",
     "if_",
-    *UNARY_MATH_NUMBER_BUILTINS,
-    *UNARY_LOGICAL_BUILTINS,
-    *UNARY_MATH_FP_BUILTINS,
-    *UNARY_MATH_FP_PREDICATE_BUILTINS,
-    *BINARY_MATH_NUMBER_BUILTINS,
-    *BINARY_MATH_COMPARISON_BUILTINS,
-    *BINARY_LOGICAL_BUILTINS,
+    *ARITHMETIC_BUILTINS,
     *TYPEBUILTINS,
 }
 
diff --git a/src/gt4py/next/iterator/type_inference.py b/src/gt4py/next/iterator/type_inference.py
index 13f7fb50..12118fbb 100644
--- a/src/gt4py/next/iterator/type_inference.py
+++ b/src/gt4py/next/iterator/type_inference.py
@@ -39,6 +39,17 @@ TYPED_IR_NODES: typing.Final = (
 )
 
 
+class UnsatisfiableConstraintsError(Exception):
+    unsatisfiable_constraints: list[tuple[Type, Type]]
+
+    def __init__(self, unsatisfiable_constraints: list[tuple[Type, Type]]):
+        self.unsatisfiable_constraints = unsatisfiable_constraints
+        msg = "Type inference failed: Can not satisfy constraints:"
+        for lhs, rhs in unsatisfiable_constraints:
+            msg += f"\n  {lhs} ≡ {rhs}"
+        super().__init__(msg)
+
+
 class EmptyTuple(Type):
     def __iter__(self) -> abc.Iterator[Type]:
         return
@@ -396,6 +407,12 @@ BUILTIN_CATEGORY_MAPPING = (
             ret=Val_T0_T1,
         ),
     ),
+    (
+        {"power"},
+        FunctionType(
+            args=Tuple.from_elems(Val_T0_T1, Val(kind=Value(), dtype=T2, size=T1)), ret=Val_T0_T1
+        ),
+    ),
     (
         ir.BINARY_MATH_NUMBER_BUILTINS,
         FunctionType(args=Tuple.from_elems(Val_T0_T1, Val_T0_T1), ret=Val_T0_T1),
@@ -582,6 +599,12 @@ class _TypeInferrer(eve.traits.VisitorWithSymbolTableTrait, eve.NodeTranslator):
         result = super().visit(node, **kwargs)
         if isinstance(node, TYPED_IR_NODES):
             assert isinstance(result, Type)
+            if not (
+                id(node) not in self.collected_types or self.collected_types[id(node)] == result
+            ):
+                # using the same node in multiple places is fine as long as the type is the same
+                # for all occurences
+                self.constraints.add((result, self.collected_types[id(node)]))
             self.collected_types[id(node)] = result
 
         return result
@@ -707,6 +730,9 @@ class _TypeInferrer(eve.traits.VisitorWithSymbolTableTrait, eve.NodeTranslator):
         if not (isinstance(node.args[0], ir.OffsetLiteral) and isinstance(node.args[0].value, str)):
             raise TypeError("The first argument to `neighbors` must be an `OffsetLiteral` tag.")
 
+        # Visit arguments such that their type is also inferred
+        self.visit(node.args, **kwargs)
+
         max_length: Type = TypeVar.fresh()
         has_skip_values: Type = TypeVar.fresh()
         if self.offset_provider:
@@ -820,7 +846,7 @@ class _TypeInferrer(eve.traits.VisitorWithSymbolTableTrait, eve.NodeTranslator):
         node: ir.FunCall,
         **kwargs,
     ) -> Type:
-        if isinstance(node.fun, ir.SymRef) and hasattr(self, f"_visit_{node.fun.id}"):
+        if isinstance(node.fun, ir.SymRef) and node.fun.id in ir.GRAMMAR_BUILTINS:
             # builtins that are treated as part of the grammar are handled in `_visit_<builtin_name>`
             return getattr(self, f"_visit_{node.fun.id}")(node, **kwargs)
 
@@ -963,10 +989,14 @@ def infer_all(
     collected_types = dict(reversed(inferrer.collected_types.items()))
 
     # Compute the most general type that satisfies all constraints
-    unified_types = unify(list(collected_types.values()), inferrer.constraints)
+    unified_types, unsatisfiable_constraints = unify(
+        list(collected_types.values()), inferrer.constraints
+    )
 
     if reindex:
-        unified_types = reindex_vars(list(unified_types))
+        unified_types, unsatisfiable_constraints = reindex_vars(
+            (unified_types, unsatisfiable_constraints)
+        )
 
     result = {
         id_: unified_type
@@ -976,6 +1006,9 @@ def infer_all(
     if save_to_annex:
         _save_types_to_annex(node, result)
 
+    if unsatisfiable_constraints:
+        raise UnsatisfiableConstraintsError(unsatisfiable_constraints)
+
     return result
 
 
diff --git a/src/gt4py/next/program_processors/codegens/gtfn/gtfn_ir.py b/src/gt4py/next/program_processors/codegens/gtfn/gtfn_ir.py
index 0f55af17..4bb54aad 100644
--- a/src/gt4py/next/program_processors/codegens/gtfn/gtfn_ir.py
+++ b/src/gt4py/next/program_processors/codegens/gtfn/gtfn_ir.py
@@ -144,24 +144,12 @@ GTFN_BUILTINS = [
     "named_range",
     "reduce",
 ]
-UNARY_MATH_NUMBER_BUILTINS = itir.UNARY_MATH_NUMBER_BUILTINS
-UNARY_LOGICAL_BUILTINS = itir.UNARY_LOGICAL_BUILTINS
-UNARY_MATH_FP_BUILTINS = itir.UNARY_MATH_FP_BUILTINS
-UNARY_MATH_FP_PREDICATE_BUILTINS = itir.UNARY_MATH_FP_PREDICATE_BUILTINS
-BINARY_MATH_NUMBER_BUILTINS = itir.BINARY_MATH_NUMBER_BUILTINS
-BINARY_MATH_COMPARISON_BUILTINS = itir.BINARY_MATH_COMPARISON_BUILTINS
-BINARY_LOGICAL_BUILTINS = itir.BINARY_LOGICAL_BUILTINS
+ARITHMETIC_BUILTINS = itir.ARITHMETIC_BUILTINS
 TYPEBUILTINS = itir.TYPEBUILTINS
 
 BUILTINS = {
     *GTFN_BUILTINS,
-    *UNARY_MATH_NUMBER_BUILTINS,
-    *UNARY_LOGICAL_BUILTINS,
-    *UNARY_MATH_FP_BUILTINS,
-    *UNARY_MATH_FP_PREDICATE_BUILTINS,
-    *BINARY_MATH_NUMBER_BUILTINS,
-    *BINARY_MATH_COMPARISON_BUILTINS,
-    *BINARY_LOGICAL_BUILTINS,
+    *ARITHMETIC_BUILTINS,
     *TYPEBUILTINS,
 }
 
diff --git a/src/gt4py/next/type_inference.py b/src/gt4py/next/type_inference.py
index f06b26e4..9b5d9070 100644
--- a/src/gt4py/next/type_inference.py
+++ b/src/gt4py/next/type_inference.py
@@ -116,27 +116,14 @@ def freshen(dtypes: list[T] | T) -> list[T] | T:
     return [_TypeVarReindexer(indexer).visit(dtype, index_map=index_map) for dtype in dtypes]
 
 
-@typing.overload
-def reindex_vars(dtypes: list[T]) -> list[T]:
-    ...
-
-
-@typing.overload
-def reindex_vars(dtypes: T) -> T:
-    ...
-
-
-def reindex_vars(dtypes: list[T] | T) -> list[T] | T:
+def reindex_vars(dtypes: typing.Any) -> typing.Any:
     """Reindex all type variables, to have nice indices starting at zero."""
-    if not isinstance(dtypes, list):
-        assert isinstance(dtypes, Type)
-        return reindex_vars([dtypes])[0]
 
     def indexer(index_map: dict[int, int]) -> int:
         return len(index_map)
 
     index_map = dict[int, int]()
-    return [_TypeVarReindexer(indexer).visit(dtype, index_map=index_map) for dtype in dtypes]
+    return _TypeVarReindexer(indexer).visit(dtypes, index_map=index_map)
 
 
 class _FreeVariables(eve.NodeVisitor):
@@ -265,20 +252,27 @@ class _Unifier:
             self._renamer.register(s)
             self._renamer.register(t)
 
-    def unify(self) -> list[Type]:
+    def unify(self) -> tuple[list[Type] | Type, list[tuple[Type, Type]]]:
         """Run the unification."""
+        unsatisfiable_constraints = []
         while self._constraints:
             constraint = self._constraints.pop()
-            handled = self._handle_constraint(constraint)
-            if not handled:
-                # Try with swapped LHS and RHS
-                handled = self._handle_constraint(constraint[::-1])
+            try:
+                handled = self._handle_constraint(constraint)
+                if not handled:
+                    # Try with swapped LHS and RHS
+                    handled = self._handle_constraint(constraint[::-1])
+            except TypeError:
+                # custom constraint handler raised an error as constraint is not satisfiable
+                # (contrary to just not handled)
+                handled = False
+
             if not handled:
-                raise TypeError(
-                    f"Can not satisfy constraint: {constraint[0].value} ≡ {constraint[1].value}"
-                )
+                unsatisfiable_constraints.append((constraint[0].value, constraint[1].value))
 
-        return [dtype.value for dtype in self._dtypes]  # unbox
+        unboxed_dtypes = [dtype.value for dtype in self._dtypes]
+
+        return unboxed_dtypes, unsatisfiable_constraints
 
     def _rename(self, x: Type, y: Type) -> None:
         """Type renaming/replacement."""
@@ -329,19 +323,30 @@ class _Unifier:
 
 
 @typing.overload
-def unify(dtypes: list[Type], constraints: set[tuple[Type, Type]]) -> list[Type]:
+def unify(
+    dtypes: list[Type], constraints: set[tuple[Type, Type]]
+) -> tuple[list[Type], list[tuple[Type, Type]]]:
     ...
 
 
 @typing.overload
-def unify(dtypes: Type, constraints: set[tuple[Type, Type]]) -> Type:
+def unify(
+    dtypes: Type, constraints: set[tuple[Type, Type]]
+) -> tuple[Type, list[tuple[Type, Type]]]:
     ...
 
 
-def unify(dtypes: list[Type] | Type, constraints: set[tuple[Type, Type]]) -> list[Type] | Type:
-    """Unify all given constraints."""
+def unify(
+    dtypes: list[Type] | Type, constraints: set[tuple[Type, Type]]
+) -> tuple[list[Type] | Type, list[tuple[Type, Type]]]:
+    """
+    Unify all given constraints.
+
+    Returns the unified types and a list of unsatisfiable constraints.
+    """
     if isinstance(dtypes, Type):
-        return unify([dtypes], constraints)[0]
+        result_types, unsatisfiable_constraints = unify([dtypes], constraints)
+        return result_types[0], unsatisfiable_constraints
 
     # Deduplicate type nodes, this can speed up later things a bit
     memo = dict[Type, Type]()
diff --git a/src/gt4py/storage/__init__.py b/src/gt4py/storage/__init__.py
index d282a782..c62c9998 100644
--- a/src/gt4py/storage/__init__.py
+++ b/src/gt4py/storage/__init__.py
@@ -14,29 +14,22 @@
 
 # flake8: noqa: F401
 
-"""GridTools storages classes."""
+"""GridTools storages utilities."""
 
-from . import interface, layout
-from .interface import empty, from_array, full, ones, zeros  # noqa: F401
-from .layout import from_name, register
-
-
-try:
-    from .interface import dace_descriptor  # noqa: F401
-except ImportError:
-    pass
+from . import cartesian
+from .cartesian import layout
+from .cartesian.interface import empty, from_array, full, ones, zeros  # noqa: F401
+from .cartesian.layout import from_name, register
 
 
 __all__ = [
-    "interface",
+    "cartesian",
     "layout",
     "empty",
     "from_array",
     "full",
     "ones",
-    "zeros" "from_name",
+    "zeros",
+    "from_name",
     "register",
 ]
-
-if "dace_descriptor" in globals():
-    __all__ += ["dace_descriptor"]
diff --git a/src/gt4py/storage/allocators.py b/src/gt4py/storage/allocators.py
new file mode 100644
index 00000000..adc45efa
--- /dev/null
+++ b/src/gt4py/storage/allocators.py
@@ -0,0 +1,408 @@
+# 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 abc
+import collections.abc
+import dataclasses
+import functools
+import math
+import operator
+
+import numpy as np
+
+from gt4py._core import definitions as core_defs
+from gt4py.eve import extended_typing as xtyping
+from gt4py.eve.extended_typing import (
+    TYPE_CHECKING,
+    Any,
+    Generic,
+    NewType,
+    Optional,
+    Protocol,
+    Sequence,
+    Tuple,
+    TypeAlias,
+    TypeGuard,
+    TypeVar,
+    Union,
+    cast,
+)
+
+
+try:
+    import cupy as cp
+except ImportError:
+    cp = None
+
+
+_ScalarT = TypeVar("_ScalarT", bound=core_defs.Scalar)
+
+
+_NDBuffer: TypeAlias = Union[
+    xtyping.ArrayInterface,
+    xtyping.CUDAArrayInterface,
+    xtyping.DLPackBuffer,
+]
+
+
+#: Tuple of positive integers encoding a permutation of the dimensions.
+BufferLayoutMap = NewType("BufferLayoutMap", Sequence[core_defs.PositiveIntegral])
+
+
+def is_valid_layout_map(value: Sequence[Any]) -> TypeGuard[BufferLayoutMap]:
+    dims = list(range(len(value)))
+    return (
+        isinstance(value, collections.abc.Sequence)
+        and len(value) == len(set(value))
+        and all(i in dims for i in value)
+    )
+
+
+@dataclasses.dataclass(frozen=True)
+class TensorBuffer(Generic[core_defs.NDArrayObjectT, _ScalarT]):
+    """
+    N-dimensional (tensor-like) memory buffer.
+
+    The actual class of the stored buffer and ndarray instances is
+    represented in the `NDBufferT` parameter and might be any n-dimensional
+    buffer-like class with a compatible buffer interface (e.g. NumPy
+    or CuPy `ndarray`.)
+
+    Attributes:
+        buffer: Raw allocated buffer.
+        memory_address: Memory address of the buffer.
+        device: Device where the buffer is allocated.
+        dtype: Data type descriptor.
+        shape: Tuple with lengths of the corresponding tensor dimensions.
+        strides: Tuple with sizes (in bytes) of the steps in each dimension.
+        layout_map: Tuple with the order of the dimensions in the buffer.
+            layout_map[i] = j means that the i-th dimension of the tensor
+            corresponds to the j-th dimension in the buffer.
+        byte_offset: Offset (in bytes) from the beginning of the buffer to
+            the first valid element.
+        byte_alignment: Alignment (in bytes) of the first valid element.
+        aligned_index: N-dimensional index of the first aligned element.
+        ndarray: N-dimensional tensor view of the allocated buffer.
+    """
+
+    buffer: _NDBuffer = dataclasses.field(hash=False)
+    memory_address: int
+    device: core_defs.Device
+    dtype: core_defs.DType[_ScalarT]
+    shape: core_defs.TensorShape
+    strides: Tuple[int, ...]
+    layout_map: BufferLayoutMap
+    byte_offset: int
+    byte_alignment: int
+    aligned_index: Tuple[int, ...]
+    ndarray: core_defs.NDArrayObjectT = dataclasses.field(hash=False)
+
+    @property
+    def ndim(self):
+        """Order of the tensor (`len(tensor_buffer.shape)`)."""
+        return len(self.shape)
+
+    def __array__(self, dtype: Optional[np.dtype] = None) -> np.ndarray:
+        if not hasattr(self.ndarray, "__array__"):
+            raise TypeError("Cannot export tensor buffer as NumPy array.")
+
+        return self.ndarray.__array__(dtype=dtype)  # type: ignore[call-overload] # TODO(egparades): figure out the mypy fix
+
+    @property
+    def __cuda_array_interface__(self) -> xtyping.CUDAArrayInterfaceTypedDict:
+        if not hasattr(self.ndarray, "__cuda_array_interface__"):
+            raise TypeError("Cannot export tensor buffer to CUDA array interface.")
+        return self.ndarray.__cuda_array_interface__
+
+    def __dlpack__(self) -> xtyping.PyCapsule:
+        if not hasattr(self.ndarray, "__dlpack__"):
+            raise TypeError("Cannot export tensor buffer to DLPack.")
+        return self.ndarray.__dlpack__()
+
+    def __dlpack_device__(self) -> xtyping.DLPackDevice:
+        if not hasattr(self.ndarray, "__dlpack_device__"):
+            raise TypeError("Cannot extract DLPack device from tensor buffer.")
+        return self.ndarray.__dlpack_device__()
+
+
+class BufferAllocator(Protocol[core_defs.NDArrayObjectT]):
+    """Protocol for buffer allocators."""
+
+    @property
+    def device_type(self) -> core_defs.DeviceType:
+        ...
+
+    def allocate(
+        self,
+        shape: Sequence[core_defs.IntegralScalar],
+        dtype: core_defs.DType[_ScalarT],
+        layout_map: BufferLayoutMap,
+        device: core_defs.Device,
+        byte_alignment: int,
+        aligned_index: Optional[Sequence[int]] = None,
+    ) -> TensorBuffer[core_defs.NDArrayObjectT, _ScalarT]:
+        """
+        Allocate a TensorBuffer with the given shape, layout and alignment settings.
+
+        Args:
+            device: Device where the buffer is allocated.
+            dtype: Data type descriptor.
+            shape: Tensor dimensions.
+            layout_map: layout of the dimensions in the buffer.
+                layout_map[i] = j means that the i-th dimension of the tensor
+                corresponds to the j-th dimension of the buffer.
+            byte_alignment: Alignment (in bytes) of the first valid element.
+            aligned_index: N-dimensional index of the first aligned element.
+        """
+        ...
+
+
+@dataclasses.dataclass(frozen=True, init=False)
+class _BaseNDArrayBufferAllocator(abc.ABC, Generic[core_defs.NDArrayObjectT]):
+    """Base class for buffer allocators using NumPy-like modules."""
+
+    def allocate(
+        self,
+        shape: Sequence[core_defs.IntegralScalar],
+        dtype: core_defs.DType[_ScalarT],
+        layout_map: BufferLayoutMap,
+        device: core_defs.Device,
+        byte_alignment: int,
+        aligned_index: Optional[Sequence[int]] = None,
+    ) -> TensorBuffer[core_defs.NDArrayObjectT, _ScalarT]:
+        if not core_defs.is_valid_tensor_shape(shape):
+            raise ValueError(f"Invalid shape {shape}")
+        ndim = len(shape)
+        if len(layout_map) != ndim or not is_valid_layout_map(layout_map):
+            raise ValueError(f"Invalid layout_map {layout_map} for shape {shape}")
+
+        # Compute number of items inside an aligned block
+        item_size = dtype.byte_size
+        if math.gcd(byte_alignment, item_size) not in (byte_alignment, item_size):
+            raise ValueError(
+                f"Incompatible 'byte_alignment' ({byte_alignment}) and 'dtype' size ({item_size})"
+            )
+        items_per_aligned_block = (byte_alignment // item_size) or 1
+
+        # Compute the padding required in the contiguous dimension to get aligned blocks
+        dims_layout = [layout_map.index(i) for i in range(len(shape))]
+        padded_shape_lst = list(shape)
+        if ndim > 0:
+            padded_shape_lst[dims_layout[-1]] = (
+                math.ceil(shape[dims_layout[-1]] / items_per_aligned_block)
+                * items_per_aligned_block
+            )
+        padded_shape = tuple(padded_shape_lst)
+        assert core_defs.is_valid_tensor_shape(padded_shape)
+        total_length = item_size * functools.reduce(operator.mul, padded_shape, 1) + (
+            byte_alignment - 1
+        )  # the worst case for misalignment is `byte_alignment - 1`
+
+        # Compute strides
+        strides_lst = [item_size] * len(shape)
+        accumulator = item_size
+        for i in range(len(shape) - 2, -1, -1):
+            accumulator = strides_lst[dims_layout[i]] = (
+                accumulator * padded_shape[dims_layout[i + 1]]
+            )
+        strides = tuple(strides_lst)
+
+        # Allocate total size
+        buffer = self.raw_alloc(total_length, device)
+        memory_address = self.array_ns.byte_bounds(buffer)[0]
+
+        # Compute final byte offset to align the requested buffer index
+        aligned_index = tuple(aligned_index or ([0] * len(shape)))
+        aligned_index_offset = (
+            (
+                items_per_aligned_block
+                * (int(math.ceil(aligned_index[dims_layout[-1]] / items_per_aligned_block)))
+                - aligned_index[dims_layout[-1]]
+            )
+            * item_size
+            if ndim > 0
+            else 0
+        )
+
+        allocation_mismatch_offset = (
+            byte_alignment - memory_address % byte_alignment
+        ) % byte_alignment
+        byte_offset = (aligned_index_offset + allocation_mismatch_offset) % byte_alignment
+
+        # Create shaped view from buffer
+        ndarray = self.tensorize(
+            buffer, dtype, shape, padded_shape, item_size, strides, byte_offset
+        )
+
+        if device.device_type == core_defs.DeviceType.ROCM:
+            # until we can rely on dlpack
+            ndarray.__hip_array_interface__ = {  # type: ignore[attr-defined]
+                "shape": ndarray.shape,  # type: ignore[union-attr]
+                "typestr": ndarray.dtype.descr[0][1],  # type: ignore[union-attr]
+                "descr": ndarray.dtype.descr,  # type: ignore[union-attr]
+                "stream": 1,
+                "version": 3,
+                "strides": ndarray.strides,  # type: ignore[union-attr, attr-defined]
+                "data": (ndarray.data.ptr, False),  # type: ignore[union-attr, attr-defined]
+            }
+
+        return TensorBuffer(
+            buffer=buffer,
+            memory_address=memory_address,
+            device=device,
+            dtype=dtype,
+            shape=shape,
+            strides=strides,
+            layout_map=layout_map,
+            byte_offset=byte_offset,
+            byte_alignment=byte_alignment,
+            aligned_index=aligned_index,
+            ndarray=ndarray,
+        )
+
+    @property
+    @abc.abstractmethod
+    def array_ns(self) -> _NumPyLikeNamespace[core_defs.NDArrayObjectT]:
+        pass
+
+    @abc.abstractmethod
+    def raw_alloc(self, length: int, device: core_defs.Device) -> _NDBuffer:
+        pass
+
+    @abc.abstractmethod
+    def tensorize(
+        self,
+        buffer: _NDBuffer,
+        dtype: core_defs.DType[_ScalarT],
+        shape: core_defs.TensorShape,
+        allocated_shape: core_defs.TensorShape,
+        item_size: int,
+        strides: Sequence[int],
+        byte_offset: int,
+    ) -> core_defs.NDArrayObjectT:
+        pass
+
+
+if TYPE_CHECKING:
+
+    class _NumPyLikeNamespace(Protocol[core_defs.NDArrayObjectT]):
+        class _NumPyLibModule(Protocol):
+            class _NumPyLibStridesModule(Protocol):
+                def as_strided(
+                    self, ndarray: core_defs.NDArrayObjectT, **kwargs: Any
+                ) -> core_defs.NDArrayObjectT:
+                    ...
+
+            stride_tricks: _NumPyLibStridesModule
+
+        lib: _NumPyLibModule
+
+        def empty(self, shape: core_defs.TensorShape, dtype: np.dtype) -> core_defs.NDArrayObjectT:
+            ...
+
+        def byte_bounds(self, ndarray: _NDBuffer) -> tuple[int, int]:
+            ...
+
+
+@dataclasses.dataclass(frozen=True)
+class NumPyLikeArrayBufferAllocator(_BaseNDArrayBufferAllocator[core_defs.NDArrayObjectT]):
+    device_type: core_defs.DeviceType
+    array_ns_ref: _NumPyLikeNamespace[core_defs.NDArrayObjectT]
+
+    @property
+    def array_ns(self) -> _NumPyLikeNamespace[core_defs.NDArrayObjectT]:
+        return self.array_ns_ref
+
+    def raw_alloc(self, length: int, device: core_defs.Device) -> _NDBuffer:
+        if device.device_type != core_defs.DeviceType.CPU and device.device_id != 0:
+            raise ValueError(f"Unsupported device {device} for memory allocation")
+
+        shape = (length,)
+        assert core_defs.is_valid_tensor_shape(shape)  # for mypy
+        return cast(
+            _NDBuffer, self.array_ns.empty(shape=shape, dtype=np.dtype(np.uint8))
+        )  # TODO(havogt): figure out how we type this properly
+
+    def tensorize(
+        self,
+        buffer: _NDBuffer,
+        dtype: core_defs.DType[_ScalarT],
+        shape: core_defs.TensorShape,
+        allocated_shape: core_defs.TensorShape,
+        item_size: int,
+        strides: Sequence[int],
+        byte_offset: int,
+    ) -> core_defs.NDArrayObjectT:
+        aligned_buffer = buffer[byte_offset : byte_offset + math.prod(allocated_shape) * item_size]  # type: ignore[index] # TODO(egparedes): should we extend `_NDBuffer`s to cover __getitem__?
+        flat_ndarray = aligned_buffer.view(dtype=np.dtype(dtype))
+        tensor_view = self.array_ns.lib.stride_tricks.as_strided(
+            flat_ndarray, shape=allocated_shape, strides=strides
+        )
+        if len(shape) and shape != allocated_shape:
+            shape_slices = tuple(slice(0, s, None) for s in shape)
+            tensor_view = tensor_view[shape_slices]
+
+        return tensor_view
+
+
+#: Registry of allocators for each device type.
+device_allocators: dict[core_defs.DeviceType, BufferAllocator] = {}
+
+device_allocators[core_defs.DeviceType.CPU] = NumPyLikeArrayBufferAllocator(
+    device_type=core_defs.DeviceType.CPU,
+    array_ns_ref=cast(_NumPyLikeNamespace, np) if TYPE_CHECKING else np,
+)
+
+if cp:
+    device_allocators[core_defs.DeviceType.CUDA] = NumPyLikeArrayBufferAllocator(
+        device_type=core_defs.DeviceType.CUDA,
+        array_ns_ref=cp,
+    )
+    device_allocators[core_defs.DeviceType.ROCM] = NumPyLikeArrayBufferAllocator(
+        device_type=core_defs.DeviceType.ROCM,
+        array_ns_ref=cp,
+    )
+
+
+def allocate(
+    shape: Sequence[core_defs.IntegralScalar],
+    dtype: core_defs.DType[_ScalarT],
+    layout_map: BufferLayoutMap,
+    *,
+    byte_alignment: int,
+    aligned_index: Optional[Sequence[int]] = None,
+    device: Optional[core_defs.Device] = None,
+    allocator: Optional[BufferAllocator] = None,
+) -> TensorBuffer:
+    """Allocate a TensorBuffer with the given settings on the given device."""
+    if device is None and allocator is None:
+        raise ValueError("No 'device' or 'allocator' specified")
+    if device is None:
+        assert allocator is not None  # for mypy
+        device = core_defs.Device(allocator.device_type, 0)
+    assert device is not None  # for mypy
+    allocator = allocator or device_allocators[device.device_type]
+    if device.device_type != allocator.device_type:
+        raise ValueError(f"Device {device} and allocator {allocator} are incompatible")
+
+    return allocator.allocate(
+        shape=shape,
+        dtype=dtype,
+        layout_map=layout_map,
+        byte_alignment=byte_alignment,
+        aligned_index=aligned_index,
+        device=device,
+    )
diff --git a/src/gt4py/storage/cartesian/__init__.py b/src/gt4py/storage/cartesian/__init__.py
new file mode 100644
index 00000000..b3689ce7
--- /dev/null
+++ b/src/gt4py/storage/cartesian/__init__.py
@@ -0,0 +1,15 @@
+# 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
diff --git a/src/gt4py/storage/interface.py b/src/gt4py/storage/cartesian/interface.py
similarity index 74%
rename from src/gt4py/storage/interface.py
rename to src/gt4py/storage/cartesian/interface.py
index 938161e3..6e19b9d7 100644
--- a/src/gt4py/storage/interface.py
+++ b/src/gt4py/storage/cartesian/interface.py
@@ -12,11 +12,16 @@
 #
 # SPDX-License-Identifier: GPL-3.0-or-later
 
+from __future__ import annotations
+
 import numbers
-from typing import Any, Optional, Sequence, Union
+from typing import Any, Optional, Protocol, Sequence, Tuple, Union
 
 import numpy as np
 
+from gt4py.storage import allocators
+from gt4py.storage.cartesian import layout, utils as storage_utils
+
 
 try:
     import cupy as cp
@@ -34,14 +39,23 @@ else:
     ArrayLike = Any  # type: ignore[misc]  # assign multiple types in both branches
     DTypeLike = Any  # type: ignore[misc]  # assign multiple types in both branches
 
-from . import layout, utils as storage_utils
+
+# Protocols
+class GTDimsInterface(Protocol):
+    __gt_dims__: Tuple[str, ...]
+
+
+class GTOriginInterface(Protocol):
+    __gt_origin__: Tuple[int, ...]
 
 
-def _error_on_invalid_backend(backend):
+# Helper functions
+def _error_on_invalid_preset(backend):
     if backend not in layout.REGISTRY:
         raise RuntimeError(f"Storage preset '{backend}' is not registered.")
 
 
+# Public interface
 def empty(
     shape: Sequence[int],
     dtype: DTypeLike = np.float64,
@@ -83,7 +97,7 @@ def empty(
         ValueError
             If illegal or inconsistent arguments are specified.
     """
-    _error_on_invalid_backend(backend)
+    _error_on_invalid_preset(backend)
     storage_info = layout.from_name(backend)
     assert storage_info is not None
     if storage_info["device"] == "gpu":
@@ -95,13 +109,15 @@ def empty(
         aligned_index, shape, dtype, dimensions
     )
 
-    _error_on_invalid_backend(backend)
+    _error_on_invalid_preset(backend)
 
     alignment = storage_info["alignment"]
     layout_map = storage_info["layout_map"](dimensions)
+    assert allocators.is_valid_layout_map(layout_map)
 
     dtype = np.dtype(dtype)
     _, res = allocate_f(shape, layout_map, dtype, alignment * dtype.itemsize, aligned_index)
+
     return res
 
 
@@ -307,11 +323,9 @@ def from_array(
         ValueError
             If illegal or inconsistent arguments are specified.
     """
-    is_cupy_array = cp is not None and isinstance(data, cp.ndarray)
-    asarray = storage_utils.as_cupy if is_cupy_array else storage_utils.as_numpy
-    shape = asarray(data).shape
+    shape = storage_utils.asarray(data).shape
     if dtype is None:
-        dtype = asarray(data).dtype
+        dtype = storage_utils.asarray(data).dtype
     dtype = np.dtype(dtype)
     if dtype.shape:
         if dtype.shape and not shape[-dtype.ndim :] == dtype.shape:
@@ -325,87 +339,8 @@ def from_array(
         dimensions=dimensions,
     )
 
-    if cp is not None and isinstance(storage, cp.ndarray):
-        storage[...] = storage_utils.as_cupy(data)
-    else:
-        storage[...] = storage_utils.as_numpy(data)
+    layout_info = layout.from_name(backend)
+    assert layout_info is not None
+    storage[...] = storage_utils.asarray(data, device=layout_info["device"])
 
     return storage
-
-
-if dace is not None:
-
-    def dace_descriptor(
-        shape: Sequence[int],
-        dtype: DTypeLike = np.float64,
-        *,
-        backend: str,
-        aligned_index: Optional[Sequence[int]] = None,
-        dimensions: Optional[Sequence[str]] = None,
-    ) -> dace.data.Array:
-        """Return a DaCe data descriptor which describes performance-optimal strides and alignment.
-
-        Parameters
-        ----------
-            shape : `Sequence` of `int`
-                The shape of the resulting data descriptor.
-            dtype :  DTypeLike, optional
-                The dtype encoded in the resulting data descriptor.
-
-        Keyword Arguments
-        -----------------
-            backend : `str`
-                The target backend for which the described allocation is optimized.
-            aligned_index: `Sequence` of `int`, optional
-                Indicate the index of the resulting array that most commonly corresponds to the origin of the compute
-                domain. If not passed, it is aligned at the data origin.
-            dimensions: `Sequence` of `str`, optional
-                Indicate the semantic meaning of the dimensions in the data descriptor. Only used for determining
-                optimal strides, the information is not stored.
-
-        Returns
-        -------
-            DaCe data descriptor
-                With strides that encode padding and aligned to provide optimal performance for the given `backend` and
-                `aligned_index`
-
-        Raises
-        -------
-            TypeError
-                If arguments of an unexpected type are specified.
-            ValueError
-                If illegal or inconsistent arguments are specified.
-        """
-        aligned_index, shape, dtype, dimensions = storage_utils.normalize_storage_spec(
-            aligned_index, shape, dtype, dimensions
-        )
-        itemsize = dtype.itemsize
-        storage_info = layout.from_name(backend)
-        assert storage_info is not None
-        layout_map = storage_info["layout_map"](dimensions)
-
-        order_idx = storage_utils.idx_from_order([i for i in layout_map if i is not None])
-        padded_shape = storage_utils.compute_padded_shape(
-            shape, storage_info["alignment"], order_idx
-        )
-
-        strides = storage_utils.strides_from_padded_shape(padded_shape, order_idx, itemsize)
-
-        storage = (
-            dace.StorageType.GPU_Global
-            if storage_info["device"] == "gpu"
-            else dace.StorageType.CPU_Heap
-        )
-        start_offset = int(np.array([aligned_index]) @ np.array([strides]).T) // itemsize
-
-        total_size = int(int(np.array([shape]) @ np.array([strides]).T) // itemsize)
-
-        start_offset = start_offset % storage_info["alignment"]
-        return dace.data.Array(
-            shape=shape,
-            strides=[s // itemsize for s in strides],
-            dtype=dace.typeclass(str(dtype)),
-            storage=storage,
-            total_size=total_size,
-            start_offset=start_offset,
-        )
diff --git a/src/gt4py/storage/layout.py b/src/gt4py/storage/cartesian/layout.py
similarity index 100%
rename from src/gt4py/storage/layout.py
rename to src/gt4py/storage/cartesian/layout.py
diff --git a/src/gt4py/storage/utils.py b/src/gt4py/storage/cartesian/utils.py
similarity index 51%
rename from src/gt4py/storage/utils.py
rename to src/gt4py/storage/cartesian/utils.py
index a68fae1f..e6060328 100644
--- a/src/gt4py/storage/utils.py
+++ b/src/gt4py/storage/cartesian/utils.py
@@ -12,15 +12,20 @@
 #
 # SPDX-License-Identifier: GPL-3.0-or-later
 
+from __future__ import annotations
+
 import collections.abc
 import math
 import numbers
-from typing import Any, Dict, Iterable, Optional, Protocol, Sequence, Tuple, Union, cast
+from typing import Any, Literal, Optional, Sequence, Tuple, Union, cast
 
 import numpy as np
 import numpy.typing as npt
 
+from gt4py._core import definitions as core_defs
 from gt4py.cartesian import config as gt_config
+from gt4py.eve.extended_typing import ArrayInterface, CUDAArrayInterface
+from gt4py.storage import allocators
 
 
 if np.lib.NumpyVersion(np.__version__) >= "1.20.0":
@@ -30,32 +35,33 @@ else:
 
 try:
     import cupy as cp
-    from cupy.lib.stride_tricks import as_strided
 except ImportError:
     cp = None
 
 
-class ArrayInterfaceType(Protocol):
-    __array_interface__: Dict[str, Any]
-
-
-class CudaArrayInterfaceType(Protocol):
-    __cuda_array_interface__: Dict[str, Any]
-
+FieldLike = Union["cp.ndarray", np.ndarray, ArrayInterface, CUDAArrayInterface]
 
-class GtDimsInterface(Protocol):
-    __gt_dims__: Tuple[str, ...]
 
-
-class GtOriginInterface(Protocol):
-    __gt_origin__: Tuple[int, ...]
+def _idx_from_order(order):
+    return list(np.argsort(order))
 
 
-FieldLike = Union["cp.ndarray", np.ndarray, ArrayInterfaceType, CudaArrayInterfaceType]
+def _compute_padded_shape(shape, items_per_alignment, order_idx):
+    padded_shape = list(shape)
+    if len(order_idx) > 0:
+        padded_shape[order_idx[-1]] = int(
+            math.ceil(padded_shape[order_idx[-1]] / items_per_alignment) * items_per_alignment
+        )
+    return padded_shape
 
 
-def idx_from_order(order):
-    return list(np.argsort(order))
+def _strides_from_padded_shape(padded_size, order_idx, itemsize):
+    stride_accumulator = 1
+    strides = [0] * len(padded_size)
+    for idx in reversed(order_idx):
+        strides[idx] = stride_accumulator * itemsize
+        stride_accumulator = stride_accumulator * padded_size[idx]
+    return list(strides)
 
 
 def dimensions_to_mask(dimensions: Tuple[str, ...]) -> Tuple[bool, ...]:
@@ -148,162 +154,6 @@ def normalize_storage_spec(
     return aligned_index, shape, dtype, dimensions
 
 
-def compute_padded_shape(shape, items_per_alignment, order_idx):
-    padded_shape = list(shape)
-    if len(order_idx) > 0:
-        padded_shape[order_idx[-1]] = int(
-            math.ceil(padded_shape[order_idx[-1]] / items_per_alignment) * items_per_alignment
-        )
-    return padded_shape
-
-
-def strides_from_padded_shape(padded_size, order_idx, itemsize):
-    stride_accumulator = 1
-    strides = [0] * len(padded_size)
-    for idx in reversed(order_idx):
-        strides[idx] = stride_accumulator * itemsize
-        stride_accumulator = stride_accumulator * padded_size[idx]
-    return list(strides)
-
-
-def allocate(aligned_index, shape, layout_map, dtype, alignment_bytes, allocate_f):
-    dtype = np.dtype(dtype)
-    if not (alignment_bytes % dtype.itemsize) == 0:
-        raise ValueError("Alignment must be a multiple of byte-width of dtype.")
-    itemsize = dtype.itemsize
-    items_per_alignment = int(alignment_bytes / itemsize)
-
-    order_idx = idx_from_order([i for i in layout_map if i is not None])
-    padded_shape = compute_padded_shape(shape, items_per_alignment, order_idx)
-
-    strides = strides_from_padded_shape(padded_shape, order_idx, itemsize)
-    if len(order_idx) > 0:
-        halo_offset = (
-            int(math.ceil(aligned_index[order_idx[-1]] / items_per_alignment)) * items_per_alignment
-            - aligned_index[order_idx[-1]]
-        )
-    else:
-        halo_offset = 0
-
-    padded_size = int(np.prod(padded_shape))
-    buffer_size = padded_size + items_per_alignment - 1
-    array, raw_buffer = allocate_f(buffer_size, dtype=dtype)
-
-    allocation_mismatch = int((array.ctypes.data % alignment_bytes) / itemsize)
-
-    alignment_offset = (halo_offset - allocation_mismatch) % items_per_alignment
-
-    field = np.reshape(array[alignment_offset : alignment_offset + padded_size], padded_shape)
-    if field.ndim > 0:
-        field.strides = strides
-        field = field[tuple(slice(0, s, None) for s in shape)]
-    return raw_buffer, field
-
-
-def allocate_gpu(
-    shape: Sequence[int],
-    layout_map: Iterable[Optional[int]],
-    dtype: DTypeLike,
-    alignment_bytes: int,
-    aligned_index: Optional[Sequence[int]],
-) -> Tuple["cp.ndarray", "cp.ndarray"]:
-    dtype = np.dtype(dtype)
-    assert (
-        alignment_bytes % dtype.itemsize
-    ) == 0, "Alignment must be a multiple of byte-width of dtype."
-    itemsize = dtype.itemsize
-    items_per_alignment = int(alignment_bytes / itemsize)
-
-    order_idx = idx_from_order([i for i in layout_map if i is not None])
-    padded_shape = compute_padded_shape(shape, items_per_alignment, order_idx)
-
-    if aligned_index is None:
-        aligned_index = [0] * len(shape)
-
-    strides = strides_from_padded_shape(padded_shape, order_idx, itemsize)
-    if len(order_idx) > 0:
-        halo_offset = (
-            int(math.ceil(aligned_index[order_idx[-1]] / items_per_alignment)) * items_per_alignment
-            - aligned_index[order_idx[-1]]
-        )
-    else:
-        halo_offset = 0
-
-    padded_size = int(np.prod(padded_shape))
-    buffer_size = padded_size + items_per_alignment - 1
-
-    device_raw_buffer = cp.empty((buffer_size,), dtype=dtype)
-
-    allocation_mismatch = int((device_raw_buffer.data.ptr % alignment_bytes) / itemsize)
-    alignment_offset = (halo_offset - allocation_mismatch) % items_per_alignment
-
-    device_field = as_strided(
-        device_raw_buffer[alignment_offset : alignment_offset + padded_size],
-        shape=padded_shape,
-        strides=strides,
-    )
-    if device_field.ndim > 0:
-        device_field = device_field[tuple(slice(0, s, None) for s in shape)]
-
-    if gt_config.GT4PY_USE_HIP:
-        # HIP/ROCm lack support for __cuda_array_interface__
-        device_field.__hip_array_interface__ = {
-            "shape": padded_shape,
-            "typestr": device_field.dtype.descr[0][1],
-            "descr": device_field.dtype.descr,
-            "stream": 1,
-            "version": 3,
-            "strides": strides,
-            "data": (device_field.data.ptr, False),
-        }
-
-    return device_raw_buffer, device_field
-
-
-def allocate_cpu(
-    shape: Sequence[int],
-    layout_map: Iterable[Optional[int]],
-    dtype: DTypeLike,
-    alignment_bytes: int,
-    aligned_index: Optional[Sequence[int]],
-) -> Tuple[np.ndarray, np.ndarray]:
-    dtype = np.dtype(dtype)
-    assert (
-        alignment_bytes % dtype.itemsize
-    ) == 0, "Alignment must be a multiple of byte-width of dtype."
-    itemsize = dtype.itemsize
-    items_per_alignment = int(alignment_bytes / itemsize)
-
-    order_idx = idx_from_order([i for i in layout_map if i is not None])
-    padded_shape = compute_padded_shape(shape, items_per_alignment, order_idx)
-
-    if aligned_index is None:
-        aligned_index = [0] * len(shape)
-
-    strides = strides_from_padded_shape(padded_shape, order_idx, itemsize)
-    if len(order_idx) > 0:
-        halo_offset = (
-            int(math.ceil(aligned_index[order_idx[-1]] / items_per_alignment)) * items_per_alignment
-            - aligned_index[order_idx[-1]]
-        )
-    else:
-        halo_offset = 0
-
-    padded_size = int(np.prod(padded_shape))
-    buffer_size = padded_size + items_per_alignment - 1
-    raw_buffer = np.empty(buffer_size, dtype=dtype)
-
-    allocation_mismatch = int((raw_buffer.ctypes.data % alignment_bytes) / itemsize)
-
-    alignment_offset = (halo_offset - allocation_mismatch) % items_per_alignment
-
-    field = np.reshape(raw_buffer[alignment_offset : alignment_offset + padded_size], padded_shape)
-    if field.ndim > 0:
-        field.strides = strides
-        field = field[tuple(slice(0, s, None) for s in shape)]
-    return raw_buffer, field
-
-
 def cpu_copy(array: Union[np.ndarray, "cp.ndarray"]) -> np.ndarray:
     if cp is not None:
         # it's not clear from the documentation if cp.asnumpy guarantees a copy.
@@ -313,23 +163,71 @@ def cpu_copy(array: Union[np.ndarray, "cp.ndarray"]) -> np.ndarray:
         return np.array(array)
 
 
-def as_numpy(array: FieldLike) -> np.ndarray:
-    return np.asarray(array)
+def asarray(
+    array: FieldLike, *, device: Literal["cpu", "gpu", None] = None
+) -> np.ndarray | cp.ndarray:
+    if device == "gpu" or (not device and hasattr(array, "__cuda_array_interface__")):
+        return cp.asarray(array)
+    if device == "cpu" or (
+        not device and (hasattr(array, "__array_interface__") or hasattr(array, "__array__"))
+    ):
+        return np.asarray(array)
 
+    if device:
+        raise ValueError(f"Invalid device: {device!s}")
 
-def as_cupy(array: FieldLike) -> "cp.ndarray":
-    return cp.asarray(array)
+    raise TypeError(f"Cannot convert {type(array)} to ndarray")
 
 
-def get_dims(obj: Union[GtDimsInterface, npt.NDArray]) -> Optional[Tuple[str, ...]]:
+def get_dims(obj: Union[core_defs.GTDimsInterface, npt.NDArray]) -> Optional[Tuple[str, ...]]:
     dims = getattr(obj, "__gt_dims__", None)
     if dims is None:
         return dims
     return tuple(str(d) for d in dims)
 
 
-def get_origin(obj: Union[GtDimsInterface, npt.NDArray]) -> Optional[Tuple[int, ...]]:
+def get_origin(obj: Union[core_defs.GTDimsInterface, npt.NDArray]) -> Optional[Tuple[int, ...]]:
     origin = getattr(obj, "__gt_origin__", None)
     if origin is None:
         return origin
     return tuple(int(o) for o in origin)
+
+
+def allocate_cpu(
+    shape: Sequence[int],
+    layout_map: allocators.BufferLayoutMap,
+    dtype: DTypeLike,
+    alignment_bytes: int,
+    aligned_index: Optional[Sequence[int]],
+) -> Tuple[allocators._NDBuffer, np.ndarray]:
+    device = core_defs.Device(core_defs.DeviceType.CPU, 0)
+    buffer = allocators.allocate(
+        shape,
+        core_defs.dtype(dtype),
+        layout_map=layout_map,
+        device=device,
+        byte_alignment=alignment_bytes,
+        aligned_index=aligned_index,
+    )
+    return buffer.buffer, buffer.ndarray
+
+
+def allocate_gpu(
+    shape: Sequence[int],
+    layout_map: allocators.BufferLayoutMap,
+    dtype: DTypeLike,
+    alignment_bytes: int,
+    aligned_index: Optional[Sequence[int]],
+) -> Tuple["cp.ndarray", "cp.ndarray"]:
+    device = core_defs.Device(
+        core_defs.DeviceType.ROCM if gt_config.GT4PY_USE_HIP else core_defs.DeviceType.CUDA, 0
+    )
+    buffer = allocators.allocate(
+        shape,
+        core_defs.dtype(dtype),
+        layout_map=layout_map,
+        device=device,
+        byte_alignment=alignment_bytes,
+        aligned_index=aligned_index,
+    )
+    return buffer.buffer, buffer.ndarray
diff --git a/tox.ini b/tox.ini
index 730d8dca..e16aaff2 100644
--- a/tox.ini
+++ b/tox.ini
@@ -36,22 +36,15 @@ extras =
     formatting
     dace: dace
     cuda: cuda
-    cuda110: cuda110
-    cuda111: cuda111
-    cuda112: cuda112
-    cuda113: cuda113
-    cuda114: cuda114
-    cuda115: cuda115
-    cuda116: cuda116
-    cuda117: cuda117
     cuda11x: cuda11x
+    cuda12x: cuda12x
 package = wheel
 wheel_build_env = .pkg
 pass_env = NUM_PROCESSES
 set_env =
     PYTHONWARNINGS = {env:PYTHONWARNINGS:ignore:Support for `[tool.setuptools]` in `pyproject.toml` is still *beta*:UserWarning}
 
-[testenv:cartesian-py{38,39,310}-{internal,dace}-{cpu,cuda,cuda110,cuda111,cuda112,cuda113,cuda114,cuda115,cuda116,cuda117,cuda11x}]
+[testenv:cartesian-py{38,39,310}-{internal,dace}-{cpu,cuda,cuda11x,cuda12x}]
 description = Run 'gt4py.cartesian' tests
 pass_env = {[testenv]pass_env}, BOOST_ROOT, BOOST_HOME, CUDA_HOME, CUDA_PATH, CXX, CC, OPENMP_CPPFLAGS, OPENMP_LDFLAGS, PIP_USER, PYTHONUSERBASE
 allowlist_externals =
@@ -62,9 +55,9 @@ allowlist_externals =
     rm
 commands =
     internal-cpu: python -m pytest --cache-clear -v -n {env:NUM_PROCESSES:1} -m "not requires_gpu and not requires_dace" {posargs} tests{/}cartesian_tests
-    internal-{cuda,cuda110,cuda111,cuda112,cuda113,cuda114,cuda115,cuda116,cuda117,cuda11x}: python -m pytest --cache-clear -v -n {env:NUM_PROCESSES:1} -m "requires_gpu and not requires_dace" {posargs} tests{/}cartesian_tests
+    internal-{cuda,cuda11x,cuda12x}: python -m pytest --cache-clear -v -n {env:NUM_PROCESSES:1} -m "requires_gpu and not requires_dace" {posargs} tests{/}cartesian_tests
     dace-cpu: python -m pytest --cache-clear -v -n {env:NUM_PROCESSES:1} -m "not requires_gpu and requires_dace" {posargs} tests{/}cartesian_tests
-    dace-{cuda,cuda110,cuda111,cuda112,cuda113,cuda114,cuda115,cuda116,cuda117,cuda11x}: python -m pytest --cache-clear -v -n {env:NUM_PROCESSES:1} -m "requires_gpu and requires_dace" {posargs} tests{/}cartesian_tests
+    dace-{cuda,cuda11x,cuda12x}: python -m pytest --cache-clear -v -n {env:NUM_PROCESSES:1} -m "requires_gpu and requires_dace" {posargs} tests{/}cartesian_tests
     python -m pytest --doctest-modules --doctest-ignore-import-errors src{/}gt4py{/}cartesian
 # commands_pre =
 #     rm -Rf tests/_reports/coverage*
@@ -92,11 +85,11 @@ commands =
     atlas: python -m pytest --cache-clear -v -n {env:NUM_PROCESSES:1} -m "requires_atlas" {posargs} tests{/}next_tests
     pytest --doctest-modules src{/}gt4py{/}next
 
-[testenv:storage-py{38,39,310}-{internal,dace}-{cpu,cuda,cuda110,cuda111,cuda112,cuda113,cuda114,cuda115,cuda116,cuda117,cuda11x}]
+[testenv:storage-py{38,39,310}-{internal,dace}-{cpu,cuda,cuda11x,cuda12x}]
 description = Run 'gt4py.storage' tests
 commands =
     cpu: python -m pytest --cache-clear -v -n {env:NUM_PROCESSES:1} -m "not requires_gpu" {posargs} tests{/}storage_tests
-    {cuda,cuda110,cuda111,cuda112,cuda113,cuda114,cuda115,cuda116,cuda117,cuda11x}: python -m pytest --cache-clear -v -n {env:NUM_PROCESSES:1} -m "requires_gpu" {posargs} tests{/}storage_tests
+    {cuda,cuda11x,cuda12x}: python -m pytest --cache-clear -v -n {env:NUM_PROCESSES:1} -m "requires_gpu" {posargs} tests{/}storage_tests
     #pytest doctest-modules {posargs} src{/}gt4py{/}storage
 
 [testenv:linters-py{38,39,310}]


Test patch:
diff --git a/src/gt4py/cartesian/testing/suites.py b/src/gt4py/cartesian/testing/suites.py
index 7e7a78c0..99ad14f8 100644
--- a/src/gt4py/cartesian/testing/suites.py
+++ b/src/gt4py/cartesian/testing/suites.py
@@ -29,7 +29,7 @@ from gt4py.cartesian import gtscript, utils as gt_utils
 from gt4py.cartesian.definitions import AccessKind, FieldInfo
 from gt4py.cartesian.gtc.definitions import Boundary, CartesianSpace, Index, Shape
 from gt4py.cartesian.stencil_object import StencilObject
-from gt4py.storage import utils as storage_utils
+from gt4py.storage.cartesian import utils as storage_utils
 
 from .input_strategies import (
     SymbolKind,
diff --git a/tests/cartesian_tests/integration_tests/multi_feature_tests/test_code_generation.py b/tests/cartesian_tests/integration_tests/multi_feature_tests/test_code_generation.py
index 51329c09..ad5f6de4 100644
--- a/tests/cartesian_tests/integration_tests/multi_feature_tests/test_code_generation.py
+++ b/tests/cartesian_tests/integration_tests/multi_feature_tests/test_code_generation.py
@@ -30,7 +30,7 @@ from gt4py.cartesian.gtscript import (
     interval,
     region,
 )
-from gt4py.storage import utils as storage_utils
+from gt4py.storage.cartesian import utils as storage_utils
 
 from cartesian_tests.definitions import ALL_BACKENDS, CPU_BACKENDS
 from cartesian_tests.integration_tests.multi_feature_tests.stencil_definitions import (
diff --git a/tests/cartesian_tests/integration_tests/multi_feature_tests/test_dace_parsing.py b/tests/cartesian_tests/integration_tests/multi_feature_tests/test_dace_parsing.py
index 7f67bc13..7705b637 100644
--- a/tests/cartesian_tests/integration_tests/multi_feature_tests/test_dace_parsing.py
+++ b/tests/cartesian_tests/integration_tests/multi_feature_tests/test_dace_parsing.py
@@ -24,7 +24,7 @@ from gt4py import cartesian as gt4pyc, storage as gt_storage
 from gt4py.cartesian import gtscript
 from gt4py.cartesian.gtscript import PARALLEL, computation, interval
 from gt4py.cartesian.stencil_builder import StencilBuilder
-from gt4py.storage import utils as storage_utils
+from gt4py.storage.cartesian import utils as storage_utils
 
 from cartesian_tests.utils import OriginWrapper
 
diff --git a/tests/next_tests/unit_tests/iterator_tests/test_type_inference.py b/tests/next_tests/unit_tests/iterator_tests/test_type_inference.py
index 5d031308..4eadb8e4 100644
--- a/tests/next_tests/unit_tests/iterator_tests/test_type_inference.py
+++ b/tests/next_tests/unit_tests/iterator_tests/test_type_inference.py
@@ -19,6 +19,31 @@ from gt4py.next.iterator import ir, ir_makers as im, type_inference as ti
 from gt4py.next.iterator.runtime import CartesianAxis
 
 
+def test_unsatisfiable_constraints():
+    a = ir.Sym(id="a", dtype=("float32", False))
+    b = ir.Sym(id="b", dtype=("int32", False))
+
+    testee = im.lambda_(a, b)(im.plus("a", "b"))
+
+    # TODO(tehrengruber): For whatever reason the ordering in the error message is not
+    #  deterministic. Ignoring for now, as we want to refactor the type inference anyway.
+    expected_error = [
+        (
+            "Type inference failed: Can not satisfy constraints:\n"
+            "  Primitive(name='int32') ≡ Primitive(name='float32')"
+        ),
+        (
+            "Type inference failed: Can not satisfy constraints:\n"
+            "  Primitive(name='float32') ≡ Primitive(name='int32')"
+        ),
+    ]
+
+    try:
+        inferred = ti.infer(testee)
+    except ti.UnsatisfiableConstraintsError as e:
+        assert str(e) in expected_error
+
+
 def test_sym_ref():
     testee = ir.SymRef(id="x")
     expected = ti.TypeVar(idx=0)
@@ -120,6 +145,14 @@ def test_plus():
     assert ti.pformat(inferred) == "(T₀¹, T₀¹) → T₀¹"
 
 
+def test_power():
+    testee = im.call("power")(im.literal_from_value(1.0), im.literal_from_value(2))
+    expected = ti.Val(kind=ti.Value(), dtype=ti.Primitive(name="float64"), size=ti.TypeVar(idx=0))
+    inferred = ti.infer(testee)
+    assert inferred == expected
+    assert ti.pformat(inferred) == "float64⁰"
+
+
 def test_eq():
     testee = ir.SymRef(id="eq")
     t = ti.Val(kind=ti.Value(), dtype=ti.TypeVar(idx=0), size=ti.TypeVar(idx=1))
diff --git a/tests/next_tests/unit_tests/test_type_inference.py b/tests/next_tests/unit_tests/test_type_inference.py
index 4b4baf5b..74178e75 100644
--- a/tests/next_tests/unit_tests/test_type_inference.py
+++ b/tests/next_tests/unit_tests/test_type_inference.py
@@ -81,5 +81,5 @@ def test_custom_type_inference():
 
     expected = Fun(arg=Fun(arg=Basic(name="int"), ret=Basic(name="int")), ret=ti.TypeVar(idx=0))
 
-    actual = ti.reindex_vars(ti.unify(dtype, constraints))
+    actual = ti.reindex_vars(ti.unify(dtype, constraints)[0])
     assert actual == expected
diff --git a/tests/storage_tests/unit_tests/test_interface.py b/tests/storage_tests/unit_tests/test_interface.py
index b1b96373..9d68f13a 100644
--- a/tests/storage_tests/unit_tests/test_interface.py
+++ b/tests/storage_tests/unit_tests/test_interface.py
@@ -12,8 +12,6 @@
 #
 # SPDX-License-Identifier: GPL-3.0-or-later
 
-from types import SimpleNamespace
-
 import hypothesis as hyp
 import hypothesis.strategies as hyp_st
 import numpy as np
@@ -26,9 +24,8 @@ except ImportError:
     cp = None
 
 import gt4py
-import gt4py.cartesian.backend as gt_backend
 from gt4py.cartesian import gtscript
-from gt4py.storage.utils import allocate_cpu, allocate_gpu, normalize_storage_spec
+from gt4py.storage.cartesian.utils import allocate_cpu, allocate_gpu, normalize_storage_spec
 
 
 CPU_LAYOUTS = [
@@ -137,13 +134,10 @@ def test_allocate_cpu(param_dict):
 
     raw_buffer, field = allocate_cpu(shape, layout_map, dtype, alignment_bytes, aligned_index)
 
-    # check that field is a view of raw_buffer
-    assert field.base is raw_buffer
-
     # check that memory of field is contained in raw_buffer
     assert (
-        field.ctypes.data >= raw_buffer.ctypes.data
-        and field[-1:].ctypes.data <= raw_buffer[-1:].ctypes.data
+        np.byte_bounds(field)[0] >= np.byte_bounds(raw_buffer)[0]
+        and np.byte_bounds(field)[1] <= np.byte_bounds(raw_buffer)[1]
     )
 
     # check if the first compute-domain point in the last dimension is aligned for 100 random "columns"
@@ -188,6 +182,9 @@ def test_allocate_cpu(param_dict):
     assert field.shape == shape
 
 
+# TODO(egparedes): keep this and the previous tests separated just to differentiate
+# between gpu and cpu tests, but factorize all common code in a single utility
+# function which gets called from the test functions bodies.
 @pytest.mark.requires_gpu
 @hyp.given(param_dict=allocation_strategy())
 def test_allocate_gpu(param_dict):
@@ -365,16 +362,19 @@ def _create_ndarray_for_test_from_array(*, dtype, aligned_index, shape, backend)
     )
 
 
-@pytest.mark.parametrize(
-    "alloc_fun",
-    [
+@pytest.fixture(
+    params=[
         gt4py.storage.empty,
         gt4py.storage.ones,
         gt4py.storage.zeros,
         _create_ndarray_for_test_fill,
         _create_ndarray_for_test_from_array,
-    ],
+    ]
 )
+def alloc_fun(request):
+    return request.param
+
+
 @pytest.mark.parametrize("backend", CPU_LAYOUTS)
 def test_cpu_constructor(alloc_fun, backend):
     stor = alloc_fun(dtype=np.float64, aligned_index=(1, 2, 3), shape=(2, 4, 6), backend=backend)
@@ -382,16 +382,13 @@ def test_cpu_constructor(alloc_fun, backend):
     assert isinstance(stor, np.ndarray)
 
 
-@pytest.mark.parametrize(
-    "alloc_fun",
-    [
-        gt4py.storage.empty,
-        gt4py.storage.ones,
-        gt4py.storage.zeros,
-        _create_ndarray_for_test_fill,
-        _create_ndarray_for_test_from_array,
-    ],
-)
+@pytest.mark.parametrize("backend", CPU_LAYOUTS)
+def test_cpu_constructor_0d(alloc_fun, backend):
+    stor = alloc_fun(shape=(), dtype=np.float64, backend=backend, aligned_index=())
+    assert stor.shape == ()
+    assert isinstance(stor, np.ndarray)
+
+
 @pytest.mark.parametrize(
     "backend",
     GPU_LAYOUTS,
@@ -402,6 +399,13 @@ def test_gpu_constructor(alloc_fun, backend):
     assert isinstance(stor, cp.ndarray)
 
 
+@pytest.mark.parametrize("backend", GPU_LAYOUTS)
+def test_gpu_constructor_0d(alloc_fun, backend):
+    stor = alloc_fun(shape=(), dtype=np.float64, backend=backend, aligned_index=())
+    assert stor.shape == ()
+    assert isinstance(stor, cp.ndarray)
+
+
 @pytest.mark.requires_gpu
 @hyp.given(param_dict=dimensions_strategy())
 def test_masked_storage_gpu(param_dict):
@@ -450,84 +454,3 @@ def test_non_existing_backend():
             aligned_index=[0, 0, 0],
             dtype=(np.float64, (3,)),
         )
-
-
-@pytest.mark.skipif(dace is None, reason="Storage __descriptor__ depends on dace.")
-class TestDescriptor:
-    @staticmethod
-    def ravel_with_padding(array):
-        if hasattr(array, "__cuda_array_interface__"):
-            interface = dict(array.__cuda_array_interface__)
-        else:
-            interface = dict(array.__array_interface__)
-
-        assert interface.get("offset", 0) == 0 and interface["data"] is not None
-
-        total_size = int(
-            int(np.array([array.shape]) @ np.array([array.strides]).T) // array.itemsize
-        )
-        interface["shape"] = (total_size,)
-        interface["strides"] = (min(array.strides),)
-
-        if hasattr(array, "__cuda_array_interface__"):
-            return cp.asarray(SimpleNamespace(__cuda_array_interface__=interface))
-        else:
-            return np.asarray(SimpleNamespace(__array_interface__=interface))
-
-    @pytest.mark.parametrize(
-        "backend",
-        [
-            "dace:cpu",
-            pytest.param("dace:gpu", marks=[pytest.mark.requires_gpu]),
-        ],
-    )
-    def test_device(self, backend):
-        backend_cls = gt_backend.from_name(backend)
-        descriptor: dace.data.Array = gt4py.storage.dace_descriptor(
-            backend=backend,
-            shape=(3, 7, 13),
-            aligned_index=(1, 2, 3),
-            dtype=np.float64,
-        )
-        if backend_cls.storage_info["device"] == "gpu":
-            assert descriptor.storage == dace.StorageType.GPU_Global
-        else:
-            assert descriptor.storage == dace.StorageType.CPU_Heap
-
-    @pytest.mark.parametrize(
-        "backend",
-        [
-            "dace:cpu",
-            pytest.param("dace:gpu", marks=[pytest.mark.requires_gpu]),
-        ],
-    )
-    def test_start_offset(self, backend):
-        backend_cls = gt_backend.from_name(backend)
-        aligned_index = (1, 2, 3)
-        stor = gt4py.storage.ones(
-            backend=backend,
-            shape=(3, 7, 13),
-            aligned_index=aligned_index,
-            dtype=np.float64,
-        )
-        descriptor: dace.data.Array = gt4py.storage.dace_descriptor(
-            backend=backend,
-            shape=(3, 7, 13),
-            aligned_index=(1, 2, 3),
-            dtype=np.float64,
-        )
-        raveled = TestDescriptor.ravel_with_padding(stor)[descriptor.start_offset :]
-        if backend_cls.storage_info["device"] == "gpu":
-            assert raveled.data.ptr % (backend_cls.storage_info["alignment"] * stor.itemsize) == 0
-            assert (
-                backend_cls.storage_info["alignment"] == 1
-                or stor.data.ptr % (backend_cls.storage_info["alignment"] * stor.itemsize) != 0
-            )
-        else:
-            assert (
-                raveled.ctypes.data % (backend_cls.storage_info["alignment"] * stor.itemsize) == 0
-            )
-            assert (
-                backend_cls.storage_info["alignment"] == 1
-                or stor.ctypes.data % (backend_cls.storage_info["alignment"] * stor.itemsize) != 0
-            )
diff --git a/tests/storage_tests/unit_tests/test_layout.py b/tests/storage_tests/unit_tests/test_layout.py
index 62270f00..a5c50981 100644
--- a/tests/storage_tests/unit_tests/test_layout.py
+++ b/tests/storage_tests/unit_tests/test_layout.py
@@ -14,7 +14,7 @@
 
 import pytest
 
-from gt4py.storage import layout as gt_layout
+from gt4py.storage.cartesian import layout as gt_layout
 
 
 @pytest.mark.parametrize(