Merge branch 'master' into reconnection

This commit is contained in:
Grzegorz Kowal 2021-10-06 16:28:08 -03:00
commit d0e97eaba4
6 changed files with 278 additions and 265 deletions

View File

@ -5,7 +5,7 @@ with open("README.md", "r", encoding="utf-8") as fh:
setuptools.setup(
name="amunpy",
version="0.9.2",
version="0.9.4",
author="Grzegorz Kowal",
author_email="grzegorz@amuncode.org",
description="Python Interface for the AMUN code's snapshots",

View File

@ -21,6 +21,6 @@ __all__ = [ 'AmunXML', 'AmunH5', 'WriteVTK', \
__author__ = "Grzegorz Kowal"
__copyright__ = "Copyright 2018-2021, Grzegorz Kowal <grzegorz@amuncode.org>"
__version__ = "0.6.2"
__version__ = "0.9.4"
__maintainer__ = "Grzegorz Kowal"
__email__ = "grzegorz@amuncode.org"

View File

@ -186,6 +186,9 @@ class Amun:
if all(v in self.variables for v in ['velx','vely','velz']):
self.variables['velocity'] = 'vvec'
self.variables['velocity magnitude'] = 'velo'
self.variables['x-velocity'] = 'velx'
self.variables['y-velocity'] = 'vely'
self.variables['z-velocity'] = 'velz'
self.variables['divergence of velocity'] = 'divv'
self.variables['x-vorticity'] = 'vorx'
self.variables['y-vorticity'] = 'vory'
@ -195,6 +198,9 @@ class Amun:
if all(v in self.variables for v in ['magx','magy','magz']):
self.variables['magnetic field'] = 'bvec'
self.variables['magnetic field magnitude'] = 'magn'
self.variables['x-magnetic field'] = 'magx'
self.variables['y-magnetic field'] = 'magy'
self.variables['z-magnetic field'] = 'magz'
self.variables['divergence of magnetic field'] = 'divb'
self.variables['x-current density'] = 'curx'
self.variables['y-current density'] = 'cury'
@ -240,18 +246,18 @@ class Amun:
"""
Get dataset array of name dataset_name from the file n.
"""
import numpy as np
import numpy
dataset = self.variables[dataset_name]
if dataset == 'mlev':
dset = np.zeros(self.chunks[chunk_number]['dims'])
dset = numpy.zeros(self.chunks[chunk_number]['dims'])
for p in range(self.chunks[chunk_number]['dblocks']):
dset[...,p] = self.chunks[chunk_number]['levels'][p]
elif dataset == 'logd':
dset = np.log10(self.__read_binary_data__('dens', chunk_number))
dset = numpy.log10(self.__read_binary_data__('dens', chunk_number))
elif dataset == 'logp':
dset = np.log10(self.__read_binary_data__('pres', chunk_number))
dset = numpy.log10(self.__read_binary_data__('pres', chunk_number))
elif dataset == 'velo':
tmp = self.__read_binary_data__('velx', chunk_number)
dset = tmp**2
@ -259,7 +265,7 @@ class Amun:
dset += tmp**2
tmp = self.__read_binary_data__('velz', chunk_number)
dset += tmp**2
dset = np.sqrt(dset)
dset = numpy.sqrt(dset)
elif dataset == 'vvec':
dset = [self.__read_binary_data__('velx', chunk_number), \
self.__read_binary_data__('vely', chunk_number), \
@ -271,7 +277,7 @@ class Amun:
dset += tmp**2
tmp = self.__read_binary_data__('magz', chunk_number)
dset += tmp**2
dset = np.sqrt(dset)
dset = numpy.sqrt(dset)
elif dataset == 'bvec':
dset = [self.__read_binary_data__('magx', chunk_number), \
self.__read_binary_data__('magy', chunk_number), \
@ -302,10 +308,10 @@ class Amun:
dset = vy * bz - vz * by
if self.attributes['resistivity'] > 0:
if self.attributes['ndims'] == 3:
tmp = (np.roll(by, 1, axis=0) - np.roll(by, -1, axis=0))
tmp += (np.roll(bz, -1, axis=1) - np.roll(bz, 1, axis=1))
tmp = (numpy.roll(by, 1, axis=0) - numpy.roll(by, -1, axis=0))
tmp += (numpy.roll(bz, -1, axis=1) - numpy.roll(bz, 1, axis=1))
else:
tmp = (np.roll(bz, -1, axis=0) - np.roll(bz, 1, axis=0))
tmp = (numpy.roll(bz, -1, axis=0) - numpy.roll(bz, 1, axis=0))
for p in range(self.chunks[chunk_number]['dblocks']):
tmp[...,p] /= self.cell_size[self.chunks[chunk_number]['levels'][p]]
@ -318,10 +324,10 @@ class Amun:
dset = vz * bx - vx * bz
if self.attributes['resistivity'] > 0:
if self.attributes['ndims'] == 3:
tmp = (np.roll(bx, -1, axis=0) - np.roll(bx, 1, axis=0))
tmp += (np.roll(bz, 1, axis=2) - np.roll(bz, -1, axis=2))
tmp = (numpy.roll(bx, -1, axis=0) - numpy.roll(bx, 1, axis=0))
tmp += (numpy.roll(bz, 1, axis=2) - numpy.roll(bz, -1, axis=2))
else:
tmp = (np.roll(bz, 1, axis=1) - np.roll(bz, -1, axis=1))
tmp = (numpy.roll(bz, 1, axis=1) - numpy.roll(bz, -1, axis=1))
for p in range(self.chunks[chunk_number]['dblocks']):
tmp[...,p] /= self.cell_size[self.chunks[chunk_number]['levels'][p]]
@ -334,11 +340,11 @@ class Amun:
dset = vx * by - vy * bx
if self.attributes['resistivity'] > 0:
if self.attributes['ndims'] == 3:
tmp = (np.roll(bx, 1, axis=1) - np.roll(bx, -1, axis=1))
tmp += (np.roll(by, -1, axis=2) - np.roll(by, 1, axis=2))
tmp = (numpy.roll(bx, 1, axis=1) - numpy.roll(bx, -1, axis=1))
tmp += (numpy.roll(by, -1, axis=2) - numpy.roll(by, 1, axis=2))
else:
tmp = (np.roll(bx, 1, axis=0) - np.roll(bx, -1, axis=0))
tmp += (np.roll(by, -1, axis=1) - np.roll(by, 1, axis=1))
tmp = (numpy.roll(bx, 1, axis=0) - numpy.roll(bx, -1, axis=0))
tmp += (numpy.roll(by, -1, axis=1) - numpy.roll(by, 1, axis=1))
for p in range(self.chunks[chunk_number]['dblocks']):
tmp[...,p] /= self.cell_size[self.chunks[chunk_number]['levels'][p]]
@ -351,10 +357,10 @@ class Amun:
dtmp = v1 * b2 - v2 * b1
if self.attributes['resistivity'] > 0:
if self.attributes['ndims'] == 3:
tmp = (np.roll(b1, 1, axis=0) - np.roll(b1, -1, axis=0))
tmp += (np.roll(b2, -1, axis=1) - np.roll(b2, 1, axis=1))
tmp = (numpy.roll(b1, 1, axis=0) - numpy.roll(b1, -1, axis=0))
tmp += (numpy.roll(b2, -1, axis=1) - numpy.roll(b2, 1, axis=1))
else:
tmp = (np.roll(b2, -1, axis=0) - np.roll(b2, 1, axis=0))
tmp = (numpy.roll(b2, -1, axis=0) - numpy.roll(b2, 1, axis=0))
for p in range(self.chunks[chunk_number]['dblocks']):
tmp[...,p] /= self.cell_size[self.chunks[chunk_number]['levels'][p]]
@ -366,10 +372,10 @@ class Amun:
dtmp = v2 * b1 - v1 * b2
if self.attributes['resistivity'] > 0:
if self.attributes['ndims'] == 3:
tmp = (np.roll(b1, -1, axis=0) - np.roll(b1, 1, axis=0))
tmp += (np.roll(b2, 1, axis=2) - np.roll(b2, -1, axis=2))
tmp = (numpy.roll(b1, -1, axis=0) - numpy.roll(b1, 1, axis=0))
tmp += (numpy.roll(b2, 1, axis=2) - numpy.roll(b2, -1, axis=2))
else:
tmp = (np.roll(b2, 1, axis=1) - np.roll(b2, -1, axis=1))
tmp = (numpy.roll(b2, 1, axis=1) - numpy.roll(b2, -1, axis=1))
for p in range(self.chunks[chunk_number]['dblocks']):
tmp[...,p] /= self.cell_size[self.chunks[chunk_number]['levels'][p]]
@ -381,17 +387,17 @@ class Amun:
dtmp = v1 * b2 - v2 * b1
if self.attributes['resistivity'] > 0:
if self.attributes['ndims'] == 3:
tmp = (np.roll(b1, 1, axis=1) - np.roll(b1, -1, axis=1))
tmp += (np.roll(b2, -1, axis=2) - np.roll(b2, 1, axis=2))
tmp = (numpy.roll(b1, 1, axis=1) - numpy.roll(b1, -1, axis=1))
tmp += (numpy.roll(b2, -1, axis=2) - numpy.roll(b2, 1, axis=2))
else:
tmp = (np.roll(b1, 1, axis=0) - np.roll(b1, -1, axis=0))
tmp += (np.roll(b2, -1, axis=1) - np.roll(b2, 1, axis=1))
tmp = (numpy.roll(b1, 1, axis=0) - numpy.roll(b1, -1, axis=0))
tmp += (numpy.roll(b2, -1, axis=1) - numpy.roll(b2, 1, axis=1))
for p in range(self.chunks[chunk_number]['dblocks']):
tmp[...,p] /= self.cell_size[self.chunks[chunk_number]['levels'][p]]
dtmp -= 0.5 * self.attributes['resistivity'] * tmp
dset += dtmp**2
dset = np.sqrt(dset)
dset = numpy.sqrt(dset)
elif dataset == 'evec':
b1 = self.__read_binary_data__('magy', chunk_number)
b2 = self.__read_binary_data__('magz', chunk_number)
@ -400,10 +406,10 @@ class Amun:
wx = v1 * b2 - v2 * b1
if self.attributes['resistivity'] > 0:
if self.attributes['ndims'] == 3:
tmp = (np.roll(b1, 1, axis=0) - np.roll(b1, -1, axis=0))
tmp += (np.roll(b2, -1, axis=1) - np.roll(b2, 1, axis=1))
tmp = (numpy.roll(b1, 1, axis=0) - numpy.roll(b1, -1, axis=0))
tmp += (numpy.roll(b2, -1, axis=1) - numpy.roll(b2, 1, axis=1))
else:
tmp = (np.roll(b2, -1, axis=0) - np.roll(b2, 1, axis=0))
tmp = (numpy.roll(b2, -1, axis=0) - numpy.roll(b2, 1, axis=0))
for p in range(self.chunks[chunk_number]['dblocks']):
tmp[...,p] /= self.cell_size[self.chunks[chunk_number]['levels'][p]]
@ -414,10 +420,10 @@ class Amun:
wy = v2 * b1 - v1 * b2
if self.attributes['resistivity'] > 0:
if self.attributes['ndims'] == 3:
tmp = (np.roll(b1, -1, axis=0) - np.roll(b1, 1, axis=0))
tmp += (np.roll(b2, 1, axis=2) - np.roll(b2, -1, axis=2))
tmp = (numpy.roll(b1, -1, axis=0) - numpy.roll(b1, 1, axis=0))
tmp += (numpy.roll(b2, 1, axis=2) - numpy.roll(b2, -1, axis=2))
else:
tmp = (np.roll(b2, 1, axis=1) - np.roll(b2, -1, axis=1))
tmp = (numpy.roll(b2, 1, axis=1) - numpy.roll(b2, -1, axis=1))
for p in range(self.chunks[chunk_number]['dblocks']):
tmp[...,p] /= self.cell_size[self.chunks[chunk_number]['levels'][p]]
@ -428,11 +434,11 @@ class Amun:
wz = v1 * b2 - v2 * b1
if self.attributes['resistivity'] > 0:
if self.attributes['ndims'] == 3:
tmp = (np.roll(b1, 1, axis=1) - np.roll(b1, -1, axis=1))
tmp += (np.roll(b2, -1, axis=2) - np.roll(b2, 1, axis=2))
tmp = (numpy.roll(b1, 1, axis=1) - numpy.roll(b1, -1, axis=1))
tmp += (numpy.roll(b2, -1, axis=2) - numpy.roll(b2, 1, axis=2))
else:
tmp = (np.roll(b1, 1, axis=0) - np.roll(b1, -1, axis=0))
tmp += (np.roll(b2, -1, axis=1) - np.roll(b2, 1, axis=1))
tmp = (numpy.roll(b1, 1, axis=0) - numpy.roll(b1, -1, axis=0))
tmp += (numpy.roll(b2, -1, axis=1) - numpy.roll(b2, 1, axis=1))
for p in range(self.chunks[chunk_number]['dblocks']):
tmp[...,p] /= self.cell_size[self.chunks[chunk_number]['levels'][p]]
@ -470,21 +476,21 @@ class Amun:
dset += tmp**2
tmp = self.__read_binary_data__('velz', chunk_number)
dset += tmp**2
dset = 1.0 / np.sqrt(1.0 - dset)
dset = 1.0 / numpy.sqrt(1.0 - dset)
elif dataset == 'divv':
p = self.attributes['ndims'] - 1
tmp = self.__read_binary_data__('velx', chunk_number)
dset = (np.roll(tmp, -1, axis=p) \
- np.roll(tmp, 1, axis=p))
dset = (numpy.roll(tmp, -1, axis=p) \
- numpy.roll(tmp, 1, axis=p))
p -= 1
tmp = self.__read_binary_data__('vely', chunk_number)
dset += (np.roll(tmp, -1, axis=p) \
- np.roll(tmp, 1, axis=p))
dset += (numpy.roll(tmp, -1, axis=p) \
- numpy.roll(tmp, 1, axis=p))
p -= 1
if p >= 0:
tmp = self.__read_binary_data__('velz', chunk_number)
dset += (np.roll(tmp, -1, axis=p) \
- np.roll(tmp, 0, axis=p))
dset += (numpy.roll(tmp, -1, axis=p) \
- numpy.roll(tmp, 0, axis=p))
dset *= 0.5
for p in range(self.chunks[chunk_number]['dblocks']):
@ -492,15 +498,15 @@ class Amun:
elif dataset == 'vorx':
if self.attributes['ndims'] == 3:
tmp = self.__read_binary_data__('vely', chunk_number)
dset = (np.roll(tmp, 1, axis=0) \
- np.roll(tmp, -1, axis=0))
dset = (numpy.roll(tmp, 1, axis=0) \
- numpy.roll(tmp, -1, axis=0))
tmp = self.__read_binary_data__('velz', chunk_number)
dset += (np.roll(tmp, -1, axis=1) \
- np.roll(tmp, 1, axis=1))
dset += (numpy.roll(tmp, -1, axis=1) \
- numpy.roll(tmp, 1, axis=1))
else:
tmp = self.__read_binary_data__('velz', chunk_number)
dset = (np.roll(tmp, -1, axis=0) \
- np.roll(tmp, 1, axis=0))
dset = (numpy.roll(tmp, -1, axis=0) \
- numpy.roll(tmp, 1, axis=0))
dset *= 0.5
for p in range(self.chunks[chunk_number]['dblocks']):
@ -508,15 +514,15 @@ class Amun:
elif dataset == 'vory':
if self.attributes['ndims'] == 3:
tmp = self.__read_binary_data__('velx', chunk_number)
dset = (np.roll(tmp, -1, axis=0) \
- np.roll(tmp, 1, axis=0))
dset = (numpy.roll(tmp, -1, axis=0) \
- numpy.roll(tmp, 1, axis=0))
tmp = self.__read_binary_data__('velz', chunk_number)
dset += (np.roll(tmp, 1, axis=2) \
- np.roll(tmp, -1, axis=2))
dset += (numpy.roll(tmp, 1, axis=2) \
- numpy.roll(tmp, -1, axis=2))
else:
tmp = self.__read_binary_data__('velz', chunk_number)
dset = (np.roll(tmp, 1, axis=1) \
- np.roll(tmp, -1, axis=1))
dset = (numpy.roll(tmp, 1, axis=1) \
- numpy.roll(tmp, -1, axis=1))
dset *= 0.5
for p in range(self.chunks[chunk_number]['dblocks']):
@ -524,18 +530,18 @@ class Amun:
elif dataset == 'vorz':
if self.attributes['ndims'] == 3:
tmp = self.__read_binary_data__('velx', chunk_number)
dset = (np.roll(tmp, 1, axis=1) \
- np.roll(tmp, -1, axis=1))
dset = (numpy.roll(tmp, 1, axis=1) \
- numpy.roll(tmp, -1, axis=1))
tmp = self.__read_binary_data__('vely', chunk_number)
dset += (np.roll(tmp, -1, axis=2) \
- np.roll(tmp, 1, axis=2))
dset += (numpy.roll(tmp, -1, axis=2) \
- numpy.roll(tmp, 1, axis=2))
else:
tmp = self.__read_binary_data__('velx', chunk_number)
dset = (np.roll(tmp, 1, axis=0) \
- np.roll(tmp, -1, axis=0))
dset = (numpy.roll(tmp, 1, axis=0) \
- numpy.roll(tmp, -1, axis=0))
tmp = self.__read_binary_data__('vely', chunk_number)
dset += (np.roll(tmp, -1, axis=1) \
- np.roll(tmp, 1, axis=1))
dset += (numpy.roll(tmp, -1, axis=1) \
- numpy.roll(tmp, 1, axis=1))
dset *= 0.5
for p in range(self.chunks[chunk_number]['dblocks']):
@ -543,32 +549,32 @@ class Amun:
elif dataset == 'wvec':
if self.attributes['ndims'] == 3:
tmp = self.__read_binary_data__('velx', chunk_number)
wy = (np.roll(tmp, -1, axis=0) \
- np.roll(tmp, 1, axis=0))
wz = (np.roll(tmp, 1, axis=1) \
- np.roll(tmp, -1, axis=1))
wy = (numpy.roll(tmp, -1, axis=0) \
- numpy.roll(tmp, 1, axis=0))
wz = (numpy.roll(tmp, 1, axis=1) \
- numpy.roll(tmp, -1, axis=1))
tmp = self.__read_binary_data__('vely', chunk_number)
wz += (np.roll(tmp, -1, axis=2) \
- np.roll(tmp, 1, axis=2))
wx = (np.roll(tmp, 1, axis=0) \
- np.roll(tmp, -1, axis=0))
wz += (numpy.roll(tmp, -1, axis=2) \
- numpy.roll(tmp, 1, axis=2))
wx = (numpy.roll(tmp, 1, axis=0) \
- numpy.roll(tmp, -1, axis=0))
tmp = self.__read_binary_data__('velz', chunk_number)
wx += (np.roll(tmp, -1, axis=1) \
- np.roll(tmp, 1, axis=1))
wy += (np.roll(tmp, 1, axis=2) \
- np.roll(tmp, -1, axis=2))
wx += (numpy.roll(tmp, -1, axis=1) \
- numpy.roll(tmp, 1, axis=1))
wy += (numpy.roll(tmp, 1, axis=2) \
- numpy.roll(tmp, -1, axis=2))
else:
tmp = self.__read_binary_data__('velx', chunk_number)
wz = (np.roll(tmp, 1, axis=0) \
- np.roll(tmp, -1, axis=0))
wz = (numpy.roll(tmp, 1, axis=0) \
- numpy.roll(tmp, -1, axis=0))
tmp = self.__read_binary_data__('vely', chunk_number)
wz += (np.roll(tmp, -1, axis=1) \
- np.roll(tmp, 1, axis=1))
wz += (numpy.roll(tmp, -1, axis=1) \
- numpy.roll(tmp, 1, axis=1))
tmp = self.__read_binary_data__('velz', chunk_number)
wx = (np.roll(tmp, -1, axis=0) \
- np.roll(tmp, 1, axis=0))
wy = (np.roll(tmp, 1, axis=1) \
- np.roll(tmp, -1, axis=1))
wx = (numpy.roll(tmp, -1, axis=0) \
- numpy.roll(tmp, 1, axis=0))
wy = (numpy.roll(tmp, 1, axis=1) \
- numpy.roll(tmp, -1, axis=1))
for p in range(self.chunks[chunk_number]['dblocks']):
h = 2 * self.cell_size[self.chunks[chunk_number]['levels'][p]]
@ -580,47 +586,47 @@ class Amun:
elif dataset == 'vort':
if self.attributes['ndims'] == 3:
tmp = self.__read_binary_data__('velx', chunk_number)
wy = (np.roll(tmp, -1, axis=0) \
- np.roll(tmp, 1, axis=0))
wz = (np.roll(tmp, 1, axis=1) \
- np.roll(tmp, -1, axis=1))
wy = (numpy.roll(tmp, -1, axis=0) \
- numpy.roll(tmp, 1, axis=0))
wz = (numpy.roll(tmp, 1, axis=1) \
- numpy.roll(tmp, -1, axis=1))
tmp = self.__read_binary_data__('vely', chunk_number)
wz += (np.roll(tmp, -1, axis=2) \
- np.roll(tmp, 1, axis=2))
wx = (np.roll(tmp, 1, axis=0) \
- np.roll(tmp, -1, axis=0))
wz += (numpy.roll(tmp, -1, axis=2) \
- numpy.roll(tmp, 1, axis=2))
wx = (numpy.roll(tmp, 1, axis=0) \
- numpy.roll(tmp, -1, axis=0))
tmp = self.__read_binary_data__('velz', chunk_number)
wx += (np.roll(tmp, -1, axis=1) \
- np.roll(tmp, 1, axis=1))
wy += (np.roll(tmp, 1, axis=2) \
- np.roll(tmp, -1, axis=2))
wx += (numpy.roll(tmp, -1, axis=1) \
- numpy.roll(tmp, 1, axis=1))
wy += (numpy.roll(tmp, 1, axis=2) \
- numpy.roll(tmp, -1, axis=2))
else:
tmp = self.__read_binary_data__('velx', chunk_number)
wz = (np.roll(tmp, 1, axis=0) \
- np.roll(tmp, -1, axis=0))
wz = (numpy.roll(tmp, 1, axis=0) \
- numpy.roll(tmp, -1, axis=0))
tmp = self.__read_binary_data__('vely', chunk_number)
wz += (np.roll(tmp, -1, axis=1) \
- np.roll(tmp, 1, axis=1))
wz += (numpy.roll(tmp, -1, axis=1) \
- numpy.roll(tmp, 1, axis=1))
tmp = self.__read_binary_data__('velz', chunk_number)
wx = (np.roll(tmp, -1, axis=0) \
- np.roll(tmp, 1, axis=0))
wy = (np.roll(tmp, 1, axis=1) \
- np.roll(tmp, -1, axis=1))
wx = (numpy.roll(tmp, -1, axis=0) \
- numpy.roll(tmp, 1, axis=0))
wy = (numpy.roll(tmp, 1, axis=1) \
- numpy.roll(tmp, -1, axis=1))
dset = 0.5 * np.sqrt(wx**2 + wy**2 + wz**2)
dset = 0.5 * numpy.sqrt(wx**2 + wy**2 + wz**2)
for p in range(self.chunks[chunk_number]['dblocks']):
dset[...,p] /= self.cell_size[self.chunks[chunk_number]['levels'][p]]
elif dataset == 'dvxdx':
p = self.attributes['ndims'] - 1
tmp = self.__read_binary_data__('velx', chunk_number)
dset = np.roll(tmp, -1, axis=p) - np.roll(tmp, 1, axis=p)
dset = numpy.roll(tmp, -1, axis=p) - numpy.roll(tmp, 1, axis=p)
dset *= 0.5
for p in range(self.chunks[chunk_number]['dblocks']):
dset[...,p] /= self.cell_size[self.chunks[chunk_number]['levels'][p]]
elif dataset == 'dvxdy':
p = self.attributes['ndims'] - 2
tmp = self.__read_binary_data__('velx', chunk_number)
dset = np.roll(tmp, -1, axis=p) - np.roll(tmp, 1, axis=p)
dset = numpy.roll(tmp, -1, axis=p) - numpy.roll(tmp, 1, axis=p)
dset *= 0.5
for p in range(self.chunks[chunk_number]['dblocks']):
dset[...,p] /= self.cell_size[self.chunks[chunk_number]['levels'][p]]
@ -628,23 +634,23 @@ class Amun:
p = self.attributes['ndims'] - 3
tmp = self.__read_binary_data__('velx', chunk_number)
if p >= 0:
dset = np.roll(tmp, -1, axis=p) - np.roll(tmp, 1, axis=p)
dset = numpy.roll(tmp, -1, axis=p) - numpy.roll(tmp, 1, axis=p)
dset *= 0.5
for p in range(self.chunks[chunk_number]['dblocks']):
dset[...,p] /= self.cell_size[self.chunks[chunk_number]['levels'][p]]
else:
dset = np.zeros_like(tmp)
dset = numpy.zeros_like(tmp)
elif dataset == 'dvydx':
p = self.attributes['ndims'] - 1
tmp = self.__read_binary_data__('vely', chunk_number)
dset = np.roll(tmp, -1, axis=p) - np.roll(tmp, 1, axis=p)
dset = numpy.roll(tmp, -1, axis=p) - numpy.roll(tmp, 1, axis=p)
dset *= 0.5
for p in range(self.chunks[chunk_number]['dblocks']):
dset[...,p] /= self.cell_size[self.chunks[chunk_number]['levels'][p]]
elif dataset == 'dvydy':
p = self.attributes['ndims'] - 2
tmp = self.__read_binary_data__('vely', chunk_number)
dset = np.roll(tmp, -1, axis=p) - np.roll(tmp, 1, axis=p)
dset = numpy.roll(tmp, -1, axis=p) - numpy.roll(tmp, 1, axis=p)
dset *= 0.5
for p in range(self.chunks[chunk_number]['dblocks']):
dset[...,p] /= self.cell_size[self.chunks[chunk_number]['levels'][p]]
@ -652,23 +658,23 @@ class Amun:
p = self.attributes['ndims'] - 3
tmp = self.__read_binary_data__('vely', chunk_number)
if p >= 0:
dset = np.roll(tmp, -1, axis=p) - np.roll(tmp, 1, axis=p)
dset = numpy.roll(tmp, -1, axis=p) - numpy.roll(tmp, 1, axis=p)
dset *= 0.5
for p in range(self.chunks[chunk_number]['dblocks']):
dset[...,p] /= self.cell_size[self.chunks[chunk_number]['levels'][p]]
else:
dset = np.zeros_like(tmp)
dset = numpy.zeros_like(tmp)
elif dataset == 'dvzdx':
p = self.attributes['ndims'] - 1
tmp = self.__read_binary_data__('velz', chunk_number)
dset = np.roll(tmp, -1, axis=p) - np.roll(tmp, 1, axis=p)
dset = numpy.roll(tmp, -1, axis=p) - numpy.roll(tmp, 1, axis=p)
dset *= 0.5
for p in range(self.chunks[chunk_number]['dblocks']):
dset[...,p] /= self.cell_size[self.chunks[chunk_number]['levels'][p]]
elif dataset == 'dvzdy':
p = self.attributes['ndims'] - 2
tmp = self.__read_binary_data__('velz', chunk_number)
dset = np.roll(tmp, -1, axis=p) - np.roll(tmp, 1, axis=p)
dset = numpy.roll(tmp, -1, axis=p) - numpy.roll(tmp, 1, axis=p)
dset *= 0.5
for p in range(self.chunks[chunk_number]['dblocks']):
dset[...,p] /= self.cell_size[self.chunks[chunk_number]['levels'][p]]
@ -676,26 +682,26 @@ class Amun:
p = self.attributes['ndims'] - 3
tmp = self.__read_binary_data__('velz', chunk_number)
if p >= 0:
dset = np.roll(tmp, -1, axis=p) - np.roll(tmp, 1, axis=p)
dset = numpy.roll(tmp, -1, axis=p) - numpy.roll(tmp, 1, axis=p)
dset *= 0.5
for p in range(self.chunks[chunk_number]['dblocks']):
dset[...,p] /= self.cell_size[self.chunks[chunk_number]['levels'][p]]
else:
dset = np.zeros_like(tmp)
dset = numpy.zeros_like(tmp)
elif dataset == 'divb':
p = self.attributes['ndims'] - 1
tmp = self.__read_binary_data__('magx', chunk_number)
dset = (np.roll(tmp, -1, axis=p) \
- np.roll(tmp, 1, axis=p))
dset = (numpy.roll(tmp, -1, axis=p) \
- numpy.roll(tmp, 1, axis=p))
p -= 1
tmp = self.__read_binary_data__('magy', chunk_number)
dset += (np.roll(tmp, -1, axis=p) \
- np.roll(tmp, 1, axis=p))
dset += (numpy.roll(tmp, -1, axis=p) \
- numpy.roll(tmp, 1, axis=p))
p -= 1
if p >= 0:
tmp = self.__read_binary_data__('magz', chunk_number)
dset += (np.roll(tmp, -1, axis=p) \
- np.roll(tmp, 0, axis=p))
dset += (numpy.roll(tmp, -1, axis=p) \
- numpy.roll(tmp, 0, axis=p))
dset *= 0.5
for p in range(self.chunks[chunk_number]['dblocks']):
@ -703,15 +709,15 @@ class Amun:
elif dataset == 'curx':
if self.attributes['ndims'] == 3:
tmp = self.__read_binary_data__('magy', chunk_number)
dset = (np.roll(tmp, 1, axis=0) \
- np.roll(tmp, -1, axis=0))
dset = (numpy.roll(tmp, 1, axis=0) \
- numpy.roll(tmp, -1, axis=0))
tmp = self.__read_binary_data__('magz', chunk_number)
dset += (np.roll(tmp, -1, axis=1) \
- np.roll(tmp, 1, axis=1))
dset += (numpy.roll(tmp, -1, axis=1) \
- numpy.roll(tmp, 1, axis=1))
else:
tmp = self.__read_binary_data__('magz', chunk_number)
dset = (np.roll(tmp, -1, axis=0) \
- np.roll(tmp, 1, axis=0))
dset = (numpy.roll(tmp, -1, axis=0) \
- numpy.roll(tmp, 1, axis=0))
dset *= 0.5
for p in range(self.chunks[chunk_number]['dblocks']):
@ -719,15 +725,15 @@ class Amun:
elif dataset == 'cury':
if self.attributes['ndims'] == 3:
tmp = self.__read_binary_data__('magx', chunk_number)
dset = (np.roll(tmp, -1, axis=0) \
- np.roll(tmp, 1, axis=0))
dset = (numpy.roll(tmp, -1, axis=0) \
- numpy.roll(tmp, 1, axis=0))
tmp = self.__read_binary_data__('magz', chunk_number)
dset += (np.roll(tmp, 1, axis=2) \
- np.roll(tmp, -1, axis=2))
dset += (numpy.roll(tmp, 1, axis=2) \
- numpy.roll(tmp, -1, axis=2))
else:
tmp = self.__read_binary_data__('magz', chunk_number)
dset = (np.roll(tmp, 1, axis=1) \
- np.roll(tmp, -1, axis=1))
dset = (numpy.roll(tmp, 1, axis=1) \
- numpy.roll(tmp, -1, axis=1))
dset *= 0.5
for p in range(self.chunks[chunk_number]['dblocks']):
@ -735,18 +741,18 @@ class Amun:
elif dataset == 'curz':
if self.attributes['ndims'] == 3:
tmp = self.__read_binary_data__('magx', chunk_number)
dset = (np.roll(tmp, 1, axis=1) \
- np.roll(tmp, -1, axis=1))
dset = (numpy.roll(tmp, 1, axis=1) \
- numpy.roll(tmp, -1, axis=1))
tmp = self.__read_binary_data__('magy', chunk_number)
dset += (np.roll(tmp, -1, axis=2) \
- np.roll(tmp, 1, axis=2))
dset += (numpy.roll(tmp, -1, axis=2) \
- numpy.roll(tmp, 1, axis=2))
else:
tmp = self.__read_binary_data__('magx', chunk_number)
dset = (np.roll(tmp, 1, axis=0) \
- np.roll(tmp, -1, axis=0))
dset = (numpy.roll(tmp, 1, axis=0) \
- numpy.roll(tmp, -1, axis=0))
tmp = self.__read_binary_data__('magy', chunk_number)
dset += (np.roll(tmp, -1, axis=1) \
- np.roll(tmp, 1, axis=1))
dset += (numpy.roll(tmp, -1, axis=1) \
- numpy.roll(tmp, 1, axis=1))
dset *= 0.5
for p in range(self.chunks[chunk_number]['dblocks']):
@ -754,32 +760,32 @@ class Amun:
elif dataset == 'jvec':
if self.attributes['ndims'] == 3:
tmp = self.__read_binary_data__('magx', chunk_number)
wy = (np.roll(tmp, -1, axis=0) \
- np.roll(tmp, 1, axis=0))
wz = (np.roll(tmp, 1, axis=1) \
- np.roll(tmp, -1, axis=1))
wy = (numpy.roll(tmp, -1, axis=0) \
- numpy.roll(tmp, 1, axis=0))
wz = (numpy.roll(tmp, 1, axis=1) \
- numpy.roll(tmp, -1, axis=1))
tmp = self.__read_binary_data__('magy', chunk_number)
wz += (np.roll(tmp, -1, axis=2) \
- np.roll(tmp, 1, axis=2))
wx = (np.roll(tmp, 1, axis=0) \
- np.roll(tmp, -1, axis=0))
wz += (numpy.roll(tmp, -1, axis=2) \
- numpy.roll(tmp, 1, axis=2))
wx = (numpy.roll(tmp, 1, axis=0) \
- numpy.roll(tmp, -1, axis=0))
tmp = self.__read_binary_data__('magz', chunk_number)
wx += (np.roll(tmp, -1, axis=1) \
- np.roll(tmp, 1, axis=1))
wy += (np.roll(tmp, 1, axis=2) \
- np.roll(tmp, -1, axis=2))
wx += (numpy.roll(tmp, -1, axis=1) \
- numpy.roll(tmp, 1, axis=1))
wy += (numpy.roll(tmp, 1, axis=2) \
- numpy.roll(tmp, -1, axis=2))
else:
tmp = self.__read_binary_data__('magx', chunk_number)
wz = (np.roll(tmp, 1, axis=0) \
- np.roll(tmp, -1, axis=0))
wz = (numpy.roll(tmp, 1, axis=0) \
- numpy.roll(tmp, -1, axis=0))
tmp = self.__read_binary_data__('magy', chunk_number)
wz += (np.roll(tmp, -1, axis=1) \
- np.roll(tmp, 1, axis=1))
wz += (numpy.roll(tmp, -1, axis=1) \
- numpy.roll(tmp, 1, axis=1))
tmp = self.__read_binary_data__('magz', chunk_number)
wx = (np.roll(tmp, -1, axis=0) \
- np.roll(tmp, 1, axis=0))
wy = (np.roll(tmp, 1, axis=1) \
- np.roll(tmp, -1, axis=1))
wx = (numpy.roll(tmp, -1, axis=0) \
- numpy.roll(tmp, 1, axis=0))
wy = (numpy.roll(tmp, 1, axis=1) \
- numpy.roll(tmp, -1, axis=1))
for p in range(self.chunks[chunk_number]['dblocks']):
h = 2 * self.cell_size[self.chunks[chunk_number]['levels'][p]]
@ -791,47 +797,47 @@ class Amun:
elif dataset == 'curr':
if self.attributes['ndims'] == 3:
tmp = self.__read_binary_data__('magx', chunk_number)
wy = (np.roll(tmp, -1, axis=0) \
- np.roll(tmp, 1, axis=0))
wz = (np.roll(tmp, 1, axis=1) \
- np.roll(tmp, -1, axis=1))
wy = (numpy.roll(tmp, -1, axis=0) \
- numpy.roll(tmp, 1, axis=0))
wz = (numpy.roll(tmp, 1, axis=1) \
- numpy.roll(tmp, -1, axis=1))
tmp = self.__read_binary_data__('magy', chunk_number)
wz += (np.roll(tmp, -1, axis=2) \
- np.roll(tmp, 1, axis=2))
wx = (np.roll(tmp, 1, axis=0) \
- np.roll(tmp, -1, axis=0))
wz += (numpy.roll(tmp, -1, axis=2) \
- numpy.roll(tmp, 1, axis=2))
wx = (numpy.roll(tmp, 1, axis=0) \
- numpy.roll(tmp, -1, axis=0))
tmp = self.__read_binary_data__('magz', chunk_number)
wx += (np.roll(tmp, -1, axis=1) \
- np.roll(tmp, 1, axis=1))
wy += (np.roll(tmp, 1, axis=2) \
- np.roll(tmp, -1, axis=2))
wx += (numpy.roll(tmp, -1, axis=1) \
- numpy.roll(tmp, 1, axis=1))
wy += (numpy.roll(tmp, 1, axis=2) \
- numpy.roll(tmp, -1, axis=2))
else:
tmp = self.__read_binary_data__('magx', chunk_number)
wz = (np.roll(tmp, 1, axis=0) \
- np.roll(tmp, -1, axis=0))
wz = (numpy.roll(tmp, 1, axis=0) \
- numpy.roll(tmp, -1, axis=0))
tmp = self.__read_binary_data__('magy', chunk_number)
wz += (np.roll(tmp, -1, axis=1) \
- np.roll(tmp, 1, axis=1))
wz += (numpy.roll(tmp, -1, axis=1) \
- numpy.roll(tmp, 1, axis=1))
tmp = self.__read_binary_data__('magz', chunk_number)
wx = (np.roll(tmp, -1, axis=0) \
- np.roll(tmp, 1, axis=0))
wy = (np.roll(tmp, 1, axis=1) \
- np.roll(tmp, -1, axis=1))
wx = (numpy.roll(tmp, -1, axis=0) \
- numpy.roll(tmp, 1, axis=0))
wy = (numpy.roll(tmp, 1, axis=1) \
- numpy.roll(tmp, -1, axis=1))
dset = 0.5 * np.sqrt(wx**2 + wy**2 + wz**2)
dset = 0.5 * numpy.sqrt(wx**2 + wy**2 + wz**2)
for p in range(self.chunks[chunk_number]['dblocks']):
dset[...,p] /= self.cell_size[self.chunks[chunk_number]['levels'][p]]
elif dataset == 'dbxdx':
p = self.attributes['ndims'] - 1
tmp = self.__read_binary_data__('magx', chunk_number)
dset = np.roll(tmp, -1, axis=p) - np.roll(tmp, 1, axis=p)
dset = numpy.roll(tmp, -1, axis=p) - numpy.roll(tmp, 1, axis=p)
dset *= 0.5
for p in range(self.chunks[chunk_number]['dblocks']):
dset[...,p] /= self.cell_size[self.chunks[chunk_number]['levels'][p]]
elif dataset == 'dbxdy':
p = self.attributes['ndims'] - 2
tmp = self.__read_binary_data__('magx', chunk_number)
dset = np.roll(tmp, -1, axis=p) - np.roll(tmp, 1, axis=p)
dset = numpy.roll(tmp, -1, axis=p) - numpy.roll(tmp, 1, axis=p)
dset *= 0.5
for p in range(self.chunks[chunk_number]['dblocks']):
dset[...,p] /= self.cell_size[self.chunks[chunk_number]['levels'][p]]
@ -839,23 +845,23 @@ class Amun:
p = self.attributes['ndims'] - 3
tmp = self.__read_binary_data__('magx', chunk_number)
if p >= 0:
dset = np.roll(tmp, -1, axis=p) - np.roll(tmp, 1, axis=p)
dset = numpy.roll(tmp, -1, axis=p) - numpy.roll(tmp, 1, axis=p)
dset *= 0.5
for p in range(self.chunks[chunk_number]['dblocks']):
dset[...,p] /= self.cell_size[self.chunks[chunk_number]['levels'][p]]
else:
dset = np.zeros_like(tmp)
dset = numpy.zeros_like(tmp)
elif dataset == 'dbydx':
p = self.attributes['ndims'] - 1
tmp = self.__read_binary_data__('magy', chunk_number)
dset = np.roll(tmp, -1, axis=p) - np.roll(tmp, 1, axis=p)
dset = numpy.roll(tmp, -1, axis=p) - numpy.roll(tmp, 1, axis=p)
dset *= 0.5
for p in range(self.chunks[chunk_number]['dblocks']):
dset[...,p] /= self.cell_size[self.chunks[chunk_number]['levels'][p]]
elif dataset == 'dbydy':
p = self.attributes['ndims'] - 2
tmp = self.__read_binary_data__('magy', chunk_number)
dset = np.roll(tmp, -1, axis=p) - np.roll(tmp, 1, axis=p)
dset = numpy.roll(tmp, -1, axis=p) - numpy.roll(tmp, 1, axis=p)
dset *= 0.5
for p in range(self.chunks[chunk_number]['dblocks']):
dset[...,p] /= self.cell_size[self.chunks[chunk_number]['levels'][p]]
@ -863,23 +869,23 @@ class Amun:
p = self.attributes['ndims'] - 3
tmp = self.__read_binary_data__('magy', chunk_number)
if p >= 0:
dset = np.roll(tmp, -1, axis=p) - np.roll(tmp, 1, axis=p)
dset = numpy.roll(tmp, -1, axis=p) - numpy.roll(tmp, 1, axis=p)
dset *= 0.5
for p in range(self.chunks[chunk_number]['dblocks']):
dset[...,p] /= self.cell_size[self.chunks[chunk_number]['levels'][p]]
else:
dset = np.zeros_like(tmp)
dset = numpy.zeros_like(tmp)
elif dataset == 'dbzdx':
p = self.attributes['ndims'] - 1
tmp = self.__read_binary_data__('magz', chunk_number)
dset = np.roll(tmp, -1, axis=p) - np.roll(tmp, 1, axis=p)
dset = numpy.roll(tmp, -1, axis=p) - numpy.roll(tmp, 1, axis=p)
dset *= 0.5
for p in range(self.chunks[chunk_number]['dblocks']):
dset[...,p] /= self.cell_size[self.chunks[chunk_number]['levels'][p]]
elif dataset == 'dbzdy':
p = self.attributes['ndims'] - 2
tmp = self.__read_binary_data__('magz', chunk_number)
dset = np.roll(tmp, -1, axis=p) - np.roll(tmp, 1, axis=p)
dset = numpy.roll(tmp, -1, axis=p) - numpy.roll(tmp, 1, axis=p)
dset *= 0.5
for p in range(self.chunks[chunk_number]['dblocks']):
dset[...,p] /= self.cell_size[self.chunks[chunk_number]['levels'][p]]
@ -887,27 +893,26 @@ class Amun:
p = self.attributes['ndims'] - 3
tmp = self.__read_binary_data__('magz', chunk_number)
if p >= 0:
dset = np.roll(tmp, -1, axis=p) - np.roll(tmp, 1, axis=p)
dset = numpy.roll(tmp, -1, axis=p) - numpy.roll(tmp, 1, axis=p)
dset *= 0.5
for p in range(self.chunks[chunk_number]['dblocks']):
dset[...,p] /= self.cell_size[self.chunks[chunk_number]['levels'][p]]
else:
dset = np.zeros_like(tmp)
dset = numpy.zeros_like(tmp)
else:
dset = self.__read_binary_data__(dataset, chunk_number)
return dset
def dataset(self, dataset_name, extent=None, maxlev=-1, shrink=1, \
def dataset(self, dataset_name, extent=None, maxlev=None, shrink=1, \
interpolation='rebin', order=3, progress=False):
"""
Function returns dataset of the requested variable resampled to
the uniform mesh.
"""
from .interpolation import interpolate
import numpy as np
import sys
import numpy, sys
if self.dataformat == None:
raise Exception("Snapshot object has not been properly initialized!")
@ -915,25 +920,26 @@ class Amun:
if not dataset_name in self.variables:
raise Exception("Dataset '{}' is not available!\nAvailable datasets: {}\n".format(dataset_name, list(self.variables.keys())))
dlo = np.array([self.attributes['xmin'], self.attributes['ymin']])
dup = np.array([self.attributes['xmax'], self.attributes['ymax']])
dln = np.array([self.attributes['xlen'], self.attributes['ylen']])
dlo = numpy.array([self.attributes['xmin'], self.attributes['ymin']])
dup = numpy.array([self.attributes['xmax'], self.attributes['ymax']])
dln = numpy.array([self.attributes['xlen'], self.attributes['ylen']])
if self.attributes['ndims'] == 3:
dlo = np.append(dlo, self.attributes['zmin'])
dup = np.append(dup, self.attributes['zmax'])
dln = np.append(dln, self.attributes['zlen'])
dlo = numpy.append(dlo, self.attributes['zmin'])
dup = numpy.append(dup, self.attributes['zmax'])
dln = numpy.append(dln, self.attributes['zlen'])
slo = np.array(dlo)
sup = np.array(dup)
slo = numpy.array(dlo)
sup = numpy.array(dup)
if extent != None:
if len(extent) != 2 * self.attributes['ndims']:
raise Exception("Wrong dimensions of the argument 'extent'!")
slo = np.array(extent)[0::2]
sup = np.array(extent)[1::2]
slo = numpy.array(extent)[0::2]
sup = numpy.array(extent)[1::2]
if any(slo > dup) or any(sup < dlo) or any (slo >= sup):
raise Exception("Wrong order of the dimensions in the argument 'extent'!")
if isinstance(maxlev,int):
if maxlev != None:
if isinstance(maxlev, (int, numpy.int32)):
if 1 <= maxlev <= self.attributes['maxlev']:
shrink = 2**(self.attributes['maxlev']-maxlev)
elif maxlev > self.attributes['maxlev']:
@ -941,21 +947,21 @@ class Amun:
else:
raise Exception("Argument 'maxlev' must be an integer between 1 and {}.\n".format(self.attributes['maxlev']))
bm = np.array([ self.attributes['ncells'] ]*self.attributes['ndims'])
bm = numpy.array([ self.attributes['ncells'] ]*self.attributes['ndims'])
if self.attributes['ndims'] == 3:
rm = np.array([self.attributes['zblocks'], self.attributes['yblocks'], self.attributes['xblocks']])
rm = numpy.array([self.attributes['zblocks'], self.attributes['yblocks'], self.attributes['xblocks']])
else:
rm = np.array([self.attributes['yblocks'], self.attributes['xblocks']])
rm = numpy.array([self.attributes['yblocks'], self.attributes['xblocks']])
dm = rm * self.attributes['ncells'] * 2**(self.attributes['toplev'] - 1) // shrink
ll = np.array(np.floor(dm[::-1] * (slo - dlo) / dln), dtype='int')
uu = np.array( np.ceil(dm[::-1] * (sup - dlo) / dln), dtype='int')
ll = numpy.array(numpy.floor(dm[::-1] * (slo - dlo) / dln), dtype='int')
uu = numpy.array( numpy.ceil(dm[::-1] * (sup - dlo) / dln), dtype='int')
dm = (uu - ll)[::-1]
if dataset_name in [ 'velocity', 'vorticity', 'magnetic field', 'current density', 'electric field']:
arr = [ np.zeros(dm[:]), np.zeros(dm[:]), np.zeros(dm[:]) ]
arr = [ numpy.zeros(dm[:]), numpy.zeros(dm[:]), numpy.zeros(dm[:]) ]
else:
arr = np.zeros(dm[:])
arr = numpy.zeros(dm[:])
if progress:
sys.stdout.write("Snapshot's path:\n '{}'\n".format(self.path))
@ -980,8 +986,8 @@ class Amun:
if all(iu[:] > ll[:]) and all(il[:] < uu[:]):
nb = il[:] - ll[:]
ne = iu[:] - ll[:]
ib = np.maximum(nb[:], 0)
ie = np.minimum(ne[:], uu[:] - ll[:])
ib = numpy.maximum(nb[:], 0)
ie = numpy.minimum(ne[:], uu[:] - ll[:])
jb = ib[:] - nb[:]
je = ie[:] - ne[:] + cm[:]
@ -1015,10 +1021,9 @@ class Amun:
"""
Function converts dataset of the requested variable to AMR VTK file.
"""
import numpy, os, sys
from .octree import OcBase, OcNode
from .vtkio import WriteVTK
import numpy as np
import os, sys
if self.dataformat == None:
raise Exception("Snapshot object has not been properly initialized!")
@ -1077,12 +1082,12 @@ class Amun:
if progress:
sys.stdout.write("Generating OverlappingAMR VTK files\n")
bm = np.array([ self.attributes['ncells'] ]*self.attributes['ndims'])
bm = numpy.array([ self.attributes['ncells'] ]*self.attributes['ndims'])
if self.attributes['ndims'] == 3:
rm = np.array([self.attributes['zblocks'], self.attributes['yblocks'], self.attributes['xblocks']])
rm = numpy.array([self.attributes['zblocks'], self.attributes['yblocks'], self.attributes['xblocks']])
else:
rm = np.array([self.attributes['yblocks'], self.attributes['xblocks']])
rm = numpy.array([self.attributes['yblocks'], self.attributes['xblocks']])
ofile = "{}_{:06d}.vthb".format(dataset_name, self.attributes['isnap'])
opath = "{}_{:06d}".format(dataset_name, self.attributes['isnap'])
@ -1108,9 +1113,9 @@ class Amun:
no = 0
for item in base.getNodesFromLevel(lv):
lo = np.array(item.index) * bm
lo = numpy.array(item.index) * bm
up = lo + bm - 1
ll = np.stack((lo,up)).T.flatten()
ll = numpy.stack((lo,up)).T.flatten()
if item.hasData:
vfile = os.path.join(opath, fmt.format(dataset_name, lv, no))
WriteVTK(vfile, label, item.data, \

View File

@ -158,6 +158,7 @@ class AmunXML(Amun):
bounds = numpy.reshape(bounds, [2, 3, self.attributes['nleafs']])
for n in range(self.attributes['nchunks']):
if self.chunks[n]['dblocks'] > 0:
ids = self.__read_binary_data('ids', n, dtype='int32')
self.chunks[n]['levels'] = numpy.array([level[ids[p]] for p in range(self.chunks[n]['dblocks'])])
@ -166,6 +167,10 @@ class AmunXML(Amun):
ii = [ index[ids[p]] for p in range(self.chunks[n]['dblocks']) ]
self.chunks[n]['bounds'] = numpy.array([bounds[:,:,ii[p]].T for p in range(self.chunks[n]['dblocks'])]).T
else:
self.chunks[n]['levels'] = numpy.zeros((1), dtype = numpy.int32)
self.chunks[n]['coords'] = numpy.zeros((3,1), dtype = numpy.int32)
self.chunks[n]['bounds'] = numpy.zeros((2,3,1), dtype = numpy.float64)
def __read_binary_data__(self, dataset_name, chunk_number):

View File

@ -29,7 +29,6 @@
--------------------------------------------------------------------------------
"""
import numpy as np
try:
from scipy.ndimage import zoom
from scipy.interpolate import splrep, splev, interp1d, pchip_interpolate
@ -43,6 +42,8 @@ def rebin(a, newshape):
Subroutine changes the size of the input array to to new shape,
by copying cells or averaging them.
'''
import numpy
assert len(a.shape) == len(newshape)
m = a.ndim - 1
@ -58,7 +59,7 @@ def rebin(a, newshape):
return a.reshape(nn).mean(3).mean(1)
else:
for n in range(a.ndim):
a = np.repeat(a, newshape[n] // a.shape[n], axis=n)
a = numpy.repeat(a, newshape[n] // a.shape[n], axis=n)
return(a)
@ -66,6 +67,8 @@ def interpolate(a, newshape, nghosts=0, method=None, order=1):
'''
Subroutine rescales the block by interpolating its values.
'''
import numpy
if method == None or method == 'rebin' or not scipy_available:
ng = nghosts
@ -85,48 +88,48 @@ def interpolate(a, newshape, nghosts=0, method=None, order=1):
elif method in [ 'monotonic', 'pchip' ]:
dims = np.arange(a.ndim)
dims = numpy.arange(a.ndim)
q = a
for n in dims:
d2 = np.roll(q,-1, axis=0) + np.roll(q, 1, axis=0) - 2.0 * q
d2 = numpy.roll(q,-1, axis=0) + numpy.roll(q, 1, axis=0) - 2.0 * q
q = q - d2 / 24.0
d = np.array(q.shape)
d = numpy.array(q.shape)
xo = (np.arange(0.5, a.shape[n]) - nghosts) / (a.shape[n] - 2 * nghosts)
xn = np.arange(0.5, newshape[n]) / newshape[n]
xo = (numpy.arange(0.5, a.shape[n]) - nghosts) / (a.shape[n] - 2 * nghosts)
xn = numpy.arange(0.5, newshape[n]) / newshape[n]
u = q.reshape([d[0], q.size // d[0]])
f = np.zeros([newshape[n], q.size // d[0]])
f = numpy.zeros([newshape[n], q.size // d[0]])
for i in range(q.size // d[0]):
f[:,i] = pchip_interpolate(xo, u[:,i], xn)
d[0] = newshape[n]
f = f.reshape(d)
q = f.transpose(np.roll(dims, -1))
q = f.transpose(numpy.roll(dims, -1))
return q
elif method == 'spline':
dims = np.arange(a.ndim)
dims = numpy.arange(a.ndim)
q = a
for n in dims:
d2 = np.roll(q,-1, axis=0) + np.roll(q, 1, axis=0) - 2.0 * q
d2 = numpy.roll(q,-1, axis=0) + numpy.roll(q, 1, axis=0) - 2.0 * q
q = q - d2 / 24.0
d = np.array(q.shape)
d = numpy.array(q.shape)
xo = (np.arange(0.5, a.shape[n]) - nghosts) / (a.shape[n] - 2 * nghosts)
xn = np.arange(0.5, newshape[n]) / newshape[n]
xo = (numpy.arange(0.5, a.shape[n]) - nghosts) / (a.shape[n] - 2 * nghosts)
xn = numpy.arange(0.5, newshape[n]) / newshape[n]
u = q.reshape([d[0], q.size // d[0]])
f = np.zeros([newshape[n], q.size // d[0]])
f = numpy.zeros([newshape[n], q.size // d[0]])
for i in range(q.size // d[0]):
t = splrep(xo, u[:,i], k=5, s=0.0)
f[:,i] = splev(xn, t)
@ -134,27 +137,27 @@ def interpolate(a, newshape, nghosts=0, method=None, order=1):
d[0] = newshape[n]
f = f.reshape(d)
q = f.transpose(np.roll(dims, -1))
q = f.transpose(numpy.roll(dims, -1))
return q
else:
dims = np.arange(a.ndim)
dims = numpy.arange(a.ndim)
q = a
for n in dims:
d2 = np.roll(q,-1, axis=0) + np.roll(q, 1, axis=0) - 2.0 * q
d2 = numpy.roll(q,-1, axis=0) + numpy.roll(q, 1, axis=0) - 2.0 * q
q = q - d2 / 24.0
d = np.array(q.shape)
d = numpy.array(q.shape)
xo = (np.arange(0.5, a.shape[n]) - nghosts) / (a.shape[n] - 2 * nghosts)
xn = np.arange(0.5, newshape[n]) / newshape[n]
xo = (numpy.arange(0.5, a.shape[n]) - nghosts) / (a.shape[n] - 2 * nghosts)
xn = numpy.arange(0.5, newshape[n]) / newshape[n]
u = q.reshape([d[0], q.size // d[0]])
f = np.zeros([newshape[n], q.size // d[0]])
f = numpy.zeros([newshape[n], q.size // d[0]])
for i in range(q.size // d[0]):
t = interp1d(xo, u[:,i], kind=method)
f[:,i] = t(xn)
@ -162,6 +165,6 @@ def interpolate(a, newshape, nghosts=0, method=None, order=1):
d[0] = newshape[n]
f = f.reshape(d)
q = f.transpose(np.roll(dims, -1))
q = f.transpose(numpy.roll(dims, -1))
return q

View File

@ -176,7 +176,7 @@ class OcNode(object):
''' Function populates all nodes at lower levels with data from higher levels '''
def populateNodeData(self):
from .interpolation import rebin
import numpy as np
import numpy
if not self.isLeaf:
for n in range(len(self.children)):
@ -191,7 +191,7 @@ class OcNode(object):
bm = comp.shape
dm = [ 2 * d for d in bm ]
arr = np.zeros(dm, dtype=comp.dtype)
arr = numpy.zeros(dm, dtype=comp.dtype)
for k, j, i in itertools.product(range(2), range(2), range(2)):
n = (k * 2 + j) * 2 + i
@ -206,7 +206,7 @@ class OcNode(object):
bm = self.children[0].data.shape
dm = [ 2 * d for d in bm ]
arr = np.zeros(dm, dtype=self.children[0].data.dtype)
arr = numpy.zeros(dm, dtype=self.children[0].data.dtype)
for k, j, i in itertools.product(range(2), range(2), range(2)):
n = (k * 2 + j) * 2 + i