Segfault with conservative regridding of tetraherdons for some(!) resolutions

[MuellerSeb]

Hey there,
I found another strange behavior when I was trying to regrid 3D unstructured meshes with conservative regridding using esmpy 8.7 from conda-forge on linux with python 3.10.

For some resolutions I get a segfault. Here is a minimal working example using PyVista grids:

import esmpy
import numpy as np
import pyvista as pv


def shp(id, dim=3):
    res = dim * [1]
    res[id] = -1
    return res


def gen_mesh(pv_mesh):
    m = esmpy.Mesh(parametric_dim=3, spatial_dim=3, coord_sys=esmpy.CoordSys.CART)
    m.add_nodes(
        node_count=pv_mesh.n_points,
        node_ids=np.arange(pv_mesh.n_points, dtype=int) + 1,
        node_coords=pv_mesh.points.ravel().astype(float),
        node_owners=np.zeros(pv_mesh.n_points, dtype=int),
    )
    m.add_elements(
        element_count=pv_mesh.n_cells,
        element_ids=np.arange(pv_mesh.n_cells, dtype=int) + 1,
        element_types=pv_mesh.celltypes,  # works with TET and HEX
        element_conn=pv_mesh.cell_connectivity,
        element_coords=pv_mesh.cell_centers().points.ravel().astype(float),
    )
    return m


def gen_grid(pv_grid):
    dims = np.array([d - 1 for d in pv_grid.dimensions], dtype=int)
    locs = [esmpy.StaggerLoc.CORNER_VFACE, esmpy.StaggerLoc.CENTER_VCENTER]
    g = esmpy.Grid(dims, staggerloc=locs, coord_sys=esmpy.CoordSys.CART)
    for i in range(3):
        ax = getattr(pv_grid, ["x", "y", "z"][i])
        grid_corner = g.get_coords(i, staggerloc=locs[0])
        grid_corner[...] = ax.reshape(*shp(i))
        grid_center = g.get_coords(i, staggerloc=locs[1])
        grid_center[...] = ((ax[1:] + ax[:-1]) / 2).reshape(*shp(i))
    return g


def add_data(pv_mesh):
    x, y, z = pv_mesh.cell_centers().points.T / 12
    pv_mesh.cell_data["test"] = np.sin(x) * np.sin(y) * np.sin(z)


# I/O resolutions
res1, res2 = 10, 10

box1 = pv.ImageData(dimensions=3 * (res1 + 1,), spacing=3 * (60 / res1,))
box1 = box1.to_tetrahedra(pass_cell_ids=False)
add_data(box1)

box2 = pv.ImageData(dimensions=3 * (res2 + 1,), spacing=3 * (60 / res2,))
box2 = box2.cast_to_rectilinear_grid()

m = gen_mesh(box1)
g = gen_grid(box2)

field1 = esmpy.Field(m, name="input", meshloc=esmpy.MeshLoc.ELEMENT)
field1.data[...] = box1.cell_data["test"]

field2 = esmpy.Field(g, name="output", staggerloc=esmpy.StaggerLoc.CENTER_VCENTER)
field2.data[...] = np.nan

regrid = esmpy.Regrid(field1, field2, regrid_method=esmpy.RegridMethod.CONSERVE)
regrid(field1, field2)

box2.cell_data["test"] = field2.data.ravel(order="F")
box1.plot()
box2.plot()

The important part is res1, res2 = 10, 10, which sets the input- and output resolution. For me, it only works when:

1. res1 <= res2 (strange, since I typically want a coarser output grid when using conservative regridding), and
2. if res1 and res2 have an integer ratio

So, res1, res2 = 10, 10 works fine, res1, res2 = 5, 10 and res1, res2 = 10, 20 also, but res1, res2 = 10, 5 or res1, res2 = 10, 12 just gives a segfault and nothing more.

Note that regrid_method=esmpy.RegridMethod.NEAREST_STOD works for all resolution combinations.

Do you have any idea what went wrong here?

Here the images for res1, res2 = 10, 10 as references:

Cheers,
Sebastian