1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
|
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
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()
@staticmethod
def backend():
return 'netcdf4'
def open(self, mode='r+'):
from netCDF4 import Dataset
self.f = Dataset(self.filename, mode=mode, parallel=True, comm=comm)
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()
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()
|
