Generation of grid points in boutdata.resize #81
Conversation
boutdata/restart.py
Outdated
| # be between 0 and 1 Calculate the old coordinates | ||
| xCoordOld = np.arange(nx - mxg + shift)*dx | ||
| yCoordOld = np.arange(ny - myg + shift)*dy | ||
| zCoordOld = np.arange(nz + shift)*dz |
There was a problem hiding this comment.
@mredenti have you tested this? It doesn't look to me like it would work:
np.arange()as used here will always return integer values- both the 'old' and 'new' coordinates will have the wrong size - i.e. different number of points from the corresponding dimension in the old/new data array.
I think you need to specify both the start and stop arguments to np.arange()...
There was a problem hiding this comment.
@mredenti have you tested this? It doesn't look to me like it would work:
np.arange()as used here will always return integer values- both the 'old' and 'new' coordinates will have the wrong size - i.e. different number of points from the corresponding dimension in the old/new data array.
I think you need to specify both the
startandstoparguments tonp.arange()...
@johnomotani I had tested the implementation in my own post-processing code but somehow forgot to copy other some parenthesis. This has now been fixed in the new commit.
Also:
- The implementation will work fine with the nearest-neighbour interpolation algorithm but for the 'linear' method I think we need to enable the user to supply
kwargs : {'fill_value' : 'extrapolate'}to thescipy regulargrid interpolatorsince the first/last guard point of the new resolution will be to the left/right of the old resolution grid points when downsizing. If this additional parameter is not passed, scipy will fill the values at these grid points with NaN. (I am currently investigating this) - This should work with multiple restart files as well now or am I missing something? I should run some tests: perhaps I could simply check that the resized data is the same whether I start with multiple restart files or a single one. This should be the case for nearest neighbour and linear interpolation
There was a problem hiding this comment.
That is a good point. In boundary or guard cells we could just not worry about the NaNs as they should be replaced by communication or boundary conditions before being used. For the z-direction though there are no boundary or guard cells, so we would end up with NaNs on grid points. As far as I remember there is no way to tell scipy's interpolation routines that some dimension is periodic. It would be possible to fix the affected points by doing the interpolation 'by hand' using periodicity, but I think just extrapolating is good enough. TL;DR I think we should just pass kwargs : {'fill_value' : 'extrapolate'} all the time.
For multiple restart files, I the function works as documented for multiple files. The only issue is that newNx, newNy and newNz are a bit inconvenient - at the moment you have to pass values that depend on how many restart files you have. I'd prefer for newNx, newNy and newNz to accept the same values that would be in the BOUT.inp for the new simulation, which would require reading the original global grid sizes and using those to calculate the per-file new grid sizes from the global inputs; would also require some tweaks to the logic for when boundary points should or should not be included (see https://bout-dev.readthedocs.io/en/latest/user_docs/bout_options.html#grids).
There was a problem hiding this comment.
That is a good point. In boundary or guard cells we could just not worry about the NaNs as they should be replaced by communication or boundary conditions before being used.
In other words, it does not really matter what values end up on the boundary guard cells and ghost guard cells since BOUT++ will enforce the correct values before any derivative evaluation?
For the z-direction though there are no boundary or guard cells, so we would end up with NaNs on grid points. As far as I remember there is no way to tell
scipy's interpolation routines that some dimension is periodic. It would be possible to fix the affected points by doing the interpolation 'by hand' using periodicity, but I think just extrapolating is good enough. TL;DR I think we should just passkwargs : {'fill_value' : 'extrapolate'}all the time.
Yes, scipy's interpolation routine does not support interpolation with periodic boundaries which is a shame really. I guess we could pad the sides with a copy of the data on the opposite side/row but it becomes quite fiddly.
For multiple restart files, I the function works as documented for multiple files. The only issue is that
newNx,newNyandnewNzare a bit inconvenient - at the moment you have to pass values that depend on how many restart files you have. I'd prefer fornewNx,newNyandnewNzto accept the same values that would be in theBOUT.inpfor the new simulation, which would require reading the original global grid sizes and using those to calculate the per-file new grid sizes from the global inputs; would also require some tweaks to the logic for when boundary points should or should not be included (see https://bout-dev.readthedocs.io/en/latest/user_docs/bout_options.html#grids).
This would be a nice feature to add. So I need to be careful with the following
nx is the number of points in X including the boundaries,
ny and nz are the number of points in Y and Z not including the boundaries
and the fact that boutdata.restart.resize accepts a MXG and MYG parameter
There was a problem hiding this comment.
@mredenti yes to all of that 👍
There was a problem hiding this comment.
The commit b98f735 takes care of extrapolation if the grid points to be interpolated are outside the range. It turns out one had to pass the argument fill_value = None to the RegularGridInterpolator. The default fill_value is np.nan.
Accounts for variables that may be staggered to CELL_XLOW, CELL_YLOW, or CELL_ZLOW locations.
When resizing from a fine grid to a coarser grid, the location of the boundary guard cells will be to the left/right of the left/right guard cell of the finer grid. This is extrapolation which by default the RegularGridInterpolator method sets to `nan`. BOUT++ might struggle to interpret this at run time so I have changed this to `0.0`
… and instead extrapolate by setting fill_value=None
This PR is in regards to the issue #78 (comment)
This patch resolves the issue of inconsistency when generating grid points for interpolation. One problem with the current implementation is that values at guard cells would end up inside the physical domain with the nearest neighbour algorithm when upsizing, i.e. increasing the resolution ].