Ошибка значения при преобразовании файлов с помощью сценария оболочки ⇐ Python

[Anonymous]

Я запускаю файл sh для преобразования разрозненных файлов в необработанные файлы. На входе указывается путь к файлам хранилища (включая имя моделирования), начальный индекс и конечный индекс интересующих файлов. Для начала решил попробовать индексы 1 и 2. Но выходит ошибка. Ниже приведен файл sh, текст ошибки, raw_reader.py и Compute_silo_derivative.py
#!/bin/bash

set -e

if (($# != 3)); then
echo "usage: $0 path/simulation_name start_index end_index"
exit 1
fi

afivo_dir=/home/erg/afivo-streamer-master/afivo

vars=('rhs' 'Je_1' 'Je_2' 'Je_3')

for ((i="$2"; i 1e-10)

if overlap:
grid_data, ix_lo, ix_hi = map_grid_data_to(
ga, gc['r_min'], gc['r_max'], gc['dr'],
args.axisymmetric, args.interpolation)

# Account for indexing offset in 'values'
g_ix = tuple([np.s_[i:j] for (i, j) in
zip(ix_lo+gc['ilo'], ix_hi+gc['ilo'])])
gc['values'][g_ix] = grid_data

return grids_c

if __name__ == '__main__':
args = p.parse_args()

all_gd = [load_file(f, args.project_dims, args.variable)
for f in args.input_files]

grids = [gd[0] for gd in all_gd]
domains = [gd[1] for gd in all_gd]
times = np.array([d['time'] for d in domains])
cycles = np.array([d['cycle'] for d in domains])

if args.deriv_type == '1st_central':
if len(grids) != 2:
raise ValueError('1st_central requires len(grids) == 2')

# Map onto last grid
grids[0] = map_grid_a_to_b(grids[0], grids[1],
args.axisymmetric, args.interpolation)

# Compute 1st order derivative
dt = times[1] - times[0]
w = np.array([-1, 1]) / dt

# Store results in grids[0]
for i in range(len(grids[0])):
grids[0]['values'] = w[0] * grids[0]['values'] + \
w[1] * grids[1]['values']

domains[0]['time'] = 0.5 * (times[0] + times[1])
domains[0]['cycle'] = int((cycles[0] + cycles[1])//2)

elif args.deriv_type == '2nd_central':
if len(grids) != 3:
raise ValueError('2nd_central requires len(grids) == 3')

# Map onto last grid
grids[0] = map_grid_a_to_b(grids[0], grids[2],
args.axisymmetric, args.interpolation)
grids[1] = map_grid_a_to_b(grids[1], grids[2],
args.axisymmetric, args.interpolation)
t_grid = times[1]

# Compute 2nd order derivative, allowing for variable dt
dt_a = times[1] - times[0]
dt_b = times[2] - times[1]
w = np.zeros(3)

w[0] = dt_b/dt_a * 2 / (dt_b**2 + dt_a * dt_b)
w[2] = 2 / (dt_b**2 + dt_a * dt_b)
w[1] = -(w[0] + w[2])

# Store results in grids[0]
for i in range(len(grids[0])):
grids[0]['values'] = w[0] * grids[0]['values'] + \
w[1] * grids[1]['values'] + \
w[2] * grids[2]['values']

domains[0]['time'] = times[1]
domains[0]['cycle'] = cycles[1]

write_raw_data(args.output, grids[0], domains[0])
print(f'Written {args.output}')

#!/usr/bin/env python3

# Module to read and convert raw data obtained by transforming Silo files.
#
# Author: Jannis Teunissen (2022)

import numpy as np
import copy
from struct import pack, unpack, calcsize
from scipy.interpolate import RegularGridInterpolator
import tempfile
import os
import subprocess
import pathlib

def load_file(fname, project_dims=None, axisymmetric=False,
variable=None, silo_to_raw=None):
"""Read a Silo or a raw file

:param fname: name of the raw file
:param project_dims: project along these dimensions
:param axisymmetric: bool, whether the data is axisymmetric
:param variable: read this variable from a Silo file
:param silo_to_raw: path to silo_to_raw executable
:returns: grids, domains
"""
extension = os.path.splitext(fname)[1]

if extension == ".silo":
if not variable:
raise ValueError('Specify variable to plot for a Silo file')

if silo_to_raw is None:
this_folder = pathlib.Path(__file__).parent.resolve()
silo_to_raw = os.path.join(this_folder, 'silo_to_raw')

if not os.path.exists(silo_to_raw):
raise ValueError(f'Could not find silo_to_raw at {silo_to_raw}')

# Convert silo to raw. Place temporary files in the same folder as
# there can otherwise be issues with a full /tmp folder
dirname = os.path.dirname(fname)
with tempfile.NamedTemporaryFile(dir=dirname) as fp:
_ = subprocess.call([silo_to_raw, fname,
variable, fp.name])
grids, domain = get_raw_data(fp.name, project_dims, axisymmetric)
else:
grids, domain = get_raw_data(fname, project_dims)

return grids, domain

def read_single_grid(f):
"""Read single grid from binary file"""
n_dims = unpack('=i', f.read(calcsize('=i')))[0]

fmt = '=' + str(n_dims) + 'i'
dims = np.array(unpack(fmt, f.read(calcsize(fmt))))

# lo and hi index range of non-phony data
ilo = np.array(unpack(fmt, f.read(calcsize(fmt))))
ihi = np.array(unpack(fmt, f.read(calcsize(fmt))))

# Mesh coordinates
coords = []
coords_cc = []
for d in dims:
fmt = '=' + str(d) + 'd'
tmp = np.array(unpack(fmt, f.read(calcsize(fmt))))
coords.append(tmp)
# Cell-centered coordinates
coords_cc.append(0.5 * (tmp[1:] + tmp[:-1]))

# Number of cell centers is one less than number of faces
n_cells = dims-1
fmt = '=' + str(np.prod(n_cells)) + 'd'
tmp = unpack(fmt, f.read(calcsize(fmt)))
vals = np.array(tmp).reshape(n_cells, order='F')

grid = {
'n_dims': n_dims,
'dims': dims,
'ilo': ilo,
'ihi': ihi,
'coords': coords,
'coords_cc': coords_cc,
'values': vals
}

return grid

def get_raw_data(fname, project_dims=None, axisymmetric=False):
"""Read raw data extracted from a Silo file

:param fname: filename of raw data
:param project_dims: int, integrate over these dimensions
:param axisymmetric: bool, whether the data is axisymmetric
:returns: grids and domain properties

"""
with open(fname, 'rb') as f:
cycle = unpack('=i', f.read(calcsize('=i')))[0]
time = unpack('=d', f.read(calcsize('=d')))[0]

grids = []
n_grids = unpack('=i', f.read(calcsize('=i')))[0]

for i in range(n_grids):
grid = read_single_grid(f)

# Add properties
props = get_grid_properties(grid)
grid.update(props)

if project_dims is not None:
grid = grid_project(grid, project_dims, axisymmetric)

grids.append(grid)

domain_props = get_domain_properties(grids)
domain_props['cycle'] = cycle
domain_props['time'] = time
return grids, domain_props

def write_single_grid(f, grid):
"""Write single grid to binary file"""
f.write(pack('=i', grid['n_dims']))

fmt = '=' + str(grid['n_dims']) + 'i'
f.write(pack(fmt, *grid['dims']))

# lo and hi index range of non-phony data
f.write(pack(fmt, *grid['ilo']))
f.write(pack(fmt, *grid['ihi']))

# Mesh coordinates
for i in range(grid['n_dims']):
fmt = '=' + str(grid['dims']) + 'd'
f.write(pack(fmt, *grid['coords']))

# Number of cell centers is one less than number of faces
n_cells = grid['dims']-1
fmt = '=' + str(np.prod(n_cells)) + 'd'
f.write(grid['values'].tobytes(order='F'))

def write_raw_data(fname, grids, domain):
"""Write grid data to a raw file again

:param fname: filename of raw data
:param grids: list of grids
:param domain: domain properties
"""
with open(fname, 'wb') as f:
f.write(pack('=i', domain['cycle']))
f.write(pack('=d', domain['time']))
f.write(pack('=i', len(grids)))

for g in grids:
write_single_grid(f, g)

def get_grid_properties(g):
"""Return some additional properties for a grid"""
# Note: ghost cells are excluded for r_min and r_max
r_min = [c for c, i in zip(g['coords'], g['ilo'])]
r_max = [c[i] for c, i in zip(g['coords'], g['ihi'])]
dr = [c[1] - c[0] for c in g['coords']]

props = {
'r_min': np.array(r_min),
'r_max': np.array(r_max),
'dr': np.array(dr)
}
return props

def get_domain_properties(grids):
"""Return properties for a domain"""
n_dims = grids[0]['n_dims']
r_min = np.full(n_dims, 1e100)
r_max = np.full(n_dims, -1e100)
dr_max = np.zeros(n_dims)
dr_min = np.full(n_dims, 1e100)

for g in grids:
r_min = np.minimum(r_min, g['r_min'])
r_max = np.maximum(r_max, g['r_max'])
dr_min = np.minimum(dr_min, g['dr'])
dr_max = np.maximum(dr_max, g['dr'])

nx_coarse = np.rint((r_max - r_min) / dr_max)

props = {
'n_dims': n_dims,
'r_min': np.array(r_min),
'r_max': np.array(r_max),
'dr_min': np.array(dr_min),
'dr_max': np.array(dr_max),
'n_cells_coarse': nx_coarse.astype(int)
}

return props

def grid_project(in_grid, project_dims, axisymmetric):
"""Integrate grid data along one or more dimensions"""
g = copy.deepcopy(in_grid)

pdims = np.sort(project_dims)

# Index range to use for values below. Along the projected dimensions,
# exclude ghost cells.
ilo = 0 * g['dims']
ihi = 1 * g['dims'] - 1 # number of faces - 1
ilo[pdims] = g['ilo'][pdims]
ihi[pdims] = g['ihi'][pdims]

# Get index range corresponding to non-ghost grid cells
valid_ix = tuple([np.s_[i:j] for (i, j) in zip(ilo, ihi)])

# Sum values along projected dimensions
if axisymmetric and 0 in pdims:
# Multiply with 2 * pi * r
r = g['coords_cc'][0][valid_ix[0]]
w = np.prod(g['dr'][pdims]) * 2 * np.pi * r

# Broadcast to have volume weight for every grid cell
w = np.broadcast_to(w[:, None], ihi-ilo)
g['values'] = (g['values'][valid_ix] * w).sum(axis=tuple(pdims))
else:
fac = np.prod(g['dr'][pdims])
g['values'] = fac * g['values'][valid_ix].sum(axis=tuple(pdims))

g['n_dims'] -= len(project_dims)
g['dims'] = np.delete(g['dims'], pdims)
g['ilo'] = np.delete(g['ilo'], pdims)
g['ihi'] = np.delete(g['ihi'], pdims)
g['r_min'] = np.delete(g['r_min'], pdims)
g['r_max'] = np.delete(g['r_max'], pdims)
g['dr'] = np.delete(g['dr'], pdims)

for d in pdims[::-1]: # Reverse order
del g['coords'][d]
del g['coords_cc'][d]

return g

def map_grid_data_to(g, r_min, r_max, dr, axisymmetric=False,
interpolation_method='linear'):
"""Map grid data to another grid with origin r_min and grid spacing dr

:param g: input grid
:param r_min: origin of new grid
:param r_max: upper location of new grid
:param dr: grid spacing of new grid
:param axisymmetric: whether to account for an axisymmetric geometry
:param interpolation_method: how to interpolate data (linear, nearest)
:returns: data and indices on new grid

"""
# Check if grids overlap
if np.any(g['r_min'] > r_max) or np.any(g['r_max'] < r_min):
Ndim = len(r_min)
# Return empty index range
return 0., np.zeros(Ndim, dtype=int), np.zeros(Ndim, dtype=int)

ratios = dr / g['dr']

if not np.allclose(ratios, ratios[0], atol=0.):
raise ValueError('Grid ratio not uniform')

# Get index range corresponding to non-ghost grid cells
valid_ix = tuple([np.s_[i:j] for (i, j) in zip(g['ilo'], g['ihi'])])

# To avoid issues due to numerical round-off errors
eps = 1e-10

# Compute coarse grid min and max index
if ratios[0] > 1 + eps:
# Fine-to-coarse, first compute index for fine grid
ix_lo_fine = np.round((g['r_min'] - r_min)/g['dr']).astype(int)
ix_hi_fine = np.round((g['r_max'] - r_min)/g['dr']).astype(int) - 1

# Then convert to coarse grid index
r = np.round(ratios).astype(int)
ix_lo, ix_hi = ix_lo_fine//r, ix_hi_fine//r
else:
# Grid is at same refinement level or coarser, so we can directly
# compute index with rounding
ix_lo = np.round((g['r_min'] - r_min)/dr).astype(int)
ix_hi = np.round((g['r_max'] - r_min)/dr).astype(int) - 1

nx = ix_hi - ix_lo + 1

# Get index range that is valid on uniform grid
nx_uniform = np.round((r_max - r_min)/dr).astype(int)
offset_lo = np.maximum(0, -ix_lo)
offset_hi = np.maximum(0, ix_hi+1 - nx_uniform)

# Index range on the grid_data
grid_lo, grid_hi = offset_lo, nx-offset_hi
d_ix = tuple([np.s_[i:j] for (i, j) in zip(grid_lo, grid_hi)])

if ratios[0] > 1 + eps:
# Reduce resolution. Determine coarse indices for each cell center
cix = []
coords_fine = []
coords_coarse = []
for d in range(g['n_dims']):
dim_coords = g['coords_cc'][d]
cc_fine = dim_coords[g['ilo'][d]:g['ihi'][d]]
ix = np.floor((cc_fine - r_min[d])/dr[d]).astype(int)
cc_coarse = r_min[d] + (ix + 0.5) * dr[d]

cix.append(ix - ix_lo[d])
coords_fine.append(cc_fine)
coords_coarse.append(cc_coarse)

# Create meshgrid of coarse indices
ixs = np.meshgrid(*cix, indexing='ij')

# Add fine grid values at coarse indices, weighted by relative volume
cdata = np.zeros(nx)

if axisymmetric:
rvolume = np.prod(g['dr']/dr) * coords_fine[0]/coords_coarse[0]

if rvolume.ndim < len(g['ihi']):
# Broadcast to have volume weight for every grid cell
rvolume = np.broadcast_to(rvolume[:, None], g['ihi']-g['ilo'])

# Important to use Fortran order here
values = rvolume.ravel() * g['values'][valid_ix].ravel()
else:
# Cartesian grid, simple averaging
rvolume = np.prod(g['dr']/dr)
values = rvolume * g['values'][valid_ix].ravel()

np.add.at(cdata, tuple(map(np.ravel, ixs)), values)

# Extract region that overlaps with uniform grid
cdata = cdata[d_ix]
elif ratios[0] < 1 - eps:
# To interpolate data, compute coordinates of new cell centers
# TODO: could maybe include axisymmetric correction here as well
r0 = g['r_min'] + (offset_lo + 0.5) * dr
r1 = g['r_max'] - (offset_hi + 0.5) * dr

c_new = [np.linspace(a, b, n) for a, b, n in
zip(r0, r1, grid_hi - grid_lo)]
mgrid = np.meshgrid(*c_new, indexing='ij')
new_coords = np.vstack(tuple(map(np.ravel, mgrid))).T

# Near physical boundaries, extrapolation will be performed when using
# linear interpolation
f_interp = RegularGridInterpolator(
tuple(g['coords_cc']), g['values'], bounds_error=False,
fill_value=None, method=interpolation_method)

cdata = f_interp(new_coords).reshape(grid_hi - grid_lo)
else:
# Can directly use available data
cdata = g['values'][valid_ix]

# Extract region that overlaps with uniform grid
cdata = cdata[d_ix]

return cdata, ix_lo+offset_lo, ix_hi+1-offset_hi


Подробнее здесь: https://stackoverflow.com/questions/798 ... ell-script