Improve interpolation of fields at receiver positions.

[prisae]

Currently, the field at receiver positions is interpolated linearly. This can lead to quite large errors as we cover many orders of magnitude in EM; the problem is particularly pronounced with stretched grids. Here an example in 2D, but it applies equally for 3D:

This shows a dipole in a homogeneous fullspace, calculated with analytical solutions (hence no modelling errors); once it is calculated on a fine 5m x 5m grid (data shown on the left); and once on a stretched grid with stretching factor 1.1, having cell widths starting at 0.6 m up to 243 m. The middle figure shows the error if interpolated linearly, the right figure if interpolated with a cubic spline. Everything red means a relative error > 1%, hence not acceptable. Dark red is a relative error of 100% and more. (Notebook.)

There are various potential approaches for this, e.g.

• Implement a custom version in discretize, as is done in the case of the linear interpolation, which lives in interputils_cython.pyx;
• Write a little wrapper using SciPy along the lines of interp_cubic_spline_3d.py.
• @banesullivan thinks that it could be possible through PyVista utilizing https://vtk.org/Wiki/VTK/Image_Interpolators;
• @leouieda thinks that it should not be too difficult to expand Verde to do that, Support n-dimensional coordinates fatiando/verde#189.

[micmitch]

Pretty sure that I am fighting with similar problems at the moment with Point_e receivers for a grounded LineCurrent src. Src is x-oriented. When testing against a numerically integrated solution in fullspace the x component measurements of the E-field are accurate but I also get reasonably large y components measurements. Directly above the src the y component should be 0. The magnitude of these numeric uncertainties varies substantially depending on where I place the receivers (i.e. edges, nodes, cell centers)

[prisae]

Possibly. You have a discretize-grid, right? Can you try this function?

def get_receiver(grid, field, rec_loc, method='cubic'):
    """Return field at receiver locations.
    Get the fields at rec_loc; dimension of field should be one of:
    - x: [nCx, nNy, nNz]
    - y: [nNx, nCy, nNz]
    - z: [nNx, nNy, nCz]
    """
    points = tuple()
    for i, coord in enumerate(['x', 'y', 'z']):
        if field.shape[i] == getattr(grid, 'nN'+coord):
            pts = (getattr(grid, 'vectorN'+coord), )
        else:
            pts = (getattr(grid, 'vectorCC'+coord), )

        # Add to points
        points += pts

        # We need at least 3 points in each direction for cubic spline.
        # This should never be an issue for a realistic 3D model.
        if pts[0].size < 4:
            method = 'linear'

    # Interpolation.
    if method == "linear":
        ifn = interpolate.RegularGridInterpolator(
                points=points, values=field, method="linear",
                bounds_error=False, fill_value=0.0)

        return ifn(xi=rec_loc)

    else:

        # Replicate the same expansion of xi as used in
        # RegularGridInterpolator, so the input xi can be quite flexible.
        xi = interpolate.interpnd._ndim_coords_from_arrays(rec_loc, ndim=3)
        xi_shape = xi.shape
        xi = xi.reshape(-1, 3)

        # map_coordinates uses the indices of the input data (field in this
        # case) as coordinates. We have therefore to transform our desired
        # output coordinates to this artificial coordinate system too.
        params = {'kind': 'cubic',
                  'bounds_error': False,
                  'fill_value': 'extrapolate'}
        x = interpolate.interp1d(
                points[0], np.arange(len(points[0])), **params)(xi[:, 0])
        y = interpolate.interp1d(
                points[1], np.arange(len(points[1])), **params)(xi[:, 1])
        z = interpolate.interp1d(
                points[2], np.arange(len(points[2])), **params)(xi[:, 2])
        coords = np.vstack([x, y, z])

        # map_coordinates only works for real data; split it up if complex.
        if 'complex' in field.dtype.name:
            real = ndimage.map_coordinates(field.real, coords, order=3)
            imag = ndimage.map_coordinates(field.imag, coords, order=3)
            result = real + 1j*imag
        else:
            result = ndimage.map_coordinates(field, coords, order=3)

    return result.reshape(xi_shape[:-1])

• grid is a discretize-grid.
• field is the field in x-, y-, or z-direction.
• rec_loc are the receiver locations, as (x, y, z). They can be floats, arrays, or even matrices. They do not have to be the same shape, but it must be possible to expand them to the same shape. E.g.:
  • x.shape=y.shape=z.shape=(n,): e.g., (np.array([1000, 2000, 3000]), np.array([0, 0, 0]), np.array([0, 0, 0]))
  • x,shape=(n,) y.shape=z.shape=(1,): e.g. (np.array(1000, 2000, 3000]), 0, 0)
    or other combinations.

[prisae]

If it helps you could also describe me the exact survey parameters, and I could calculate the responses with the analytical solution in empymod for comparison.

[micmitch]

Thanks! I'm trying to put together a simple example now. Takes a few minutes to calculate the fields. I'll send you the notebook when I get it running!

[prisae]

Cool. This is the same result as above, but for the y-field at the receiver; x-directed point-dipole.

Once I have your notebook I'll make the same comparison.

[micmitch]

I forgot to mention that I am doing all of this on a octree mesh.... This is going to complicate things.

[jcapriot]

On the rectilinear type meshes, implementing the higher order interpolation shouldn’t be too much of an issue, just that it’ll take longer to construct those matrices.

However, @micmitch, right now it is pretty hard so do more than nearest, or linear on a TreeMesh in general. This is still something I have been thinking on for a long while, just haven’t come up with a good thought yet. I have been experimenting with triangulating them, and I’ve found a decent way to eliminate degenerate simplicies in the scipy triangulation, but not quite sure if it’s worth it for the time.

[micmitch]

Took a stab at rewriting the function for octree but didn't really get anywhere... Not sure that this approach as feasible with an octree mesh.

[micmitch]

That's fair. Thanks for letting me know before I spend too much more time trying to sort through this @jcapriot.

[prisae]

If it is just in modelling or at the end of an inversion to get the signal at receiver position then I think it is not time critical, so it wouldn't matter if the interpolation takes a second or so.

With an octree the approaches mentioned by @banesullivan with PyVista or by @leouieda with Verde might be the better approach.

[micmitch]

Certainly something that is worth looking into when we have the time!

[banesullivan]

I'm running low on time to experiment with this - but it's on my list!

I'm thinking this will be totally doable with VTK via PyVista, however the implementation will likely be as follows:

There will be an input dataset (the fine mesh) and a source dataset (the coarse mesh with data). Then the input dataset (triangular, octree, rectilinear, or whatever) will have a uniform mesh of equivalent/greater coverage made. It's on this new uniform mesh that we could run a spline interpolation from the source dataset. Then the input dataset (the triangular or octree) would probe/sample the uniform mesh for data values. This is a bit cumbersome but it should have compelling (and fast) results for 2- and 3-D problems on all mesh types (regular, tree, or tetrahedral).

What I will likely end up doing is creating a new filter in PyVista called interpolate_bspline which will have similar usage to interpolate to perform exactly this. Note that the interpolate filter has this same workflow, but uses a radial Gaussian kernel for the interpolation - would there be interest from this community to add linear, ellipsoidal-gaussian, and/or probabilistic Voronoi kernels for interpolations?

[prisae]

@lheagy , you said on Slack that in inversions this is a bit less of a concern, as you usually have designed the mesh so that all receivers are in the fine region. However, I think even then it could be beneficial. So if we improve the interpolation, then you would not have to put all receivers in the fine region, so you can end up with a grid with fewer cells and therefore faster inversions. I did some more testing and am convinced that for EM we should not use linear interpolation, not even with constant spacing.