Does change_grid function work with nonuniform grids?

[malamast]

Hi all,

I’m trying to use the change_grid function to remap restart files onto a higher-resolution grid:
https://github.com/boutproject/boutdata/blob/master/src/boutdata/restart.py#L1013

However, I’m finding that the interpolated solution differs significantly from the solution.

It seems that the interpolation is performed in a normalized index space, e.g.:
https://github.com/boutproject/boutdata/blob/master/src/boutdata/restart.py#L1176
which I think assumes a uniform grid in x within each region. In my case, the grid was generated with Hypnotoad and is non-uniform in the radial direction.

I’m planning to modify the function to account for non-uniform grids, but I wanted to ask:
What coordinate would be most appropriate for interpolation in this case? Perhaps the physical coordinates such as Rxy and Zxy?

Any guidance would be appreciated.
Cheers!

[bendudson]

Hey @malamast ! Yes I think that's right: It just interpolates in index space, which is ok in some cases but not in general. I think probably poloidal flux space (psixy) might be best. I'm not sure what you should do if the grids don't fully overlap. I guess you have to extrapolate.

[johnomotani]

This kind of issue was one of the things boutproject/xBOUT#231 was meant to help with. If you're doing 2D transport simulations that PR might be usable as-is, although it could still use some work from someone to test/document (and add toroidal interpolation for 3D simulations).

[malamast]

Thanks @johnomotani. I appreciate that. How do you usually call it?

I try to interpolate a 2D solution (either from dump or restart files) from Hermes-3 to a higher resolution grid and then save them as restart files.

I do something like that:

    ds_1 = xbout.open_boutdataset(datapath=file_1+"/BOUT.dmp.*.nc", geometry="toroidal",
                                  gridfilepath=gridFile_1, inputfilepath=file_1+"/BOUT.inp")

    ds_2 = ds_1.bout.interpolate_to_new_grid(
        variables=['Pe', 'Pd+'],
        new_gridfile=gridFile_2,
        method="cubic",
    )

    ds_2.to_restart(savepath=file_2, nxpe=1, nype=16)        

but it gives me this error:
ValueError: Dimensions {'theta'} do not exist. Expected one or more of ('x', 'y')
in line:
https://github.com/boutproject/xBOUT/blob/new-grid-interpolation/xbout/boutdataarray.py#L742

[johnomotani]

@malamast not sure why, sorry I made the PR for a student and have never used it myself! I was never a big fan of the dimension-renaming feature, and it looks like that's what's going wrong. Apparently at the line you reference, new_gridfile has dimensions {"x","y"}, rather than {"psi_poloidal","theta"}. Maybe there is some special handling for grid files now? Grid files contain a theta variable, that conflicts with renaming the 'y' dimension to 'theta', so maybe the renaming has been disabled for grid files.

Possibly the cleanest simple fix (if it works) would be to define

        newgrid_xcoord = new_gridfile.metadata["bout_xdim"]
        newgrid_ycoord = new_gridfile.metadata["bout_ydim"]
        newgrid_zcoord = new_gridfile.metadata["bout_zdim"]

and then use those instead of xcoord/ycoord/zcoord anywhere they apply to new_gridfile? Long-term, my preference would be to get rid of the dimension renaming, stick with {"t", "x", "y", "z"} everywhere and get rid of all the messing around with "bout_xdim", etc., but that's probably too big a job for it to be likely to happen, even if the current xbout developers thought it was a good idea, and it would probably also break many users' code.