Irregular grid size support

[SquirrelKnight]

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...

[bernard-giroux]

Hi,

You don't have the cubic cell restriction with the SPM or DSPM (although I did not test it much in terms of accuracy).

I think that pygmsh has been abandoned since python support has been incorporated in gmsh (I should update the example...).

[SquirrelKnight]

Hi @bernard-giroux . Thanks for the quick response. I tried out DSPM. It looks like I don't need cubic cells any longer, however it still doesn't work with differently spaced Z cells (e.g. 0, 5, 10, 20, 30, 50). The program assumes that the difference in Z is the difference between the first pair of Z values (in this case 5), and so raytracing for a point at a Z of 45 comes up with an 'outside grid' error. In this case, I believe the program would think the grid ends at 25. I'm not sure if this would be too complicated to implement, or even very applicable to other cases.

FSM is also significantly faster when comparing computations on identical grids, but I suppose that's to be expected.

[bernard-giroux]

Can you share code excerpts that allows reproducing this? I might have some time in the upcoming days to look at it.

[SquirrelKnight]

@bernard-giroux Here is the code (3d_raytrace_ttcrpy.txt) I have written to test with some of my data. I have also included the data files as attachments. Note that you will need Pandas for loading the CSV files with the way the code is currently setup.

The variables x_int, y_int, and z_int define the (meter) spacing between nodes in their respective directions. Note that the code is currently set to interpolate depth and velocities from a 1-D model for the Z-axis, but this can be changed in lines 21 and 22. In line 101, you can change the Grid3d method.

The rest of the velocity model dimensions are determined from the initcfg input file. The final line reports the travel-time difference between the observed and estimated travel times for the record being examined.

GNSSA-03.2022_1-obs_sub.csv GNSSA-03.2022_1-syn-svp.csv GNSSA-03.2022_1-initcfg.csv 3d_raytrace_ttcrpy.txt

[bernard-giroux]

Hi,

Actually, non constant spacing is not currently implemented for grids, sorry. What is possible is to use a different dy, dy, or dz... I'd like to implement non constant spacing, but it will take some time.

The short term solution is to use tetrahedral meshes. Example 2 has been updated to use gmsh, perhaps you can have a look at this one to see the sequence of commands to create the mesh.

[SquirrelKnight]

Hi @bernard-giroux. I have followed your advice and am currently working on adapting your code for example 2 to work with my velocity structure. Because I am using a 3d model, I have hit a bit of a snag. I understand we want to have the 4 nodes forming a tetrahedra, as opposed to the 3 nodes forming a triangle, which you have as:

        for n in range(len(elemTags[0])):
            t = elemNodeTags[0][3*n:(3*n+3)]
            triangles.append(t)

For this snippet of code, would you simply change it to this?

        for n in range(len(elemTags[0])):
            t = elemNodeTags[0][3*n:(3*n+4)]
            tetra.append(t)

or should the other 3s be replaced with 4s as well?

Thanks for you patience!

[bernard-giroux]

I think that the other 3s should be replaced as well, best way to know it to test it... I'm on vacation until July 25. I'll have more time when I get back. Cheers

[SquirrelKnight]

@bernard-giroux Thanks! I appreciate your time. I have figured out that the 4s are indeed the best way to do things. I am going to upload to you the script I have made so far; it uses the same input files as the previous example.

The takes the layers of my 3d velocity model and extrapolates them across 100x100 meters in X and Y directions. I use the mesh recombination algorithm from gmsh to create tetrahedra (as opposed to triangles), but I'm not sure if this step is even really necessary. Whether perform recombination or not, the program seems to hang indefinitely on the raypath step, so I am likely still doing something wrong. I have noticed that the mesh.get_grid_traveltimes() does not seem to work with Mesh3d, even in much simpler velocity models

Have a good vacation; I'll be doing so myself right about when you are getting back! ray_gmsh.txt

[SquirrelKnight]

@bernard-giroux I've just started looking into things again. It looks like a major issue is created by building the mesh from several rectangles (in 2D) or cubes (in 3D).

For example, my code uses the addCube function for each layer of the velocity model. When raytracing within the first velocity layer, I have no difficulty getting an output traveltime or raypath. When setting either the source or the receiver in within a depth range of the second layer, the program either hangs or crashes.

I have subsequently created a singular cube model, that linearly interpolates velocity for each node based on the depth. Using this method, I have not come across any issues, so I have some nicely functional code, at the moment, and can compute the two-way travel-time for over 10,000 points in about a minute.

So, I have identified the source of the problem, but not really why. I'm a bit perplexed because the input to Mesh2d or 3d simply uses the 2d numpy arrays. I've checked to see that the triangle values correspond to node indices, and there doesn't seem to be any problem there. Could the source of the problem have to do with abrupt reconfigurations of the triangles for each layer within the grid?

What's your opinion?

[bernard-giroux]

@SquirrelKnight Sorry for the long silence. It's hard to see what could be the problem, I'm not sure what addCube function you refer to...

[SquirrelKnight]

Hi @bernard-giroux. I think I mistakenly referred to gmsh.model.occ.addBox as addCube! I've found other ways of working around the issue (manually creating layered boxes and assigning interpolated velocities to each node), but it is interesting that nodes generated for gmsh cubes don't seem to be fully compatible with ttcrpy.