[v4-dev] UXArray.Grid.get_spatial_hash segmentation fault with python 3.13

[fluidnumerics-joe]

Parcels version

v4-dev

Description

When working with uxarray datasets for parcels.Field the initialization fails during the creation of the Field within Fieldset.__init__.

This is related to https://github.com/UXARRAY/uxarray/issues/1212

This only has been observed with python v3.13 . Currently, the workaround is to use python v3.11 or v3.12

Code sample

def stommel_fieldset_uxarray(xdim=200, ydim=200):
    """Simulate a periodic current along a western boundary, with significantly
    larger velocities along the western edge than the rest of the region

    The original test description can be found in: N. Fabbroni, 2009,
    Numerical Simulation of Passive tracers dispersion in the sea,
    Ph.D. dissertation, University of Bologna
    http://amsdottorato.unibo.it/1733/1/Fabbroni_Nicoletta_Tesi.pdf
    """
    import math

    import numpy as np
    import pandas as pd
    import uxarray as ux

    a = b = 66666 * 1e3/ 1111111.111111111 # Convert to rough lat/lon
    scalefac = 0.00025  # to scale for physically meaningful velocities

    # Coordinates of the test fieldset
    # Crowd points to the west edge of the domain
    # using a polyonmial map on x-direction
    lon, lat = np.meshgrid(np.linspace(0, a, ydim, dtype=np.float32),
                           np.linspace(0, b, ydim, dtype=np.float32))
    points = (lon.flatten(), lat.flatten())

    # Create the grid
    uxgrid = ux.Grid.from_points(points, method="regional_delaunay")
    uxgrid.construct_face_centers()
    uxgrid.get_spatial_hash()

    # Define arrays U (zonal), V (meridional) and P (sea surface height)
    U = np.zeros((1, 1, lat.size), dtype=np.float32)
    V = np.zeros((1, 1, lat.size), dtype=np.float32)
    P = np.zeros((1, 1, lat.size), dtype=np.float32)

    beta = 2e-11
    r = 1 / (11.6 * 86400)
    es = r / (beta * a)

    i = 0
    for x, y in zip(lon.flatten(), lat.flatten()):
        xi = x / a
        yi = y / b
        P[0, 0, i] = (
            (1 - math.exp(-xi / es) - xi) * math.pi * np.sin(math.pi * yi) * scalefac
        )
        U[0, 0, i] = (
            -(1 - math.exp(-xi / es) - xi)
            * math.pi**2
            * np.cos(math.pi * yi)
            * scalefac
        )
        V[0, 0, i] = (
            (math.exp(-xi / es) / es - 1) * math.pi * np.sin(math.pi * yi) * scalefac
        )
        i += 1

    u = ux.UxDataArray(
        data=U,
        name="u",
        uxgrid=uxgrid,
        dims=["time", "nz1", "n_node"],
        coords=dict(
            time=(["time"], pd.to_datetime(["2000-01-01"])),
            nz1=(["nz1"], [0]),
        ),
        attrs=dict(
            description="zonal velocity",
            units="m/s",
            location="node",
            mesh="delaunay",
        ),
    )
    v = ux.UxDataArray(
        data=V,
        name="v",
        uxgrid=uxgrid,
        dims=["time", "nz1", "n_node"],
        coords=dict(
            time=(["time"], pd.to_datetime(["2000-01-01"])),
            nz1=(["nz1"], [0]),
        ),
        attrs=dict(
            description="meridional velocity",
            units="m/s",
            location="node",
            mesh="delaunay",
        ),
    )
    p = ux.UxDataArray(
        data=P,
        name="p",
        uxgrid=uxgrid,
        dims=["time", "nz1", "n_node"],
        coords=dict(
            time=(["time"], pd.to_datetime(["2000-01-01"])),
            nz1=(["nz1"], [0]),
        ),
        attrs=dict(
            description="pressure",
            units="N/m^2",
            location="node",
            mesh="delaunay",
        ),
    )

    return ux.UxDataset({"U": u, "V": v, "p": p}, uxgrid=uxgrid)


from parcels import Fieldset

uxds = stommel_fieldset_uxarray(10, 10)
fieldset = FieldSet([uxds])```


```text
Segmentation fault (core dumped)

[VeckoTheGecko]

Tests are passing in #2327 , and from our discussions @fluidnumerics-joe it sounds like this can be closed now?