Volume averaging interpolation

[prisae]

As a follow up of questions I asked on Slack (Simpeg#discretize and Swung#python). EM fields are arguably best interpolated by a cubic spline (#175). Up- or downscaling model parameters is, however, often best done with a method which keeps the total value constant (X*Y*Z*value = constant), something which is not given with linear nor cubic spline interpolation. Although there are more sophisticated solutions (based on effective medium theory) the volume averaging technique is a simple method to ensure the total value of the property remains unchanged.

I implemented the technique in emg3d, following Plessix et al., 2007, Geophysics. I post it here if you ever want to implement it in discretize. This is a numba-version, but it could be re-written for a Cython (unless discretize makes numba a dependency):

import numba as nb
import numpy as np


@nb.njit(nogil=True, fastmath=True, cache=True)
def volume_average(edges_x, edges_y, edges_z, values,
                   new_edges_x, new_edges_y, new_edges_z, new_values):
    """Interpolation using the volume averaging technique.

    The result is added to new_values.


    Parameters
    ----------
    edges_[x, y, z] : ndarray
        The edges in x-, y-, and z-directions for the original grid.

    values : ndarray
        Values corresponding to ``grid``.

    new_edges_[x, y, z] : ndarray
        The edges in x-, y-, and z-directions for the new grid.

    new_values : ndarray
        Array where values corresponding to ``new_grid`` will be added.

    """

    # Get cell indices.
    # First and last edges ignored => first and last cells extend to +/- infty.
    ix_l = np.searchsorted(edges_x[1:-1], new_edges_x, 'left')
    ix_r = np.searchsorted(edges_x[1:-1], new_edges_x, 'right')
    iy_l = np.searchsorted(edges_y[1:-1], new_edges_y, 'left')
    iy_r = np.searchsorted(edges_y[1:-1], new_edges_y, 'right')
    iz_l = np.searchsorted(edges_z[1:-1], new_edges_z, 'left')
    iz_r = np.searchsorted(edges_z[1:-1], new_edges_z, 'right')

    # Get number of cells.
    ncx = len(new_edges_x)-1
    ncy = len(new_edges_y)-1
    ncz = len(new_edges_z)-1

    # Working arrays for edges.
    x_edges = np.empty(len(edges_x)+2)
    y_edges = np.empty(len(edges_y)+2)
    z_edges = np.empty(len(edges_z)+2)

    # Loop over new_grid cells.
    for iz in range(ncz):
        hz = np.diff(new_edges_z[iz:iz+2])[0]  # To calc. current cell volume.

        for iy in range(ncy):
            hyz = hz*np.diff(new_edges_y[iy:iy+2])[0]  # " "

            for ix in range(ncx):
                hxyz = hyz*np.diff(new_edges_x[ix:ix+2])[0]  # " "

                # Get start edge and number of cells of original grid involved.
                s_cx = ix_r[ix]
                n_cx = ix_l[ix+1] - s_cx

                s_cy = iy_r[iy]
                n_cy = iy_l[iy+1] - s_cy

                s_cz = iz_r[iz]
                n_cz = iz_l[iz+1] - s_cz

                # Get the involved original grid edges for this cell.
                x_edges[0] = new_edges_x[ix]
                for i in range(n_cx):
                    x_edges[i+1] = edges_x[s_cx+i+1]
                x_edges[n_cx+1] = new_edges_x[ix+1]

                y_edges[0] = new_edges_y[iy]
                for j in range(n_cy):
                    y_edges[j+1] = edges_y[s_cy+j+1]
                y_edges[n_cy+1] = new_edges_y[iy+1]

                z_edges[0] = new_edges_z[iz]
                for k in range(n_cz):
                    z_edges[k+1] = edges_z[s_cz+k+1]
                z_edges[n_cz+1] = new_edges_z[iz+1]

                # Loop over each (partial) cell of the original grid which
                # contributes to the current cell of the new grid and add its
                # (partial) value.
                for k in range(n_cz+1):
                    dz = np.diff(z_edges[k:k+2])[0]
                    k += s_cz

                    for j in range(n_cy+1):
                        dyz = dz*np.diff(y_edges[j:j+2])[0]
                        j += s_cy

                        for i in range(n_cx+1):
                            dxyz = dyz*np.diff(x_edges[i:i+2])[0]
                            i += s_cx

                            # Add this cell's contribution.
                            new_values[ix, iy, iz] += values[i, j, k]*dxyz

                # Normalize by new_grid-cell volume.
                new_values[ix, iy, iz] /= hxyz

[prisae]

Also, this is strictly for 3D, but for discretize it probably would have to be generalized to also work for 1D and 2D.

[prisae]

As with #175 , there might also be a way within PyVista to do the same, @banesullivan ?

[prisae]

Numba can, AFAIK, not access structure attributes, that is why the above function call is more complicated then necessary. But with a small wrapper this can be solved:

def grid2grid(grid, values, new_grid):
    """Interpolate ``values`` located on ``grid`` to ``new_grid``.

    Points on ``new_grid`` which are outside of ``grid`` are filled by the
    nearest value.


    Parameters
    ----------
    grid, new_grid : TensorMesh
        Input and output model grids; ``TensorMesh``-instances.

    values : ndarray
        Values corresponding to ``grid``.

    Returns
    -------
    new_values : ndarray
        Values corresponding to ``new_grid``.

    """

    points = (grid.vectorNx, grid.vectorNy, grid.vectorNz)
    new_points = (new_grid.vectorNx, new_grid.vectorNy, new_grid.vectorNz)
    
    # Pre-allocate output array.
    new_values = np.zeros(new_grid.vnC, dtype=values.dtype)

    # Calculate new values.
    volume_average(*points, values, *new_points, new_values)

    return new_values

And it can be used then like

new_values = grid2grid(grid, values, new_grid)

[prisae]

Note: This volume-averaging is massively slower than simple linear interpolation. (There might be smarter ways to implement it though...)

So in the current state it should not be used within an inversion or something, just to build models.

[lheagy]

Thanks @prisae: this is a good point! We have a few of the pieces for this - within discretize, I would see this being implemented via the construction of a sparse-matrix but same idea

[prisae]

Yes, that makes sense, I would love to see a sparse-matrix volume averaging implementation!

[jcapriot]

I bet we can get a nice quick cython implementations

[prisae]

Unless you change the algorithm per se I'd be surprised if a cython implementation is significantly faster than numba.

[prisae]

It would still be ideal though for discretize, because discretize has already the Cython dependency.

[banesullivan]

As with #175 , there might also be a way within PyVista to do the same, @banesullivan ?

Sorry, I just saw this, @prisae - can you explain this algorithm? There's likely a way to do it in PyVista/VTK - but I want to make sure I understand this correctly.

At a glance, this looks pretty similar to how PyVista's .interpolate() filter behaves - it might be worth trying that out - it needs a bit of improvement but it should work well for 3D Tensor Meshes

[prisae]

It simply means that the volume average of a property stays constant. So you don't just interpolate, say, from old cell centres to new cell centres. But you multiply the values with their dx*dy*dz to get the cell averaged value, and ensure the sum of those stays constant (for the same extent). Could very well be that this is how PyVista's interpolate works.

You could make a very small example in 2D:

    _ _ _ _ _
1. |_ _ _ _|_|
x: 0 1 2 3 4 5
    _ _ _ _ _
2. |_ _|_ _ _|

Say, in the first case, the left cell has value 100 (at cell centre x=2), and the right 10 (at cell centre x=4.5), so the total value is 4*1*100+1*1*10=410, so the interpolated values in the second case must also fulfil 2*1*c1+3*1*c2=410, where c1=100 and c2=70 therefore. Linear interpolation of old to new cell centres (new centres are x={1.5, 3.5}) would yield c2=46, and for c1 it would depend if you do extrapolation or what kind of rules you have for extrapolation for new values outside of your old values.

[jcapriot]

So, from what I'm seeing in the code essentially the operation is in three steps:

1. Multiply each value by the current mesh cells' volume
2. Distribute those values to every cell in the new mesh that overlaps each old mesh cell
   • for cells that partly overlap, distribute only the overlapping part?
3. Divide by the volume of the new mesh cells'

Correct me if i'm wrong

[prisae]

Correct @jcapriot, that is exactly what it does. And yes, it only contributes the overlapping part (your point with the '?').

[prisae]

Volume averaging was implemented in #212