Interpolation methods

[johnomotani]

It would be nice to have some generic interpolation method to convert a Dataset to a higher (or lower) resolution grid. We could then use to_restart() to create restart files for a simulation to continue at a different resolution. This would be a more general version of BoutDataset.interpolate_parallel().

Things to remember/consider:

• update dx, dy and dz for the new resolution
• update the metadata (nx, nx, nz, topology indices, etc.) for the new resolution
• parallel and perpendicular (radial/toroidal) interpolation will need to be done separately, probably there should be an argument to the method to specify which should be done first
  • is it better to do the radial/toroidal interpolation as a 2d interpolation or as two 1d interpolations?
• would be nice to have a switch to choose whether to do the radial interpolation in standard or in field-aligned coordinates. Radial interpolation in standard coordinates is probably more accurate (doesn't have sheared radial/toroidal grids), but interpolation in field-aligned coordinates may produce structures that are more field-aligned in the result and reduce initial transients in restarted simulations.
• should interpolation in the toroidal direction be polynomial or fft-based? Does fft interpolation work (efficiently?) if the new nz is not a multiple of or a divisor of the old nz?
• cf Dataset interpolation efficiency #164
• an obvious option would be to call interpolate_parallel within this new interpolation method. However, it might be more efficient to re-implement the functionality (or to convert interpolate_parallel to a special case of this method) - in particular if we are operating region-by-region it would probably be best to split into regions once, do all the necessary interpolations, and then re-combine once, rather than splitting/recombining for the interpolate_parallel part and then splitting/recombining again for the perpendicular interpolation.
• for simulations using a grid file, it would be useful to be able to interpolate onto a grid from a new grid file
  • do we require (or assume) that the grids use the same equilbrium so psi(R,Z) is the same and we can interpolate in psi?
  • what do we do if the boundaries (radial and/or poloidal) are not in the same place? If they have moved inside the old grid we might be OK, but otherwise we'd need to extrapolate - maybe a choice of extrapolation methods (e.g. nearest-neighbour to start with)? An alternative might be that we just interpolate in some normalised x- and y-coordinates so the old grid gets stretched onto the new one?

[bendudson]

Agreed, this would be really nice functionality. Very useful for many different projects.

[johnomotani]

Another comment:
x-direction interpolation should be done region-by-region (so similar code to interpolate_parallel(), but without the shifting to/from field-aligned).
z-direction interpolation doesn't need to split into regions since there are no region boundaries in the z-direction, so should be simpler (although there's a question of whether we use FFT-based interpolation or polynomial interpolation).