search

SciTools-incubator / iris-esmf-regrid

A collection of structured and unstructured ESMF regridding schemes for Iris.
BSD 3-Clause "New" or "Revised" License
17 stars 16 forks source link

Segmentation fault when regridding a face located cube with the bilinear scheme #331

Closed jamessmith1metoffice closed 5 months ago

jamessmith1metoffice commented 6 months ago

Apologies in advance for the long listings. I thought it would be helpful to show my workings in full.

  1. I create a mesh programmatically:
def create_mesh():
    x_node_aux_coord = AuxCoord(
        points=[0.0, 5.0, 0.0, 5.0, 8.0],
        units="degrees",
        attributes={"location": Location.NODE},
        standard_name="longitude",
    )
    y_node_aux_coord = AuxCoord(
        points=[3.0, 3.0, 0.0, 0.0, 0.0],
        units="degrees",
        attributes={"location": Location.NODE},
        standard_name="latitude",
    )

    x_edge_aux_coord = AuxCoord(
        points=[2.5, 0, 5, 6.5, 2.5, 6.5],
        units="degrees",
        attributes={"location": Location.EDGE},
        standard_name="longitude",
    )
    y_edge_aux_coord = AuxCoord(
        points=[3, 1.5, 1.5, 1.5, 0, 0],
        units="degrees",
        attributes={"location": Location.EDGE},
        standard_name="latitude",
    )

    x_face_aux_coord = AuxCoord(
        points=[2.0, 6.0],
        units="degrees",
        attributes={"location": Location.FACE},
        standard_name="longitude",
    )
    y_face_aux_coord = AuxCoord(
        points=[1.0, 1.0],
        units="degrees",
        attributes={"location": Location.FACE},
        standard_name="latitude",
    )

    edge_node_connectivity = Connectivity(
        indices=[[0, 1], [0, 2], [1, 3], [1, 4], [2, 3], [3, 4]],
        cf_role="edge_node_connectivity",
    )

    face_node_connectivity = Connectivity(
        indices=np.ma.masked_equal([[0, 1, 3, 2], [1, 4, 3, 999]], 999),
        cf_role="face_node_connectivity",
    )

    mesh = Mesh(
        long_name="mesh",
        connectivities=[face_node_connectivity, edge_node_connectivity],
        topology_dimension=2,
        node_coords_and_axes=[(x_node_aux_coord, "x"), (y_node_aux_coord, "y")],
        edge_coords_and_axes=[(x_edge_aux_coord, "x"), (y_edge_aux_coord, "y")],
        face_coords_and_axes=[(x_face_aux_coord, "x"), (y_face_aux_coord, "y")],
    )

    return mesh
  1. I use this mesh to create a cube:
def create_mesh_cube(name, location=Location.FACE):
    mesh = create_mesh()

    data = _create_data(mesh, location)

    height_dim_coords = DimCoord([0], "height")

    x_mesh_coords, y_mesh_coords = mesh.to_MeshCoords(location.value)

    mesh_cube = Cube(
        long_name=name,
        units="K",
        data=data,
        dim_coords_and_dims=[(height_dim_coords, 1)],
        aux_coords_and_dims=[(x_mesh_coords, 0), (y_mesh_coords, 0)],
    )

    return mesh_cube

Note that I will invoke this function with location set to Face.

  1. The data I create by querying the length of the auxiliary coordinates. Note that this is why I set attributes on them when I created the mesh:
def _create_data(mesh, location):
    for aux_coord in mesh.all_coords:
        if aux_coord.attributes["location"] == location:
            data_length = aux_coord.shape[0]
            break

    data_height = 1

    data_shape = (data_length, data_height)

    data = np.arange(np.prod(data_shape)).reshape(data_shape)

    return data

When I create this cube and invoke the regridder...

    mesh_regridder = MeshToGridESMFRegridder(
        source_mesh_cube, target_grid_cube, method=method
    )

...I get the following segmentation fault:

Thread 0x00007fa950f47700 (most recent call first):
  File "/home/h02/jsmith1/.conda/envs/stage_ngms_development/lib/python3.10/concurrent/futures/thread.py", line 81 in _worker
  File "/home/h02/jsmith1/.conda/envs/stage_ngms_development/lib/python3.10/threading.py", line 953 in run
  File "/home/h02/jsmith1/.conda/envs/stage_ngms_development/lib/python3.10/threading.py", line 1016 in _bootstrap_inner
  File "/home/h02/jsmith1/.conda/envs/stage_ngms_development/lib/python3.10/threading.py", line 973 in _bootstrap
  File "/usr/share/java/pycharm-community/plugins/python-ce/helpers/pydev/_pydev_bundle/pydev_monkey.py", line 795 in __call__

Thread 0x00007fa98506d700 (most recent call first):
  File "/home/h02/jsmith1/.conda/envs/stage_ngms_development/lib/python3.10/threading.py", line 324 in wait
  File "/home/h02/jsmith1/.conda/envs/stage_ngms_development/lib/python3.10/threading.py", line 607 in wait
  File "/usr/share/java/pycharm-community/plugins/python-ce/helpers/pydev/pydevd.py", line 150 in _on_run
  File "/usr/share/java/pycharm-community/plugins/python-ce/helpers/pydev/_pydevd_bundle/pydevd_comm.py", line 219 in run
  File "/home/h02/jsmith1/.conda/envs/stage_ngms_development/lib/python3.10/threading.py", line 1016 in _bootstrap_inner
  File "/home/h02/jsmith1/.conda/envs/stage_ngms_development/lib/python3.10/threading.py", line 973 in _bootstrap

Thread 0x00007fa98586e700 (most recent call first):
  File "/usr/share/java/pycharm-community/plugins/python-ce/helpers/pydev/_pydevd_bundle/pydevd_comm.py", line 293 in _on_run
  File "/usr/share/java/pycharm-community/plugins/python-ce/helpers/pydev/_pydevd_bundle/pydevd_comm.py", line 219 in run
  File "/home/h02/jsmith1/.conda/envs/stage_ngms_development/lib/python3.10/threading.py", line 1016 in _bootstrap_inner
  File "/home/h02/jsmith1/.conda/envs/stage_ngms_development/lib/python3.10/threading.py", line 973 in _bootstrap

Thread 0x00007fa98606f700 (most recent call first):
  File "/home/h02/jsmith1/.conda/envs/stage_ngms_development/lib/python3.10/threading.py", line 324 in wait
  File "/home/h02/jsmith1/.conda/envs/stage_ngms_development/lib/python3.10/queue.py", line 180 in get
  File "/usr/share/java/pycharm-community/plugins/python-ce/helpers/pydev/_pydevd_bundle/pydevd_comm.py", line 368 in _on_run
  File "/usr/share/java/pycharm-community/plugins/python-ce/helpers/pydev/_pydevd_bundle/pydevd_comm.py", line 219 in run
  File "/home/h02/jsmith1/.conda/envs/stage_ngms_development/lib/python3.10/threading.py", line 1016 in _bootstrap_inner
  File "/home/h02/jsmith1/.conda/envs/stage_ngms_development/lib/python3.10/threading.py", line 973 in _bootstrap

Current thread 0x00007fa98e7ce740 (most recent call first):
  File "/home/h02/jsmith1/.conda/envs/stage_ngms_development/lib/python3.10/site-packages/ESMF/interface/cbindings.py", line 2089 in ESMP_FieldRegridStore
  File "/home/h02/jsmith1/.conda/envs/stage_ngms_development/lib/python3.10/site-packages/ESMF/api/regrid.py", line 184 in __init__
  File "/home/h02/jsmith1/.conda/envs/stage_ngms_development/lib/python3.10/site-packages/ESMF/util/decorators.py", line 81 in new_func
  File "/home/h02/jsmith1/.conda/envs/stage_ngms_development/lib/python3.10/site-packages/esmf_regrid/esmf_regridder.py", line 21 in _get_regrid_weights_dict
  File "/home/h02/jsmith1/.conda/envs/stage_ngms_development/lib/python3.10/site-packages/esmf_regrid/esmf_regridder.py", line 107 in __init__
  File "/home/h02/jsmith1/.conda/envs/stage_ngms_development/lib/python3.10/site-packages/esmf_regrid/schemes.py", line 624 in _regrid_unstructured_to_rectilinear__prepare
  File "/home/h02/jsmith1/.conda/envs/stage_ngms_development/lib/python3.10/site-packages/esmf_regrid/schemes.py", line 1344 in __init__
  File "/home/h02/jsmith1/.conda/envs/stage_ngms_development/lib/python3.10/site-packages/esmf_regrid/experimental/unstructured_scheme.py", line 169 in __init__
  File "/home/h02/jsmith1/Projects/stage_ngms/lib/stage_ngms/plugins/regrid.py", line 51 in CreateMeshRegridder
  File "/net/home/h02/jsmith1/Projects/stage/lib/stage/plugins/core.py", line 180 in result
  File "/net/home/h02/jsmith1/Projects/stage/lib/stage/plugins/core.py", line 112 in operation
  File "/home/h02/jsmith1/Projects/stage_ngms/lib/stage_ngms/tests/plugins/regrid_scheme.py", line 31 in test_regrid_cube_with_face_location_and_bilinear_method_does_throw_exception
  File "/home/h02/jsmith1/.conda/envs/stage_ngms_development/lib/python3.10/site-packages/_pytest/python.py", line 194 in pytest_pyfunc_call
  File "/home/h02/jsmith1/.conda/envs/stage_ngms_development/lib/python3.10/site-packages/pluggy/_callers.py", line 39 in _multicall
  File "/home/h02/jsmith1/.conda/envs/stage_ngms_development/lib/python3.10/site-packages/pluggy/_manager.py", line 80 in _hookexec
  File "/home/h02/jsmith1/.conda/envs/stage_ngms_development/lib/python3.10/site-packages/pluggy/_hooks.py", line 265 in __call__
  File "/home/h02/jsmith1/.conda/envs/stage_ngms_development/lib/python3.10/site-packages/_pytest/python.py", line 1788 in runtest
  File "/home/h02/jsmith1/.conda/envs/stage_ngms_development/lib/python3.10/site-packages/_pytest/runner.py", line 169 in pytest_runtest_call
  File "/home/h02/jsmith1/.conda/envs/stage_ngms_development/lib/python3.10/site-packages/pluggy/_callers.py", line 39 in _multicall
  File "/home/h02/jsmith1/.conda/envs/stage_ngms_development/lib/python3.10/site-packages/pluggy/_manager.py", line 80 in _hookexec
  File "/home/h02/jsmith1/.conda/envs/stage_ngms_development/lib/python3.10/site-packages/pluggy/_hooks.py", line 265 in __call__
  File "/home/h02/jsmith1/.conda/envs/stage_ngms_development/lib/python3.10/site-packages/_pytest/runner.py", line 262 in <lambda>
  File "/home/h02/jsmith1/.conda/envs/stage_ngms_development/lib/python3.10/site-packages/_pytest/runner.py", line 341 in from_call
  File "/home/h02/jsmith1/.conda/envs/stage_ngms_development/lib/python3.10/site-packages/_pytest/runner.py", line 261 in call_runtest_hook
  File "/home/h02/jsmith1/.conda/envs/stage_ngms_development/lib/python3.10/site-packages/_pytest/runner.py", line 222 in call_and_report
  File "/home/h02/jsmith1/.conda/envs/stage_ngms_development/lib/python3.10/site-packages/_pytest/runner.py", line 133 in runtestprotocol
  File "/home/h02/jsmith1/.conda/envs/stage_ngms_development/lib/python3.10/site-packages/_pytest/runner.py", line 114 in pytest_runtest_protocol
  File "/home/h02/jsmith1/.conda/envs/stage_ngms_development/lib/python3.10/site-packages/pluggy/_callers.py", line 39 in _multicall
  File "/home/h02/jsmith1/.conda/envs/stage_ngms_development/lib/python3.10/site-packages/pluggy/_manager.py", line 80 in _hookexec
  File "/home/h02/jsmith1/.conda/envs/stage_ngms_development/lib/python3.10/site-packages/pluggy/_hooks.py", line 265 in __call__
  File "/home/h02/jsmith1/.conda/envs/stage_ngms_development/lib/python3.10/site-packages/_pytest/main.py", line 349 in pytest_runtestloop
  File "/home/h02/jsmith1/.conda/envs/stage_ngms_development/lib/python3.10/site-packages/pluggy/_callers.py", line 39 in _multicall
  File "/home/h02/jsmith1/.conda/envs/stage_ngms_development/lib/python3.10/site-packages/pluggy/_manager.py", line 80 in _hookexec
  File "/home/h02/jsmith1/.conda/envs/stage_ngms_development/lib/python3.10/site-packages/pluggy/_hooks.py", line 265 in __call__
  File "/home/h02/jsmith1/.conda/envs/stage_ngms_development/lib/python3.10/site-packages/_pytest/main.py", line 324 in _main
  File "/home/h02/jsmith1/.conda/envs/stage_ngms_development/lib/python3.10/site-packages/_pytest/main.py", line 270 in wrap_session
  File "/home/h02/jsmith1/.conda/envs/stage_ngms_development/lib/python3.10/site-packages/_pytest/main.py", line 317 in pytest_cmdline_main
  File "/home/h02/jsmith1/.conda/envs/stage_ngms_development/lib/python3.10/site-packages/pluggy/_callers.py", line 39 in _multicall
  File "/home/h02/jsmith1/.conda/envs/stage_ngms_development/lib/python3.10/site-packages/pluggy/_manager.py", line 80 in _hookexec
  File "/home/h02/jsmith1/.conda/envs/stage_ngms_development/lib/python3.10/site-packages/pluggy/_hooks.py", line 265 in __call__
  File "/home/h02/jsmith1/.conda/envs/stage_ngms_development/lib/python3.10/site-packages/_pytest/config/__init__.py", line 166 in main
  File "/usr/share/java/pycharm-community/plugins/python-ce/helpers/pycharm/_jb_pytest_runner.py", line 51 in <module>
  File "/usr/share/java/pycharm-community/plugins/python-ce/helpers/pydev/_pydev_imps/_pydev_execfile.py", line 18 in execfile
  File "/usr/share/java/pycharm-community/plugins/python-ce/helpers/pydev/pydevd.py", line 1496 in _exec
  File "/usr/share/java/pycharm-community/plugins/python-ce/helpers/pydev/pydevd.py", line 1489 in run
  File "/usr/share/java/pycharm-community/plugins/python-ce/helpers/pydev/pydevd.py", line 2177 in main
  File "/usr/share/java/pycharm-community/plugins/python-ce/helpers/pydev/pydevd.py", line 2195 in <module>

Extension modules: numpy.core._multiarray_umath, numpy.core._multiarray_tests, numpy.linalg._umath_linalg, numpy.fft._pocketfft_internal, numpy.random._common, numpy.random.bit_generator, numpy.random._bounded_integers, numpy.random._mt19937, numpy.random.mtrand, numpy.random._philox, numpy.random._pcg64, numpy.random._sfc64, numpy.random._generator, cftime._cftime, cf_units._udunits2, yaml._yaml, cytoolz.utils, cytoolz.itertoolz, cytoolz.functoolz, cytoolz.dicttoolz, cytoolz.recipes, xxhash._xxhash, psutil._psutil_linux, psutil._psutil_posix, markupsafe._speedups, scipy._lib._ccallback_c, scipy.sparse._sparsetools, _csparsetools, scipy.sparse._csparsetools, scipy.sparse.linalg._isolve._iterative, scipy.linalg._fblas, scipy.linalg._flapack, scipy.linalg.cython_lapack, scipy.linalg._cythonized_array_utils, scipy.linalg._solve_toeplitz, scipy.linalg._decomp_lu_cython, scipy.linalg._matfuncs_sqrtm_triu, scipy.linalg.cython_blas, scipy.linalg._matfuncs_expm, scipy.linalg._decomp_update, scipy.linalg._flinalg, scipy.sparse.linalg._dsolve._superlu, scipy.sparse.linalg._eigen.arpack._arpack, scipy.sparse.csgraph._tools, scipy.sparse.csgraph._shortest_path, scipy.sparse.csgraph._traversal, scipy.sparse.csgraph._min_spanning_tree, scipy.sparse.csgraph._flow, scipy.sparse.csgraph._matching, scipy.sparse.csgraph._reordering, scipy._lib._uarray._uarray, scipy.special._ufuncs_cxx, scipy.special._ufuncs, scipy.special._specfun, scipy.special._comb, scipy.special._ellip_harm_2, scipy.fftpack.convolve, scipy.interpolate._fitpack, scipy.interpolate.dfitpack, scipy.optimize._minpack2, scipy.optimize._group_columns, scipy._lib.messagestream, scipy.optimize._trlib._trlib, scipy.optimize._lbfgsb, _moduleTNC, scipy.optimize._moduleTNC, scipy.optimize._cobyla, scipy.optimize._slsqp, scipy.optimize._minpack, scipy.optimize._lsq.givens_elimination, scipy.optimize._zeros, scipy.optimize.__nnls, scipy.optimize._highs.cython.src._highs_wrapper, scipy.optimize._highs._highs_wrapper, scipy.optimize._highs.cython.src._highs_constants, scipy.optimize._highs._highs_constants, scipy.linalg._interpolative, scipy.optimize._bglu_dense, scipy.optimize._lsap, scipy.spatial._ckdtree, scipy.spatial._qhull, scipy.spatial._voronoi, scipy.spatial._distance_wrap, scipy.spatial._hausdorff, scipy.spatial.transform._rotation, scipy.optimize._direct, scipy.interpolate._bspl, scipy.interpolate._ppoly, scipy.interpolate.interpnd, scipy.interpolate._rbfinterp_pythran, scipy.interpolate._rgi_cython, scipy.ndimage._nd_image, _ni_label, scipy.ndimage._ni_label, scipy.integrate._odepack, scipy.integrate._quadpack, scipy.integrate._vode, scipy.integrate._dop, scipy.integrate._lsoda, scipy.special.cython_special, scipy.stats._stats, scipy.stats.beta_ufunc, scipy.stats._boost.beta_ufunc, scipy.stats.binom_ufunc, scipy.stats._boost.binom_ufunc, scipy.stats.nbinom_ufunc, scipy.stats._boost.nbinom_ufunc, scipy.stats.hypergeom_ufunc, scipy.stats._boost.hypergeom_ufunc, scipy.stats.ncf_ufunc, scipy.stats._boost.ncf_ufunc, scipy.stats.ncx2_ufunc, scipy.stats._boost.ncx2_ufunc, scipy.stats.nct_ufunc, scipy.stats._boost.nct_ufunc, scipy.stats.skewnorm_ufunc, scipy.stats._boost.skewnorm_ufunc, scipy.stats.invgauss_ufunc, scipy.stats._boost.invgauss_ufunc, scipy.stats._biasedurn, scipy.stats._levy_stable.levyst, scipy.stats._stats_pythran, scipy.stats._statlib, scipy.stats._sobol, scipy.stats._qmc_cy, scipy.stats._mvn, scipy.stats._rcont.rcont, shapely.lib, shapely._geos, shapely._geometry_helpers, pyproj._compat, pyproj._datadir, pyproj._network, pyproj._geod, pyproj.list, pyproj._crs, pyproj._transformer, pyproj.database, pyproj._sync, cartopy.trace, netCDF4._netCDF4, tornado.speedups, msgpack._cmsgpack, lz4._version, lz4.block._block (total: 145)

I have little idea of where to start to make sense of this. I am hoping it will make sense to you!

I would expect this either to throw a helpful error informing me that I cannot regrid a face-located cube with the bilinear scheme, or to run without exceptions.

The version is 0.9.0. I am running on Red Hat Enterprise Linux Client 7.9.

stephenworsley commented 6 months ago

I suspect this may be due to the fact that your mesh contains only two faces. From what I can tell about how bilinear regridding works in ESMF, ESMF will try and create "cells" out of the face centers based on how the faces are connected. The data on target points will then be a weighted average the data on the of the face centers which make up the nodes of that "cell". In order to make a geometrically proper cell, it seems like you would need at least three face centers, so I imagine this might be what is causing ESMF to segfault. If you were to create a simple mesh out of three faces, each of which was connected to each other, I suspect this might work.

jamessmith1metoffice commented 6 months ago

Okay, understood. I will try that. But ESMF should not segfault, no? Should I open an issue somewhere else, perhaps?

jamessmith1metoffice commented 6 months ago

I can confirm that it is my fault, so to speak. A cube with three faces works fine.

pp-mo commented 5 months ago

@scitools/peloton it sounds like the only way to "fix" this would be to handle within ESMF itself, or just possibly in ESMPy, since there is generally no way to catch this in Python. It seems like checking whether data is well-conditioned for this operation isn't really practical at the Python (ESMPy) level either. You could create a issue https://github.com/esmf-org/esmf

jamessmith1metoffice commented 5 months ago

I have created an issue here:

https://github.com/esmf-org/esmf/issues/215