I am facing problem for computing the ray paths and L matrix in the 3D slowness-based tetrahedral mesh using SPM method.

[BhaskarIlla]

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][4*n:(4*n+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"csrmatrix{i}": matrix for i, matrix in enumerate(csr_matrix_list)} np.savez(filename, **csr_matrix_dict)

[bernard-giroux]

What is the problem? You get an error meeeage?

[BhaskarIlla]

Please see the attached image. The DSPM/FSM methods generate an accurate ray path distribution (in yellow).) But the number of ray paths was significantly reduced when I used the SPM method (in white). So, I am unable to generate a reliable L matrix and theoretical travel times.

[BhaskarIlla]

Regarding this, I did not get any error messages.

[bernard-giroux]

Are you sure that you have less ray paths? With the SPM, rays are passing through nodes, so it might look like you have less rays because, for coarse meshes the ray paths can be superimposed for most of the trajectory.

[BhaskarIlla]

Professor, thank you very much for the clarification, and you're right. Just now, I verified that the number of SPM method-related theoretical travel times matches the number of receivers. The RMS between the SPM and DSPM theoretical travel times is 0.15 sec.

[bernard-giroux]

The tt you get with the SPM are higher because the ray paths are longer due to the constraint that they pass through nodes. You can increase accuray by using the option tt_from_rp, but there would be no real gain over the DSPM & FSM.

[BhaskarIlla]

Professor, thank you very much for the valuable information, and I will go through it.