4 Support functions for reading DISPATCH2 data. 8 from __future__
import print_function
10 if sys.version_info[0]!=2:
11 from importlib
import reload
14 from abc
import abstractproperty, ABCMeta
20 import stagger_utils
as su
23 """An abstract class to describe tasks.""" 24 __metaclass__ = ABCMeta
29 """An iterable that yields the data associated with the object.""" 33 """A derived class for patches.""" 34 __metaclass__ = ABCMeta
35 def __init__(self, filename, suffix='.txt', read_derivs=False, overlap=0.01, verbose=False, fd=None):
36 self.filename = filename
42 fp = open(filename+suffix,
'r') 43 legacy_reader.single_patch (self, fp, verbose=False, fd=fd)
44 self.active_corners[1] = self.active_corners[1] + self.ds*overlap
45 self.corners[1] = self.corners[1] + self.ds*overlap
46 if self.kind ==
'poisson':
50 if read_derivs: self.nvar = 2 * self.nvar
52 self.offset=(self.id-1)*np.product(self.nwrite)*self.nvar*4
55 self.varidx = self.variable_indices(read_derivs)
58 for k
in self.varidx.keys():
59 self.varnames.append(k)
60 for k
in self.varidx.keys():
67 The actual data for a patch is accessed through a numpy 'memory map', 68 which uses on-demand access rather than reading the whole file in. 69 Note that a file handle is opened when this is invoked for the first time 70 and, if there are a very large number of patches, python will raise 71 an exception when it reaches its file handler limit. 77 def _read_legacy_data(self, filename, dims, nvar, nbytes):
78 """Return a memory map pointing to MHD data.""" 79 myshape = dims + (nvar,)
85 return np.memmap(filename, dtype=thistype, offset=self.
offset, mode=
'r', order='F', shape=myshape)
87 def _determine_corners(self, active_only=False):
89 For a given patch, calculate the corners of the cuboid. 91 `corners` are defined as the absolute lowest coords. of a patch (i.e, 92 the first or last interface). `active_corners` are instead defined to 93 only include active cells. 98 for o,n,d,ng
in zip(self.origin, self.n, self.dx, self.nghost):
99 c_lo.append(o - ng * d)
100 c_hi.append(o + (n + ng) * d)
102 for o,n,d
in zip(self.origin, self.n, self.dx):
104 c_hi.append(o + n * d)
106 c_lo = np.array(c_lo)
107 c_hi = np.array(c_hi)
110 def stripv (self, v):
111 l=np.array(self.nghost)
115 return v[l[0]:u[0],l[1]:u[1],l[2]:u[2]]
117 if self.
kind[0:9]==
'rt_solver':
118 self.
t = self.
t.transpose(tr)
119 self.
Q = self.
Q.transpose(tr)
120 elif self.
kind[0:7]==
'poisson':
121 self.
d = self.
d.transpose(tr)
122 self.
phi = self.
phi.transpose(tr)
124 '''Remove guard cells from cached arrays''' 125 if not 'd' in dir(self):
127 self.
d = self.
d.transpose(tr)
128 self.
e = self.
e.transpose(tr)
129 self.
s = self.
s.transpose(tr)
130 self.
sd = self.
sd.transpose(tr)
132 self.
pg = self.
pg.transpose(tr)
133 self.
u1 = self.
u1.transpose(tr)
134 self.
u2 = self.
u2.transpose(tr)
135 self.
u3 = self.
u3.transpose(tr)
136 self.
p1 = self.
p1.transpose(tr)
137 self.
p2 = self.
p2.transpose(tr)
138 self.
p3 = self.
p3.transpose(tr)
140 self.
phi= self.
phi.transpose(tr)
142 self.
b1 = self.
b1.transpose(tr)
143 self.
b2 = self.
b2.transpose(tr)
144 self.
b3 = self.
b3.transpose(tr)
146 t1.append(self.
x); t1.append(self.
y); t1.append(self.
z)
147 self.
x=t1[tr[0]]; self.
y=t1[tr[1]]; self.
z=t1[tr[2]]
148 ac=self.active_corners
149 self.active_corners[0]=[ac[0][tr[0]],ac[0][tr[1]],ac[0][tr[2]]]
150 self.active_corners[1]=[ac[1][tr[0]],ac[1][tr[1]],ac[1][tr[2]]]
152 l=self.nghost[0]; u=l+self.n[0]; self.
x = self.
x[l:u]; self.
xs = self.
xs[l:u]
153 l=self.nghost[1]; u=l+self.n[1]; self.
y = self.
y[l:u]; self.
ys = self.
ys[l:u]
154 l=self.nghost[2]; u=l+self.n[2]; self.
z = self.
z[l:u]; self.
zs = self.
zs[l:u]
157 if self.
kind[0:9]==
'rt_solver':
160 elif self.
kind[0:7]==
'poisson':
164 '''Remove guard cells from cached arrays''' 165 if not 'd' in dir(self):
185 l=self.nghost; u=l+self.n
188 self.
x = self.
x[l[0]:u[0]]; self.
xs = self.
xs[l[0]:u[0]]
189 self.
y = self.
y[l[1]:u[1]]; self.
ys = self.
ys[l[1]:u[1]]
190 self.
z = self.
z[l[2]:u[2]]; self.
zs = self.
zs[l[2]:u[2]]
195 if self.
kind[0:9]==
'rt_solver':
199 elif self.
kind[0:7]==
'poisson':
216 if self.
kind[0:7]==
'stagger':
217 logd=np.log10(self.
d)
219 self.
u1=self.
p1/np.exp(su.xdn(ln10*logd))
220 self.
u2=self.
p2/np.exp(su.ydn(ln10*logd))
221 self.
u3=self.
p3/np.exp(su.zdn(ln10*logd))
222 if self.
kind[0:9]==
'stagger2_':
224 self.
e =np.exp((self.
s/self.
d+ln10*logd)*g1)/g1
225 self.
pg=self.
d*self.
e*g1
226 elif self.
kind[0:9]==
'stagger2e':
228 self.
pg=self.
d*self.
e*g1
229 self.
s=np.log(self.
pg/self.
d**self.gamma)/g1
230 elif self.
kind[0:6]==
'ramses':
232 self.
u1=self.
p1/self.
d 233 self.
u2=self.
p2/self.
d 234 self.
u3=self.
p3/self.
d 235 self.
e=self.
e-0.5*(self.
u1**2+self.
u2**2+self.
u3**2)
236 self.
pg=self.
d*self.
e*g1
237 self.
s=np.log(self.
pg/self.
d**self.gamma)/g1
238 self.
sd =self.
s/self.
d 243 '''Make arrays available in a .var dict''' 244 if not 'var' in dir(self):
246 if self.
kind==
'rt_solver':
247 self.
var={
'Q':self.
Q,\
251 self.
var={
'd':self.
d,\
261 'logd':np.log10(self.
d),\
262 'loge':np.log10(self.
e)}
264 self.
var[
'b1']=self.
b1 265 self.
var[
'b2']=self.
b2 266 self.
var[
'b3']=self.
b3 271 ''' array of corner radiii ''' 272 c=self.active_corners
274 for x
in (c[0,0],c[1,0]):
275 for y
in (c[0,1],c[1,1]):
276 for z
in (c[0,2],c[1,2]):
277 pt.append(x**2+y**2+z**2)
278 return np.sqrt(np.array(pt))
282 A simple function to provide the indices of the variable stored in the binary 285 assert isinstance(self.
kind, str),
"expecting a string for `self.kind`" 286 solver = self.
kind.lower()
287 if solver[0:7]==
'stagger':
289 'p1': 2,
'p2': 3,
'p3': 4}
290 if solver[7:9] ==
'2e':
309 'dp1': 10,
'dp2': 11,
'dp3': 12,
310 'db1': 13,
'db2': 14,
'db3': 15,}
311 indices.update(indicesd)
312 elif solver[0:6] ==
'ramses':
315 'p1': 1,
'p2': 2,
'p3': 3,
316 'b1': 5,
'b2': 6,
'b3': 7,
317 'v1': 8,
'v2': 9,
'v3':10}
321 'dp1': 10,
'dp2': 11,
'dp3': 12,
322 'db1': 13,
'db2': 14,
'db3': 15,}
323 indices.update(indicesd)
324 elif solver[0:4] ==
'zeus':
327 'p1': 3,
'p2': 4,
'p3': 5,
328 'b1': 6,
'b2': 7,
'b3': 8,
330 elif solver.lower() ==
'rt_solver':
334 elif solver.lower() ==
'poisson':
335 indices = {
'd': 0,
'phi': 1}
337 indicesd = {
'dd': 0,
'dphi': 1}
338 indices.updates(indicesd)
340 raise Exception(
"Unrecognised solver type! => "+solver_type.lower())
343 def Search(fd, ip, verbose=False):
344 line=fd.readline().split()
345 ioformat=int(line[0])
350 line=fd.readline().split()
359 nskip=np.int(line[1])
360 for i
in range(nskip):
367 def Patches(file, overlap=0.01, verbose=False):
368 with open(file+
'.txt',
'r') as fd: 370 line=fd.readline().split() 371 ioformat=int(line[0]) 375 id,nskip=fd.readline().split()
378 patch=
Patch(file,fd=fd,overlap=overlap,verbose=verbose)
379 patches.append(patch)
383 def __init__(self, iout=0, run='', data='../data', overlap=0.01, verbose=0):
386 dir=data+
'/{:05d}/'.format(iout)
388 dir=data+
'/'+run+
'/{:05d}/'.format(iout)
390 self.
rankfiles=[dir+f.replace(
'.txt',
'')
for f
in os.listdir(dir)
if f.endswith(
'.rank.txt')]
397 patches=Patches(file,overlap=overlap,verbose=verbose)
398 for patch
in patches:
402 file=self.
directory+
'{:05d}'.format(p.id)
403 self.
files.append(file)
409 self.
files=[dir+f.replace(suffix,
'')
for f
in os.listdir(dir)
if f.endswith(suffix)]
410 if len(self.
files) == 0:
413 self.
files=[dir+f.replace(suffix,
'')
for f
in os.listdir(dir)
if f.endswith(suffix)]
414 for file
in self.
files:
415 self.
patches.append (
Patch(file,suffix=suffix,overlap=overlap))
417 print(
'{} snapshot.patches ({:.2f} s)'.format(np.size(self.
patches),time()-start))
def _read_legacy_data(self, filename, dims, nvar, nbytes)
def variable_indices(self, read_derivs=False)