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

[anntsay]

Thanks for letting us know. we will take a look and get back to you. Please note that because of backlogs, we may take a bit longer than usual to do so. Thanks for the patience.

[billsacks]

Do you get a back trace with your segmentation fault? In particular, can you tell if it's occurring in the creation of one of the Field objects? I'm wondering if it's similar to the error I just encountered and am looking into: #326 (comment)

[billsacks]

Do you get a back trace with your segmentation fault? In particular, can you tell if it's occurring in the creation of one of the Field objects? I'm wondering if it's similar to the error I just encountered and am looking into: #326 (comment)

Actually, I'm guessing your problem has a different cause, because the problem I encountered in #326 turned out not to be directly due to the creation of a Field object. The fix in 6d5fdfe does fix some Segfault issues with ESMPy, but at a glance at your code, it doesn't look like it would be likely to fix your issue. If it's easy for you to apply that patch and retry, you could do so, but I don't have high hopes.