Hi there,
I am interested in using this code. It looks awesome as it can be run with python. I looked at the two examples. I was wondering whether it can trace only refractions. Is it possible to trace reflections as well?

Thanks
Pasan

This has been implemented for 3D but not 2D

Bernard,

I get a ValueError when running your example 5 inside Spyder.

Traceback (most recent call last):

File ~\AppData\Local\Programs\Anaconda3\Lib\site-packages\spyder_kernels\py3compat.py:356 in compat_exec
exec(code, globals, locals)

File c:\users\bjorn\documents\projekt\ambiguity\example5.py:199
mesh = Mesh2d(nodes, triangles.astype(int), method='SPM', n_secondary=10, aniso='weakly_anelliptical')

File src\ttcrpy\tmesh.pyx:1228 in ttcrpy.tmesh.Mesh2d.cinit

ValueError: Buffer dtype mismatch, expected 'int64_t' but got 'long'

I am confused as Jupyter works fine???

Pylint also reports an Error

E0611(no-name-in-module) ... No name 'Mesh2d' in module 'tt.crpy.tmesh'

but that could be due to a missing __init__.py' in the directory scr.

Björn

Dear sir, I got an error when using function "data_kernel_straight_rays" of Mesh3d. I had read the help of "Mesh3d.data_kernel_straight_rays", and it did not tell me what type of data I should use if 'numpy.ndarray' was not allowed.
Here is the error:

TypeError Traceback (most recent call last)
in
----> 1 C = Mesh3d.data_kernel_straight_rays(src, rcv)

TypeError: descriptor 'data_kernel_straight_rays' for 'ttcrpy.tmesh.Mesh3d' objects doesn't apply to a 'numpy.ndarray' object

Here is the help of "Mesh3d.data_kernel_straight_rays":

Help on method_descriptor:

data_kernel_straight_rays(...)
data_kernel_straight_rays(Tx, Rx) -> L

Raytracing with straight rays in 3D

Parameters
----------
Tx : np.ndarray
    source coordinates, nTx by 3
        - 1st column contains X coordinates,
        - 2nd contains Y coordinates
        - 3rd contains Z coordinates
Rx : np.ndarray
    receiver coordinates, nTx by 3
        - 1st column contains X coordinates,
        - 2nd contains Y coordinates
        - 3rd contains Z coordinates
Returns
-------
L : scipy csr_matrix
    data kernel matrix (tt = L*slowness)

Note
----
Tx and Rx should contain the same number of rows, each row corresponding
to a source-receiver pair

Thank you very much!

I am writing to seek your expertise and assistance in resolving a problem I've encountered while working with ttcrpy.

The issue I'm facing pertains to the ray path going beyond the preset grid. After examining the code, specifically the "getTraveltimeFromRaypath" function, I have identified that the problem occurs when the algorithm tries to determine the next point in the ray path, and this next point, referred to as "curr_pt," falls outside the boundaries of the grid.

I have attached my code for reference. In this code, I have selected random source and receiver coordinates within the range of the grid space. However, I have observed that there is a probability of triggering this error when running the code. Strangely, when I change the selection method to 'DSPM,' this error is not reported.

I would greatly appreciate your insights and guidance in resolving this issue. If you could review my code and provide any suggestions or modifications that could help ensure the ray path remains within the preset grid, it would be immensely helpful.

Thank you in advance for your time and assistance. I look forward to hearing from you 。
demo.txt

Hi, have you ever needed to take gradient of the "raytrace" function? Can you please give some advice about taking the gradient? I am struggling to use JAX package for applying "automatic differentiation" to "raytrace" function. However, since the slowness model is not a variable of the raytrace function (>> grid.set_slowness(model) ), I am having some problems. Or is the numeric method the only way to take the derivative of the raytrace function? Thank you

One more question. I cannot find an answer in the documentation.

Is Z positive down or up?

I just saw this issue, I am interested in making a ray tracing for 3D seismic acquisition purposes, something like this, and could you give an example of the code, or the steps, to carry out the reflections, it is not for a complex model, my idea is to do it on a flat layered model.

Thanks in advance.

python version:3.11.3 ttcrpy version:1.2.1 numpy:1.24.3 scipy: 1.10.1

src: 50 3D coordinates ,rcv: 50 3D coordinates
They are both generated in the following way:
src = np.array([[4356.383252, 1003.381716, 500.0514687].
[4356.383252, 1003.381716, 500.0514687].
................. ,])
xmin, xmax = 4200.00, 4600.00
ymin, ymax = 600.00, 1200.00
zmin, zmax = 500.00, 1000.00
#Gridding parameters
XGridlength = 25.0
YGridlength = 40.0
ZGridlength = 50.0
x = np.arange(xmin, xmax, XGridlength)
y = np.arange(ymin, ymax, YGridlength)
z = np.arange(zmin, zmax, ZGridlength)
#grid
grid = rg.Grid3d(x, y, z, cell_slowness=False, method='SPM',n_threads=12)

I made sure I entered the slowness correctly
The problematic line of code is as follows:
tt, rays = grid.raytrace(src, rcv, return_rays=True)
I encounter a warning message: [WinError 10054] The remote host forced an existing connection to close.

Additionally, if I attempt to solve only for the time-to-parameter (tt) using the following line of code:
tt = grid.raytrace(src, rcv)
I encounter a warning message:
"This evaluation is really slow and this is expected. In this case, it's possible to silence this error by raising the timeout, setting the PYDEVD_WARN_EVALUATION_TIMEOUT environment variable to a bigger value.

• The evaluation may need other threads running while it's running. In this case, it's possible to set the PYDEVD_UNBLOCK_THREADS_TIMEOUT environment variable so that if after a given timeout an evaluation doesn't finish, other threads are unblocked, or you can manually resume all threads."
  Thank you for reading

Hi

I just installed ttcrpy using pip and when I try to run the simple examples notebook, I get an error while importing tmesh.

I hope someone can help me with this.

selected.xlsx
for help.txt
Dear Author,
Recently when I use the 'DSPM' method of 'ttcrpy' in ipynb to get the 'L', I often encountered the following Kernel' crash.Also, when I used the 'SPM' method, it was even worse and my computer would just seize up. I have tried reducing the amount of src's, it doesn't work.

The error message is as follows：DeprecationWarning: Conversion of an array with ndim > 0 to a scalar is deprecated, and will error in future. Ensure you extract a single element from your array before performing this operation. (Deprecated NumPy 1.25.)tt_line,M_line= grid.raytrace(src, rcv,slowness=s0_f, compute_L=True)

My code and files have been uploaded and I hope you take the time to answer questions.
Thank you very much.

I'm trying to install the python wrapper (Python 3.6.6, Cython 0.29, clang-1000.10.44.4 on Mac Mojave); compilation of the 3D mesh code fails with the following error.

Mesh3Dttcr.cpp:47:29: error: no matching constructor for initialization of
      'ttcr::mesh' (aka 'Grid3Dunfs<double, unsigned int>')
        mesh_instance = new mesh(no, tet, eps, maxit, refPts, 2, rp, nt);
                            ^    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
../ttcr/Grid3Dunfs.h:42:9: note: candidate constructor not viable: no known
      conversion from 'std::vector<sxyz<double> >' to 'const bool' for 5th
      argument
        Grid3Dunfs(const std::vector<sxyz<T1>>& no,
        ^
../ttcr/Grid3Dunfs.h:52:9: note: candidate constructor not viable: requires at
      least 10 arguments, but 8 were provided
        Grid3Dunfs(const std::vector<sxyz<T1>>& no,
        ^
../ttcr/Grid3Dunfs.h:40:11: note: candidate constructor (the implicit copy
      constructor) not viable: requires 1 argument, but 8 were provided
    class Grid3Dunfs : public Grid3Dun<T1,T2,Node3Dn<T1,T2>> {

Hello,

I am trying to use TTCR to get travel-time tables from source locations along a rectangular grid. I have been able to get raypaths and sensitivity kernel using the matlab codes provided. I am not too familiar with Python and C++; I am wondering if there is a way from Matlab to get extract that information.

I have the following issue with example2.py, example1.py runs fine. This seems to be a gmsh issue. I have the latest gmsh installed.

python example2.py
Info : Clearing all models and views...
Info : Done clearing all models and views
Info : Meshing 1D...
Info : [ 0%] Meshing curve 10 (Nurb)
Info : [ 20%] Meshing curve 11 (Line)
Info : [ 30%] Meshing curve 12 (Nurb)
Info : [ 50%] Meshing curve 13 (Line)
Info : [ 60%] Meshing curve 14 (Line)
Info : [ 80%] Meshing curve 15 (Line)
Info : [ 90%] Meshing curve 16 (Line)
Info : Done meshing 1D (Wall 0.00973982s, CPU 0.009607s)
Info : Meshing 2D...
Info : [ 0%] Meshing surface 31 (Plane, Frontal-Delaunay)
Info : [ 50%] Meshing surface 32 (Plane, Frontal-Delaunay)
Info : Done meshing 2D (Wall 0.0669393s, CPU 0.065739s)
Info : 2390 nodes 4929 elements
Traceback (most recent call last):
File "example2.py", line 94, in
node, _ = gmsh.model.mesh.getNode(tag)
ValueError: too many values to unpack (expected 2)

Hi, This is a discussion rather than an issue. Iw ant to know if it is possible to implement the code for circular geometry. Is any example available?

Any answers and suggestions would be appreciated.

Regards
Abhi

Hi,
I saw that ttcr can perform ray tracing in complex topography. I'm wondering if there is any example available to perform tomography to obtain the inverted velocity model? Furthermore, if there is any example for reflection, refraction and joint tomography? Kindly advise!

Looking forward for a response.

Thanks
Abhi

Dear Prof. Giroux,

I am facing the problem for computing the L matrix using SPM. I tried couple of approaches but not succeed. Please find the attached code and generated errors information below.

1. "ValueError: Slowness vector has wrong size" occurred if I keep the cell_slowness as a 'False'

2."NotImplementedError: compute_L not implemented for mesh with slowness defined in cells" occurred if I keep the cell_slowness as a 'True'

code:
import gmsh
import numpy as np
from ttcrpy.tmesh import Mesh3d,Mesh2d
from ttcrpy import tmesh
import vtk
from vtk.util.numpy_support import vtk_to_numpy
import math
from scipy.sparse import csr_matrix, save_npz

gmsh.initialize()
mesh_file = "smesh_3d_new.msh"
gmsh.merge(mesh_file)
gmsh.model.mesh.renumberNodes()
gmsh.model.mesh.renumberElements()
num_elements = gmsh.model.mesh.getElements()
num_nodes = gmsh.model.mesh.getNodes()[1].shape[0]
physical_groups = gmsh.model.getPhysicalGroups()
for dim, tag in physical_groups:
name = gmsh.model.getPhysicalName(dim, tag)
print(f"Physical group (Dim {dim}, Tag {tag}): {name}")
entities = gmsh.model.getEntities()

slowness,mgrids=[],[]

for dim, tag in gmsh.model.getEntities():
elemTypes, elemTags, elemNodeTags = gmsh.model.mesh.getElements(dim, tag)
physicalTags = gmsh.model.getPhysicalGroupsForEntity(dim, tag)
if dim == 3:
elemTypes, elemTags, elemNodeTags = gmsh.model.mesh.getElements(dim, tag)
physicalTags = gmsh.model.getPhysicalGroupsForEntity(dim, tag)
name = gmsh.model.getPhysicalName(dim, physicalTags[0])
#print(name)
for n in range(len(elemTags[0])):
t = elemNodeTags[0][4n:(4n+4)]
mgrids.append(t)
if name == "volume7":
slowness.append(1/3.5)
elif name == "volume6":
slowness.append(1/3.65)
elif name == "volume5":
slowness.append(1/4.0)
elif name == "volume4":
slowness.append(1/4.63)
elif name == "volume3":
slowness.append(1/3.2)
elif name == "volume2":
slowness.append(1/3.9)
elif name == "volume1":
slowness.append(1/4.51)
print("shape\n")
slowness = np.array(slowness)
mgrids = np.array(mgrids)
uniqueTags = np.unique(mgrids)
equiv = np.empty((int(1+uniqueTags.max()),))
nodes = []
for n, tag in enumerate(uniqueTags):
equiv[tag] = n
node = gmsh.model.mesh.getNode(tag)
nodes.append(node[0])

print(mgrids.shape[0])
for n1 in range(mgrids.shape[0]):
for n2 in range(mgrids.shape[1]):
mgrids[n1, n2] = equiv[mgrids[n1, n2]]
nodes = np.array(nodes)
gmsh.finalize()

Source and receiver information

tx,ty,tz,rx,ry,rz,otime=np.genfromtxt("pick.dat",dtype=float,unpack=True,usecols=[1,3,2,4,6,5,7],delimiter=",") #z is lat,x is lon, and y is depth
#receivers
Rx=np.column_stack((rx,ry,rz))

source

Tx=np.column_stack((tx,ty,tz))

print("slowness shape ",np.shape(slowness))
print("nodes shape ",np.shape(nodes))
print("mgrids shape ",np.shape(mgrids))

#Error 1
mesh = Mesh3d(nodes, mgrids.astype(int), method='SPM',n_threads=8,cell_slowness=False,gradient_method=1,maxit=30,n_secondary=2,n_tertiary=2,translate_grid=False)
tt, L = mesh.raytrace(Tx,Rx,slowness=slowness_array,compute_L=True,return_rays=False)

#Error 2
mesh = Mesh3d(nodes, mgrids.astype(int), method='SPM',n_threads=8,cell_slowness=True,gradient_method=1,maxit=30,n_secondary=2,n_tertiary=2,translate_grid=False)
tt, L = mesh.raytrace(Tx,Rx,slowness=slowness_array,compute_L=True,return_rays=False)

Hi,

I was trying to run your example 1 and got an error message:

ModuleNotFoundError Traceback (most recent call last)
Cell In[1], line 9
6 plt.style.use('seaborn-paper')
7 plt.rc('font', **{'family': 'sans-serif', 'sans-serif': ['Helvetica']})
----> 9 from ttcrpy.rgrid import Grid3d

File src\ttcrpy\rgrid.pyx:28, in init ttcrpy.rgrid()

ModuleNotFoundError: No module named 'vtk'

I run the code on Jupyter Notebook, and version of ttcrypy is 1.2.2, numpy is 1.24.3.

Hi,

I've been using this package for a while now. Today I noticed an inconsistency between the results of Grid2d.raytrace. Here is an example, where the left figure shows the rays in case of using version 1.1.8 and the right figure shows the rays for version 1.2.1.

Here is the code for the slowness model and how I called Grid2d.raytrace

n_rows_cols = 8

xn=np.linspace(-150, 1750, n_rows_cols+1)
zn=xn.copy()

xfast = 1037.5
zfast = 1037.5
xslow = 562.5
zslow = 562.5

vbackground = 2900
vtrue = np.full((n_rows_cols, n_rows_cols), vbackground)
for i in range(n_rows_cols):
    for j in range(n_rows_cols):
        vtrue[i,j] *= (1 + 0.10 * np.exp(-.00004 * ((xn[i] - xfast)**2 + (zn[j] - zfast)**2)))
        vtrue[i,j] *= (1 - 0.15 * np.exp(-.00004 * ((xn[i] - xslow)**2 + (zn[j] - zslow)**2)))
mtrue = 1 / vtrue.reshape(-1,1)

grid = Grid2d(xn, zn)
tt, rays, J = grid.raytrace(sc, rc, slowness=mtrue.flatten(), 
                            return_rays=True, compute_L=True)

Respected Prof. Giroux,

I am working on 3D finite-element first-arrival tomography. For this research, I am using tetrahedral meshes with a 2 km size for the entire study region. This complex mesh generation has been generated in the GMSH software. I kept the source and receiver locations in these meshes, and I would like to compute the theoretical travel times, raypaths, and L matrix using this package. But I am facing problems with the 'SPM' method only, and the other methods ('FSM' and 'DSPM') are working fine.

Could you please suggest how to fix this problem?

Also, please see below the code.

import gmsh
import numpy as np
from ttcrpy.tmesh import Mesh3d,Mesh2d
from ttcrpy import tmesh
import vtk
from vtk.util.numpy_support import vtk_to_numpy
import math
from scipy.sparse import csr_matrix, save_npz

gmsh.initialize()
mesh_file = "smesh_3d_new.msh"
gmsh.merge(mesh_file)
gmsh.model.mesh.renumberNodes()
gmsh.model.mesh.renumberElements()
num_elements = gmsh.model.mesh.getElements()
num_nodes = gmsh.model.mesh.getNodes()[1].shape[0]
physical_groups = gmsh.model.getPhysicalGroups()
for dim, tag in physical_groups:
name = gmsh.model.getPhysicalName(dim, tag)
print(f"Physical group (Dim {dim}, Tag {tag}): {name}")
entities = gmsh.model.getEntities()

slowness,mgrids=[],[]

for dim, tag in gmsh.model.getEntities():
elemTypes, elemTags, elemNodeTags = gmsh.model.mesh.getElements(dim, tag)
physicalTags = gmsh.model.getPhysicalGroupsForEntity(dim, tag)
if dim == 3:
elemTypes, elemTags, elemNodeTags = gmsh.model.mesh.getElements(dim, tag)
physicalTags = gmsh.model.getPhysicalGroupsForEntity(dim, tag)
name = gmsh.model.getPhysicalName(dim, physicalTags[0])
#print(name)
for n in range(len(elemTags[0])):
t = elemNodeTags[0][4n:(4n+4)]
mgrids.append(t)
if name == "volume7":
slowness.append(1/3.5)
elif name == "volume6":
slowness.append(1/3.65)
elif name == "volume5":
slowness.append(1/4.0)
elif name == "volume4":
slowness.append(1/4.63)
elif name == "volume3":
slowness.append(1/3.2)
elif name == "volume2":
slowness.append(1/3.9)
elif name == "volume1":
slowness.append(1/4.51)
print("shape\n")
slowness = np.array(slowness)
mgrids = np.array(mgrids)
uniqueTags = np.unique(mgrids)
equiv = np.empty((int(1+uniqueTags.max()),))
nodes = []
for n, tag in enumerate(uniqueTags):
equiv[tag] = n
node = gmsh.model.mesh.getNode(tag)
nodes.append(node[0])

print(mgrids.shape[0])
for n1 in range(mgrids.shape[0]):
for n2 in range(mgrids.shape[1]):
mgrids[n1, n2] = equiv[mgrids[n1, n2]]
nodes = np.array(nodes)
gmsh.finalize()

Source and receiver information

tx,ty,tz,rx,ry,rz,otime=np.genfromtxt("pick.dat",dtype=float,unpack=True,usecols=[1,3,2,4,6,5,7],delimiter=",") #z-lat,x-lon, and y-depth
#receivers
Rx=np.column_stack((rx,ry,rz))

source

Tx=np.column_stack((tx,ty,tz))

#SPM method
mesh = Mesh3d(nodes, mgrids.astype(int), method='SPM',n_threads=10,cell_slowness=True,gradient_method=1,maxit=30,n_secondary=2,n_tertiary=2,translate_grid=False)
mesh.set_slowness(slowness)
mesh.to_vtk({'slowness': slowness}, 'smesh_3d_sl_spm_par')
ref=time.time()
tt, rays = mesh.raytrace(Tx, Rx, return_rays=True)
print(tt)
print("run time \t:",time.time()-ref)
mesh.to_vtk({"ray paths":rays},"smesh_3d_sl_rays_fsm_par")
RMSE = math.sqrt(np.square(np.subtract(otime,tt)).mean())
print("SPM \t",RMSE)

#save
np.savetxt("traveltimes_spm.dat",tt)
csr_matrix_list = [L]
filename = "L_matrix_spm.npz"
csr_matrix_dict = {f"csr_matrix_{i}": matrix for i, matrix in enumerate(csr_matrix_list)}
np.savez(filename, **csr_matrix_dict)

Hello,

I was wondering if you can/will add any functionality for irregular grid sizes to Grid3d. For example, I am working with data that has a finer vertical velocity which are homogeneous horizontally across the layers. The 1-D vertical spacing starts out at 10 (meters), which increases with depth until there is 200 meter spacing. For the horizontal spacing, I would like to set it at a higher interval, such as 50 meters. These nodes reflect changes in speed within the water column, and I am computing the travel-time of acoustic waves from a transducer to a transponder.

When trying to use an irregular grid with the current code I get this message: "FSM: Grid cells must be cubic".

I wasn't sure if this would be a more suitable problem to the Mesh3d module, but the example provided does not seem to work with the current version of pygmsh, and I am at a bit of a loss for where to start...

Hi,
I wish to use this code. But When I compile it, I got this error:

ld: library not found for -lvtkCommon
clang: error: linker command failed with exit code 1 (use -v to see invocation)
make: *** [ttcr3d] Error 1

I have installed VTK, eigen and boost. Below is my Makefile:

#!/bin/csh

CXX = g++
CXXOPTIONS = -O3 -DVTK -std=c++11
CXXFLAGS = $(CXXOPTIONS) -I/usr/local/include/vtk-8.1 -I/usr/local/include/eigen3 -I/usr/local/include/boost -I./ttcr
LFLAGS = -L/usr/local/lib
LIBS = -lvtkCommon -lvtkFilters -lvtkIO

ttcr3d: ./ttcr/ttcr3d.cpp ./ttcr/ttcr_io.h ./ttcr/structs_ttcr.h ./ttcr/utils.h ./ttcr/Rcv.h ./ttcr/Src.h ttcr_io.o
	$(CXX) $(CXXFLAGS) $(LFLAGS) $(LIBS) -Ittcr ttcr_io.o ./ttcr/ttcr3d.cpp -o ttcr3d

ttcr_io.o: ./ttcr/ttcr_io.h ./ttcr/ttcr_io.cpp ./ttcr/structs_ttcr.h
	$(CXX) $(CXXFLAGS) -c ./ttcr/ttcr_io.cpp

ttcr2d: ./ttcr/ttcr2d.cpp ./ttcr/ttcr_io.h ./ttcr/structs_ttcr.h ./ttcr/utils_ttcr.h ./ttcr/Rcv2D.h ./ttcr/Src2D.h ttcr_io.o
	$(CXX) $(CXXFLAGS) $(LFLAGS) $(LIBS) ttcr_io.o ./ttcr/ttcr2d.cpp -o ttcr2d

Hi!
I a graduate student at UT austin. I have couples of questions about this work.

1. Can it deal with a steep velocity change (such as 1.5μs/mm and 4μs/mm)? Reclinear grids with second order don't work in this circumstance. I saw your papers, your meshing can be refined at the boundary of two medium.
2. Does this matlab with c++ scripts are the same as python scripts?

I am trying to run example3.ipynb.

During ray tracing I get the following error. I have directly copied the notebook and have not changed anything.

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
/home/dmi/GIT/20220087-RUA/example3.ipynb Cell 5 in <cell line: 15>()
     [12] src = np.kron(src, np.ones((3, 1)))
     [13] rcv = np.kron(np.ones((2, 1)), rcv)
---> [15] tt, rays = grid.raytrace(src, rcv, return_rays=True)
     [17]) # Save raypaths
     [18]) grid.to_vtk({'raypaths for shot no 1': rays[:3],
     [19]             'raypaths for shot no 2': rays[3:]}, 'example3_rays')

File ttcrpy/rgrid.pyx:1064, in ttcrpy.rgrid.Grid3d.raytrace()

RuntimeError: Error while computing raypaths: going outside grid 
                Rx: 222.2 3.4 1.2
                Tx: 9 4.3 124.2

I ran into a similar error when doing the 2D version for my own problem. This is the first time I have used this code. Do you have any thoughts on the issue?

Hello, thank you for sharing this package.
I've installed numpy, python,vtk, cython. when I built the package, it failed because lack of definition of a macro. the terminal print is like this:
ttcrpy/Grid2Dttcr.cpp:535:9: error: use of undeclared identifier 'NUMPY_IMPORT_ARRAY_RETVAL'
import_array_throw(); // to use PyArray_SimpleNewFromData
I am fresh with python. What should I do?
Much appreciated!

Hi!

I got an error when trying your example. The error occurred when the jupyter notebook ran from ttcrpy.rgrid import Grid3d.

Here is the error:

---

ValueError Traceback (most recent call last)
/var/folders/3q/yps71nyj565b2g31vrgbwr50fztv4w/T/ipykernel_18831/2332482423.py in
7 plt.rc('font', **{'family': 'sans-serif', 'sans-serif': ['Helvetica']})
8
----> 9 from ttcrpy.rgrid import Grid3d

ttcrpy/rgrid.pyx in init ttcrpy.rgrid()

ValueError: numpy.ndarray size changed, may indicate binary incompatibility. Expected 88 from C header, got 80 from PyObject

Hope to hear from you soon. Thank you!

Dear author,

I encountered an error and crash when using cell_slowness in RayInvr grid computation. My grid settings are as follows:
########
python

#grid setting
XGridlength=500
YGridlength=700
ZGridlength=10
x = np.arange(0.00, 5000.00+XGridlength,XGridlength)
y = np.arange(0.00, 7000.00+YGridlength,YGridlength)
z = np.arange(50.00, 850.00+ZGridlength,ZGridlength)
#method setting
cpu_threads = multiprocessing.cpu_count()
grid=rg.Grid3d(x,y,z,cell_slowness=True,method='SPM',n_threads=cpu_threads)
grid.set_slowness(1.0/V)
tt, rays = grid.raytrace(src, rcv, return_rays=True)
########
The velocity model V has the same shape (10,10,80) as the grid. But this crash only happens when using cell_slowness. Node slowness works fine under the same condition.

I have tried checking the validity of V, reducing grid size, using single thread, but the problem remains. Please advise what could be the potential causes and how I can fix this. Any suggestions would be greatly appreciated.

Let me know if any other details are needed to reproduce the problem. Thanks for the great library!

Best regards,