import os
import numpy as np
from mpi4py import MPI
from .file_base import FileBase
# https://github.com/Unidata/netcdf4-python/blob/master/examples/mpi_example.py
# Note. Not using groups because Visit does not understand it
__all__ = ('NCFile',)
comm = MPI.COMM_WORLD
[docs]
class NCFile(FileBase):
"""Class for writing data to NetCDF4 format
Parameters
----------
ncname : str
Name of netcdf file to be created
domain : Sequence, optional
An optional spatial mesh or domain to go with the data.
Sequence of either
- 2-tuples, where each 2-tuple contains the (origin, length)
of each dimension, e.g., (0, 2*pi).
- Arrays of coordinates, e.g., np.linspace(0, 2*pi, N). One
array per dimension.
mode : str
``r``, ``w`` or ``a`` for read, write or append. Default is ``a``.
clobber : bool, optional
If True (default), opening a file with mode='w' will clobber an
existing file with the same name. If False, an exception will be
raised if a file with the same name already exists.
kw : dict, optional
Optional additional keyword arguments used when creating the file
used to store data.
Note
----
Each class instance creates one unique NetCDF4-file, with one step-counter.
It is possible to store multiple fields in each file, but all snapshots of
the fields must be taken at the same time. If you want one field stored
every 10th timestep and another every 20th timestep, then use two different
class instances with two different filenames ``ncname``.
"""
def __init__(self, ncname, domain=None, mode='a', clobber=True, **kw):
FileBase.__init__(self, ncname, domain=domain)
from netCDF4 import Dataset
# netCDF4 does not seem to handle 'a' if the file does not already exist
if mode == 'a' and not os.path.exists(ncname):
mode = 'w'
self.f = Dataset(ncname, mode=mode, clobber=clobber, parallel=True,
comm=comm, **kw)
self.dims = None
if 'time' not in self.f.variables:
self.f.createDimension('time', None)
self.f.createVariable('time', float, ('time'))
self.close()
def _check_domain(self, group, field):
N = field.global_shape[field.rank:]
if self.domain is None:
self.domain = []
for i in range(field.dimensions):
self.domain.append(np.linspace(0, 2*np.pi, N[i]))
assert len(self.domain) == field.dimensions
if len(self.domain[0]) == 2:
d = self.domain
self.domain = []
for i in range(field.dimensions):
self.domain.append(np.linspace(d[i][0], d[i][1], N[i]))
self.dims = ['time']
for i in range(field.rank):
ind = 'ijk'[i]
self.dims.append(ind)
if not ind in self.f.variables:
self.f.createDimension(ind, field.dimensions)
n = self.f.createVariable(ind, float, (ind))
n[:] = np.arange(field.dimensions)
for i in range(field.dimensions):
xyz = 'xyzrst'[i]
self.dims.append(xyz)
if not xyz in self.f.variables:
self.f.createDimension(xyz, N[i])
nc_xyz = self.f.createVariable(xyz, float, (xyz))
nc_xyz[:] = self.domain[i]
self.f.sync()
[docs]
@staticmethod
def backend():
return 'netcdf4'
[docs]
def open(self, mode='r+'):
from netCDF4 import Dataset
self.f = Dataset(self.filename, mode=mode, parallel=True, comm=comm)
[docs]
def write(self, step, fields, **kw):
"""Write snapshot ``step`` of ``fields`` to NetCDF4 file
Parameters
----------
step : int
Index of snapshot.
fields : dict
The fields to be dumped to file. (key, value) pairs are group name
and either arrays or 2-tuples, respectively. The arrays are complete
arrays to be stored, whereas 2-tuples are arrays with associated
*global* slices.
as_scalar : boolean, optional
Whether to store rank > 0 arrays as scalars. Default is False.
Example
-------
>>> from mpi4py import MPI
>>> from mpi4py_fft import PFFT, NCFile, newDistArray
>>> comm = MPI.COMM_WORLD
>>> T = PFFT(comm, (15, 16, 17))
>>> u = newDistArray(T, forward_output=False, val=1)
>>> v = newDistArray(T, forward_output=False, val=2)
>>> f = NCFile('ncfilename.nc', mode='w')
>>> f.write(0, {'u': [u, (u, [slice(None), 4, slice(None)])],
... 'v': [v, (v, [slice(None), 5, 5])]})
>>> f.write(1, {'u': [u, (u, [slice(None), 4, slice(None)])],
... 'v': [v, (v, [slice(None), 5, 5])]})
This stores the following datasets to the file ``ncfilename.nc``.
Using in a terminal 'ncdump -h ncfilename.nc', one gets::
netcdf ncfilename {
dimensions:
time = UNLIMITED ; // (2 currently)
x = 15 ;
y = 16 ;
z = 17 ;
variables:
double time(time) ;
double x(x) ;
double y(y) ;
double z(z) ;
double u(time, x, y, z) ;
double u_slice_4_slice(time, x, z) ;
double v(time, x, y, z) ;
double v_slice_5_5(time, x) ;
}
"""
self.open()
nc_t = self.f.variables.get('time')
nc_t.set_collective(True)
it = nc_t.size
if step in nc_t.__array__(): # If already stored at this step previously
it = np.argwhere(nc_t.__array__() == step)[0][0]
else:
nc_t[it] = step
FileBase.write(self, it, fields, **kw)
self.close()
[docs]
def read(self, u, name, **kw):
step = kw.get('step', 0)
self.open()
s = u.local_slice()
s = (step,) + s
u[:] = self.f[name][s]
self.close()
def _write_slice_step(self, name, step, slices, field, **kw):
assert name not in self.dims # Crashes if user tries to name fields x, y, z, .
rank = field.rank
slices = list((slice(None),)*rank + tuple(slices))
slname = self._get_slice_name(slices[rank:])
s = field.local_slice()
slices, inside = self._get_local_slices(slices, s)
sp = np.nonzero([isinstance(x, slice) for x in slices])[0]
sf = np.take(s, sp)
sdims = ['time'] + list(np.take(self.dims, np.array(sp)+1))
fname = "_".join((name, slname))
if fname not in self.f.variables:
h = self.f.createVariable(fname, field.dtype, sdims)
else:
h = self.f.variables[fname]
h.set_collective(True)
h[step] = 0 # collectively create dataset
h.set_collective(False)
sf = tuple([step] + list(sf))
sl = tuple(slices)
if inside:
h[sf] = field[sl]
h.set_collective(True)
self.f.sync()
def _write_group(self, name, u, step, **kw):
assert name not in self.dims # Crashes if user tries to name fields x, y, z, .
s = u.local_slice()
if name not in self.f.variables:
h = self.f.createVariable(name, u.dtype, self.dims)
else:
h = self.f.variables[name]
h.set_collective(True)
s = (step,) + s
h[s] = u
self.f.sync()