segfault in ScalarField3D_value

[filefolder]

Hi Malcom,

I had been experiencing this issue for quite some time, I suspect caused by an edge effect somewhere as it mostly happens with shallower models.

2023260T11:18:19::INFO::0000:: Tracing rays. (using within the context of pyvoro as usual!)
[gadi-cpu-clx-2555:2519722:0:2519722] Caught signal 11 (Segmentation fault: address not mapped to object at address 0x813db60)
==== backtrace (tid:2519722) ====
 0 0x0000000000012cf0 __funlockfile()  :0
 1 0x000000000004e41a __pyx_f_7pykonal_6fields_13ScalarField3D_value()  /home/158/rp6923/pykonal/pykonal/fields.cpp:10512
 2 0x0000000000057f26 __pyx_f_7pykonal_6fields_13ScalarField3D_trace_ray()  /home/158/rp6923/pykonal/pykonal/fields.cpp:9622
 3 0x000000000005a736 __pyx_pf_7pykonal_6fields_13ScalarField3D_4trace_ray()  /home/158/rp6923/pykonal/pykonal/fields.cpp:10149
 4 0x000000000005a736 __pyx_pw_7pykonal_6fields_13ScalarField3D_5trace_ray()  /home/158/rp6923/pykonal/pykonal/fields.cpp:10132

adding the last bit within the else clause in fields.pyx around line 500 seems to be a valid workaround... FWIW to anyone else and curious if this is a thing you've ran into yourself or might have a better way to address.

        for iax in range(3):
            if (
                (
                    point[iax] < self.cy_min_coords[iax]
                    or point[iax] > self.cy_max_coords[iax]
                )
                and not self.cy_iax_isperiodic[iax]
                and not self.cy_iax_isnull[iax]
            ):
                return (null)
            idx[iax]   = (point[iax] - self.cy_min_coords[iax]) / self.cy_node_intervals[iax]
            if self.cy_iax_isnull[iax]:
                ii[iax][0] = 0
                ii[iax][1] = 0
            else:
                ii[iax][0]  = <Py_ssize_t>idx[iax]
                ii[iax][1]  = <Py_ssize_t>(ii[iax][0]+1) % self.npts[iax]

                # **** double check still in bounds
                if ii[iax][0] < 0 or ii[iax][0] >= self.npts[iax] or ii[iax][1] < 0 or ii[iax][1] >= self.npts[iax]:
                    return null

[jobh]

Hi @filefolder , from inspecting the code you posted above my guess is that this can only happen if self.cy_iax_isperiodic[iax] is True; otherwise, it would either have returned early in the first if test (as out of bounds) or it would not have entered the else branch in the second.

And if it is periodic (and out of bounds), the assigment to ii[iax][0] should probably be wrapped into bounds, i.e.,

                ii[iax][0]  = <Py_ssize_t>idx[iax] % self.npts[iax]

Only guesswork from my side.

[filefolder]

Thank you... that's interesting. I am not quite sure why that parameter would be periodic necessarily, but the model is type pykonal.fields.ScalarField3D(coord_sys='spherical') so perhaps that is the default?

In practice, despite the spherical coordinate system I am dealing only with local areas, so I would never want any wrap-around.

Been a while since I thought about this issue. I tend to think my hack is probably not very sound, but I can say that I haven't had any issues since, for whatever that's worth.

[jobh]

For a spherical model the periodicity would be on the angle (e.g., 181° = -179°).

[filefolder]

Ah, a NZ model I was building could be affected then. Perhaps that was the issue?

I could try removing my hack and replacing

ii[iax][0] = <Py_ssize_t>idx[iax] with
ii[iax][0] = <Py_ssize_t>idx[iax] % self.npts[iax] as you suggest

at some point (but not any time soon unfortunately)

[jobh]

I'm not familiar enough with this code to say for sure, maybe @malcolmw ? A sanity check on the result might be to transform all input coordinates by 180° and compare, if that is an option.

[jobh]

After re-reading your comment, I was mistaken: If your area is local, there should be no periodicity. Sorry for adding confusion here.