Smarter grid discretization for regular grids.

[thurinj]

Following a recent email thread (@carltape, @ammcpherson, @rmodrak), we've identified a flaw in the current design of the regular grid definition. While it follows uniformity requirements in each direction/axes of the grid, the current definition based on a constant number of point per axis (npts_per_axis) does not take into account the dimensionality of each axes.

For example, in the orientation space (which is also relevant for the spherical domain of the point force grid):

• strike spans 360°, while dip and rake each span 180°

So ideally, to achieve a consistent density of points in the moment tensor space, strike should have ~2× more points than dip or rake. Contrarily to the uniform random grid, which ensures that the density of samples accross the space is uniform, the samples density is not preserved with a constant npts_per_axis.

@carltape proposed to add the ability to request a "base" number of points per axis, and the returned grid would account for a required number of samples appropriate for each direction in space:

So a user would ask for 10 "base" pts per axis and then get something like these numbers of axes: force phi 20 force theta 10 strike 20 dip 10 rake 10 v 10 w 35

which makes a lot of sense. We would have to ensure that the deviatoric and DC solution are always returned within the grid. A potential limitation might be faced with smaller sized grid, as low resolution might not allow to honor accurately the axes extent, but rather give a lose approximation. There might be backward compatibility issues/considerations in implementing this solution, as pointed out by @rmodrak .

A proposal solution would be to add a fully abstracted logic following @carltape 's proposition to work from a "base" level, and hide away all of the grid generation logic, or let the user define a fully custom grid extent as a dictionary: ex: npts_per_axis={'strike': 20, 'dip': 10, 'rake': 10, 'v': 10: 'w' 35} , which I also like as it allows having completely custom grid, but I recognize it might be super confusing for users.

What I propose would be something along these lines:

• Allow npts_per_axis to be passed as a dictionary to define per-axis resolution explicitly.
• If a single integer is passed, treat it as a base resolution and internally scale the number of points per axis using sensible defaults based on physical extent (e.g., 2× for strike, 3.5× for w) -- that would allow to include a level of abstraction that would be suitable for most users, and the other options would be for when you know what you're doing with the grid generation.
• Ensure that key physically meaningful configurations (e.g., deviatoric source with w=0, or DC with v=w=0) are always included in the returned grid.

And we will have to accept that for low values of the "base" resolution, the scaling might only be approximate... I think this would make sure most script are backward compatible (although the results from running inversion with and without this grid update would be different, by virtue of having different grids).

(I have left out the discussion related to the "intermediate" meshing/gridding strategies that are between the fully regular Cartesian grids and the uniform random grid.)