Problem with xBOUT-created restart files and restart.resize()

[johnomotani]

From the Slack:
Amy Krause 3:28 PM
Sadly that doesn't seem to work. I can resize the restart files from my simulation run. If I use to_restart() on the output files it produces the restart files, however I get an error when using restart.resize() on those
ValueError: cannot reshape array of size 5120 into shape (1,)

John Omotani 3:32 PM
that's odd. Do the restart files produced by xBOUT look OK (e.g. if you look at them with ncdump, or open them in Python with netCDF4.Dataset?

Amy Krause 3:37 PM
the files open ok in python

John Omotani 3:37 PM
not sure what's going on then. Sorry, I haven't used restart.resize() myself. Anyone else have any ideas?

Amy Krause 3:40 PM
Writing vort
Traceback (most recent call last):
File "src/netCDF4/_netCDF4.pyx", line 4916, in netCDF4._netCDF4.Variable.setitem
ValueError: cannot reshape array of size 5120 into shape (1,)

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
File "/work/d175/d175/akexcml/envs/xbout/lib/python3.8/site-packages/boututils/datafile.py", line 666, in write
var.assignValue(data)
File "src/netCDF4/_netCDF4.pyx", line 4963, in netCDF4._netCDF4.Variable.assignValue
File "src/netCDF4/_netCDF4.pyx", line 4918, in netCDF4._netCDF4.Variable.setitem
File "<array_function internals>", line 5, in broadcast_to
File "/work/d175/d175/akexcml/envs/xbout/lib/python3.8/site-packages/numpy/lib/stride_tricks.py", line 411, in broadcast_to
return _broadcast_to(array, shape, subok=subok, readonly=True)
File "/work/d175/d175/akexcml/envs/xbout/lib/python3.8/site-packages/numpy/lib/stride_tricks.py", line 348, in _broadcast_to
it = np.nditer(
ValueError: input operand has more dimensions than allowed by the axis remapping
3:41
That's the full error, it happens when it tries to write the variable (in our case vort) - I tried specifying one variable or just leaving the default (all time-evolving vars)

John Omotani 3:58 PM
Could you post the output of ncdump -h BOUT.restart.0.nc (for one of the files produced by xBOUT) so I can check the variable attributes please?
I'm slightly suspicious of the use multiprocessing (don't see anything obviously wrong, but it's often tricky...). Might be worth trying to pass the argument maxProc=1 to restart.resize() so that it only uses 1 process - might be slightly safer...
Otherwise, you could try to check the dimensions of new and newData before the write here

boutdata/boutdata/restart.py

Line 251 in 07a2473

new.write(var, newData)

e.g.
print("variable dims", newData.shape)
if len(new.list()) > 0:
var0 = new.list()[0]
print("new file dims", new.size(var0))
(hoping that var0 happens to be a 3d variable...)

restart.py
new.write(var, newData)
boutproject/boutdata | Added by GitHub

Amy Krause 4:13 PM
netcdf BOUT.restart.0 {
dimensions:
x = 9 ;
y = 1 ;
z = 128 ;
variables:
double psi(x, y, z) ;
psi:_FillValue = NaN ;
double U_x(x, y, z) ;
U_x:_FillValue = NaN ;
double U_y(x, y, z) ;
U_y:_FillValue = NaN ;
double U_z(x, y, z) ;
U_z:_FillValue = NaN ;
double divU(x, y, z) ;
divU:_FillValue = NaN ;
int ncalls ;
int ncalls_e ;
int ncalls_i ;
double vort(x, y, z) ;
vort:_FillValue = NaN ;
double dx(x, y) ;
dx:_FillValue = NaN ;
double dy(x, y) ;
dy:_FillValue = NaN ;
double dz ;
dz:_FillValue = NaN ;
double g11(x, y) ;
g11:_FillValue = NaN ;
double g22(x, y) ;
g22:_FillValue = NaN ;
double g33(x, y) ;
g33:_FillValue = NaN ;
double g12(x, y) ;
g12:_FillValue = NaN ;
double g13(x, y) ;
g13:_FillValue = NaN ;
double g23(x, y) ;
g23:_FillValue = NaN ;
double g_11(x, y) ;
g_11:_FillValue = NaN ;
double g_22(x, y) ;
g_22:_FillValue = NaN ;
double g_33(x, y) ;
g_33:_FillValue = NaN ;
double g_12(x, y) ;
g_12:_FillValue = NaN ;
double g_13(x, y) ;
g_13:_FillValue = NaN ;
double g_23(x, y) ;
g_23:_FillValue = NaN ;
double J(x, y) ;
J:_FillValue = NaN ;
int64 zperiod ;
int64 MZSUB ;
int64 MXG ;
int64 MYG ;
int64 MZG ;
int64 nx ;
int64 ny ;
int64 nz ;
int64 MZ ;
int64 NZPE ;
int64 ixseps1 ;
int64 ixseps2 ;
int64 jyseps1_1 ;
int64 jyseps2_1 ;
int64 jyseps1_2 ;
int64 jyseps2_2 ;
int64 ny_inner ;
double ZMAX ;
ZMAX:_FillValue = NaN ;
double ZMIN ;
ZMIN:_FillValue = NaN ;
double BOUT_VERSION ;
BOUT_VERSION:_FillValue = NaN ;
int64 MXSUB ;
int64 MYSUB ;
int64 NXPE ;
int64 NYPE ;
int64 hist_hi ;
double tt ;
tt:_FillValue = NaN ;
}

John Omotani 4:14 PM
Thanks @amy Krause - that ncdump looks OK to me. Don't know what's causing the error!

Amy Krause 4:15 PM
Thanks - yes the only thing that looks different is that there are a few more attributes in the original restart files. Which shouldn't matter.
4:19
like this:
double dx(x, y) ;
dx:cell_location = "CELL_CENTRE" ;
dx:direction_y = "Standard" ;
dx:direction_z = "Average" ;
4:20
actually, this is different too (does it matter?):
dimensions:
x = 9 ;
y = 1 ;
z = 128 ;
t = UNLIMITED ; // (0 currently)
4:21
vs xarray produced
dimensions:
x = 9 ;
y = 1 ;
z = 128 ;

John Omotani 4:25 PM
I don't think restart.resize() is using any of the attributes (cell_location, etc.) (pretty sure it was written before they existed!) so indeed that shouldn't be the problem.
I wouldn't have thought that the zero-size time dimension would be an issue, but maybe it is somehow. Do any of the variables in the original restart files have a t-dimension? I'd expect them not to, but maybe I'm wrong and that's confusing DataFile somehow...

Amy Krause 4:28 PM
no, resize pulls out only variables with 3 dimensions...
4:28
I'll see if I can dig in deeper
Today
New

Amy Krause 9:24 AM
I took out multiprocessing jobs and copied the rest of the code and it works fine (maxProc=1 didn't seem to help).

John Omotani 9:37 AM
@amy Krause thanks! glad you found a solution 

[mredenti]

@johnomotani I have run into the same problem when using the xBOUT to create the restart files and the BOUTDATA library to resize the restart files (increasing/decreasing resolution).

Could the problem have something to do with the following?

When creating the restart files with the boutdata.restart.create(...) method we have that

from boututils.datafile import DataFile 
df = DataFile('BOUT.restart.0.nc'), write=False)
print(df.dimensions('vort')) # ('t', 'x', 'y')

whereas if I use the XBOUT to_restart(...) method the dimensions of the variables are ('x', 'y', 'z')

from boututils.datafile import DataFile 
df = DataFile('BOUT.restart.0.nc'), write=False)
print(df.dimensions('vort')) # ('x', 'y', 'z')

The error message when to trying to resize with boutdata.restart.resize(...) contains the following

Traceback (most recent call last):
  File "/mnt/lustre/a2fs-work3/work/ta070/ta070/mich/miniconda3/envs/post-processing/lib/python3.9/site-packages/boututils/datafile.py", line 674, in write
    var.assignValue(data)
  File "src/netCDF4/_netCDF4.pyx", line 4963, in netCDF4._netCDF4.Variable.assignValue
  File "src/netCDF4/_netCDF4.pyx", line 4918, in netCDF4._netCDF4.Variable.__setitem__
  File "<__array_function__ internals>", line 180, in broadcast_to
  File "/mnt/lustre/a2fs-work3/work/ta070/ta070/mich/miniconda3/envs/post-processing/lib/python3.9/site-packages/numpy/lib/stride_tricks.py", line 412, in broadcast_to
    return _broadcast_to(array, shape, subok=subok, readonly=True)
  File "/mnt/lustre/a2fs-work3/work/ta070/ta070/mich/miniconda3/envs/post-processing/lib/python3.9/site-packages/numpy/lib/stride_tricks.py", line 348, in _broadcast_to
    it = np.nditer(
ValueError: input operand has more dimensions than allowed by the axis remapping

which makes me think the 'z' dimension is the "extra" dimension that boutdata.restart.resize is not expecting.

In the meantime, is there a function in the XBOUT library that increases/decreases the resolution from a BOUT.restart file? (On a grid)

[johnomotani]

I don't know - the fact that the routine apparently worked on xBOUT-produced restart files if the multiprocessing stuff was removed makes me think maybe the dimension names aren't the problem. No idea what the issue is though.

I started making some interpolation routines for xBOUT here boutproject/xBOUT#231. They're intended to be a lot more general than restart.resize(), and aren't finished yet. If anyone wants to continue, they're welcome! Otherwise that PR probably won't get finished unless/until I need it for something...

[mredenti]

the fact that the routine apparently worked on xBOUT-produced restart files if the multiprocessing stuff was removed makes me think maybe the dimension names aren't the problem.

How do you remove the multiprocessing stuff?

Also, could I use https://github.com/boutproject/xBOUT/blob/fe0e8eecb2966a1e7ef380a2c168891ef17f2000/xbout/boutdataset.py#L181-L188

to increasing the resolution of a (x=Nx, y=1, z=Nz) cartesian grid?

[johnomotani]

How do you remove the multiprocessing stuff?

I guess call resize3DField() directly instead of using the multiprocessing.Pool.

Also, could I use https://github.com/boutproject/xBOUT/blob/fe0e8eecb2966a1e7ef380a2c168891ef17f2000/xbout/boutdataset.py#L181-L188
to increasing the resolution of a (x=Nx, y=1, z=Nz) cartesian grid?

No, that interpolates in the parallel direction (i.e. y-dimension).