|
#!/usr/bin/env python |
|
# Copyright (c) 2018 Parsa Amini |
|
# Copyright (c) 2018 Patrick Diehl |
|
# |
|
# Distributed under the Boost Software License, Version 1.0. (See accompanying |
|
# file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) |
|
|
|
# ## Synopsis |
|
# ``` |
|
# usage: octoviz_hdf5.py [-h] hdf5_file vtu_file |
|
# |
|
# Generate VTU from an Octo-tiger Silo file |
|
# |
|
# positional arguments: |
|
# hdf5_file Octo-tiger .silo file to read |
|
# vtu_file Path to write the .vtu file |
|
# |
|
# optional arguments: |
|
# -h, --help show this help message and exit |
|
# ``` |
|
|
|
from __future__ import print_function |
|
import argparse |
|
import os |
|
import sys |
|
|
|
|
|
def silo_to_vtu(hdf5_file, vtu_file): |
|
print('Importing h5py, numpy, pandas, and vtk...') |
|
try: |
|
import h5py |
|
import numpy as np |
|
import pandas as pd |
|
import vtk |
|
except ImportError: |
|
print( |
|
'h5py, numpy, pandas, and vtk modules are required to run this script.', |
|
file=sys.stderr |
|
) |
|
sys.exit(1) |
|
|
|
db = h5py.File(hdf5_file, 'r') |
|
print('File {} opened.'.format(hdf5_file)) |
|
|
|
# ## 1. Field Variables |
|
# `n_species` Tells me how many $\rho$ fractions there are |
|
rho_fractions = db['n_species'].value[0] |
|
print('n_species: {}'.format(rho_fractions)) |
|
rhos = ['rho_' + str(x) for x in range(1, rho_fractions + 1)] |
|
|
|
groups = [] |
|
datasets = [] |
|
for k, v in db.iteritems(): |
|
if isinstance(v, h5py.Group): |
|
groups.append(v) |
|
elif isinstance(v, h5py.Dataset): |
|
datasets.append(v) |
|
elif isinstance(v, h5py.Datatype): |
|
pass |
|
else: |
|
raise Exception('Bad key: ' + k) |
|
print('Number of directories: {}'.format((len(groups) - 2))) |
|
|
|
# Symbol | SILO variable name | Description |
|
# --- | --- | --- |
|
# $\rho$ | `rho_i` | $i$-th density fraction |
|
# $s_i$ | `sx, sy, sz` | inertial frame momentum |
|
# $z_i$ | `zx, zy, zz` | inertial fram spin momentum relative to cell center |
|
# $g_i$ | `gx, gy, gz` | gravitational acceleration |
|
# $\phi$ | `phi` | gravitational potential |
|
# $E$ | `egas` | kinetic + internal energy density |
|
# $\tau$ | `tau` | entropy tracer |
|
print('Reading directories...') |
|
col = [] |
|
|
|
for g in groups[2:]: |
|
gf = pd.DataFrame({ |
|
#'group': g.name, |
|
'E': db[g['egas'].attrs['silo'][1]].value.flatten(), |
|
'g_x': db[g['gx'].attrs['silo'][1]].value.flatten(), |
|
'g_y': db[g['gy'].attrs['silo'][1]].value.flatten(), |
|
'g_z': db[g['gz'].attrs['silo'][1]].value.flatten(), |
|
'phi': db[g['phi'].attrs['silo'][1]].value.flatten(), |
|
's_x': db[g['sx'].attrs['silo'][1]].value.flatten(), |
|
's_y': db[g['sy'].attrs['silo'][1]].value.flatten(), |
|
's_z': db[g['sz'].attrs['silo'][1]].value.flatten(), |
|
'tau': db[g['tau'].attrs['silo'][1]].value.flatten(), |
|
'z_x': db[g['zx'].attrs['silo'][1]].value.flatten(), |
|
'z_y': db[g['zy'].attrs['silo'][1]].value.flatten(), |
|
'z_z': db[g['zz'].attrs['silo'][1]].value.flatten(), |
|
}) |
|
for i in range(rho_fractions): |
|
gf[rhos[i]] = db[g[rhos[i]].attrs['silo'][1]].value.flatten() |
|
|
|
gm = np.vstack(( |
|
db[g['quadmesh'].attrs['silo'][0]].value, |
|
db[g['quadmesh'].attrs['silo'][1]].value, |
|
db[g['quadmesh'].attrs['silo'][2]].value, |
|
)) |
|
gp = .5 * (gm[:, :-1] + gm[:, 1:]) |
|
|
|
gf['z'] = gp[2].repeat(gp.shape[1] ** 2) |
|
gf['y'] = (np.tile(gp[1], gp.shape[1])).repeat(gp.shape[1]) |
|
gf['x'] = np.tile(gp[0], gp.shape[1] ** 2) |
|
|
|
col.append(gf) |
|
|
|
df = pd.concat(col).reset_index(drop=True) |
|
|
|
# ## 2. Other variables |
|
|
|
# Symbol | SILO variable name | Description |
|
# --- | --- | --- |
|
# $\mu_i$ | `A[i-1]` | atomic mass number of $i$-th fraction |
|
# $Z_i$ | `Z[i-1]` | electron number of $i$-th fraction |
|
# $\omega$ | `omega` | rotation frequency of the grid |
|
atomic_mass = np.array(db['atomic_mass'].value) |
|
atomic_number = np.array(db['atomic_number'].value) |
|
omega = db['omega'].value[0] |
|
|
|
CYCLE = db['cycle'].value |
|
TIME = db['dtime'].value |
|
|
|
# ## 3. Derived Expressions |
|
|
|
# total mass density |
|
# \rho := \Sigma_{i=1}^N \rho_i |
|
print('Calculating mass density...') |
|
df['rho'] = 0 |
|
|
|
for rho in rhos: |
|
df.rho += df[rho] |
|
|
|
# x - velocity relative to grid |
|
# v_x := \frac{s_x}{\rho} + y \Omega |
|
# v_x := \frac{s_{x_{ijk}}}{\rho_{ijk}} + \frac{1}{2}(y_j + y_{j+1}) \Omega |
|
print('Calculating v_x velocity relative to grid...') |
|
df['v_x'] = (df.s_x / df.rho) + df.y * omega |
|
|
|
# y - velocity relative to grid |
|
# v_y := \frac{s_y}{\rho} - x \Omega |
|
# v_y := \frac{s_{y_{ijk}}}{\rho_{ijk}} - \frac{1}{2}(x_i + x_{i+1}) \Omega |
|
print('Calculating v_y velocity relative to grid...') |
|
df['v_y'] = (df.s_y / df.rho) + df.x * omega |
|
|
|
# z - velocity relative to grid |
|
# v_z := \frac{s_z}{\rho} |
|
print('Calculating v_z velocity relative to grid...') |
|
df['v_z'] = df.s_z / df.rho |
|
|
|
# \gamma = \frac{5}{3} |
|
gamma = 5. / 3. |
|
|
|
# internal gas energy density |
|
# e := \begin{cases} |
|
# E - \frac{1}{2} \frac{s^2}{\rho} & \text{if } |
|
# E - \frac{1}{2} \frac{s^2}{\rho} > 0.001 E \\ |
|
# \tau^\gamma & \text{else} \\ |
|
# \end{cases} |
|
print('Calculating internal gas energy density...') |
|
df['e'] = df.E - .5 * ((df.s_x ** 2 + df.s_y **2 + df.s_z ** 2) / df.rho) |
|
df.e = np.where(df.e > .001 * df.E, df.e, df.tau ** gamma) |
|
|
|
# $m_H$ (mass of hydrogen) = $1.6733 \times 10^{24}$ |
|
mH = 1.6733e-24 |
|
|
|
# number density of ith fraction (ions + electrons) |
|
# n_i := \frac{\rho_i}{\mu_i m_H} \left ( 1 + Z_i \right ) |
|
print('Calculating number density of ith fraction...') |
|
for i in range(len(atomic_mass)): |
|
df['n_i' + str(i + 1)] = (df.rho / (atomic_mass[i] * mH)) * (1. + atomic_number[i]) |
|
|
|
# total number density |
|
# n := \Sigma_{i=1}^N n_i |
|
print('Calculating total number density...') |
|
df['n'] = 0 |
|
for i in range(len(atomic_mass)): |
|
df.n += df['n_i' + str(i + 1)] |
|
|
|
# $K_b := 1.380658 \times 10^{-16}$ |
|
K_b = 1.380658e-16 |
|
|
|
# temperature |
|
# T := \frac{1}{\gamma-1} \frac{e}{n} |
|
df['T'] = ((1. / (gamma - 1.)) * (df.e / df.n)) / K_b |
|
|
|
# pressure |
|
# P := \left(\gamma-1\right) e |
|
print('Calculating pressure...') |
|
df['P'] = (gamma - 1.) * df.e |
|
|
|
db.close() |
|
|
|
print('Number of particles: {}'.format(df.shape[0])) |
|
|
|
# 3. Write to VTU file |
|
|
|
# Add points and fields |
|
points = vtk.vtkPoints() |
|
points.SetNumberOfPoints(df.shape[0]) |
|
|
|
points.SetDataTypeToDouble() |
|
|
|
for idx, x_i, y_i, z_i in zip(range(df.shape[0]), df.x, df.y, df.z): |
|
points.InsertPoint(idx, x_i, y_i, z_i) |
|
|
|
# Configure the grid |
|
grid = vtk.vtkUnstructuredGrid() |
|
grid.SetPoints(points) |
|
|
|
field_data = grid.GetFieldData() |
|
|
|
atomic_mass_array = vtk.vtkDoubleArray() |
|
atomic_mass_array.SetName('atomic_mass') |
|
atomic_mass_array.SetNumberOfComponents(len(atomic_mass)) |
|
atomic_mass_array.SetNumberOfTuples(1) |
|
atomic_mass_array.SetTuple(0, atomic_mass) |
|
|
|
field_data.AddArray(atomic_mass_array) |
|
|
|
atomic_number_array = vtk.vtkDoubleArray() |
|
atomic_number_array.SetName('atomic_number') |
|
atomic_number_array.SetNumberOfComponents(len(atomic_number)) |
|
atomic_number_array.SetNumberOfTuples(1) |
|
atomic_number_array.SetTuple(0, atomic_number) |
|
|
|
field_data.AddArray(atomic_number_array) |
|
|
|
omega_array = vtk.vtkDoubleArray() |
|
omega_array.SetName('omega') |
|
omega_array.SetNumberOfComponents(1) |
|
omega_array.SetNumberOfTuples(1) |
|
omega_array.SetTuple1(0, omega) |
|
|
|
field_data.AddArray(omega_array) |
|
|
|
TIME_array = vtk.vtkDoubleArray() |
|
TIME_array.SetName('TIME') |
|
TIME_array.SetNumberOfComponents(1) |
|
TIME_array.SetNumberOfTuples(1) |
|
TIME_array.SetTuple1(0, TIME) |
|
|
|
field_data.AddArray(TIME_array) |
|
|
|
CYCLE_array = vtk.vtkDoubleArray() |
|
CYCLE_array.SetName('CYCLE') |
|
CYCLE_array.SetNumberOfComponents(1) |
|
CYCLE_array.SetNumberOfTuples(1) |
|
CYCLE_array.SetTuple1(0, CYCLE) |
|
|
|
field_data.AddArray(CYCLE_array) |
|
|
|
data_out = grid.GetPointData() |
|
|
|
keys = list(set(df.keys()).difference({'x', 'y', 'z', 'group'})) |
|
keys.sort() |
|
|
|
for k in keys: |
|
print('Adding key \'{}\'...'.format(k)) |
|
dis_array = vtk.vtkDoubleArray() |
|
dis_array.SetName(k) |
|
dis_array.SetNumberOfComponents(1) |
|
dis_array.SetNumberOfTuples(df.shape[0]) |
|
|
|
for idx, val in zip(range(df.shape[0]), df[k]): |
|
dis_array.SetTuple1(idx, val) |
|
|
|
data_out.AddArray(dis_array) |
|
|
|
# Configure the VTU writer |
|
writer = vtk.vtkXMLUnstructuredGridWriter() |
|
writer.SetFileName(vtu_file) |
|
|
|
writer.SetInputData(grid) |
|
|
|
# Writer Data modes: |
|
# 1. Binary (with compression), saves space |
|
writer.SetDataModeToBinary() |
|
# 2. ASCII - Plain text, human readable |
|
# writer.SetDataModeToAscii() |
|
|
|
print('Writing VTU file {}...'.format(vtu_file)) |
|
writer.Write() |
|
|
|
|
|
# Determine if a path is actually a readable file |
|
def existing_file(path): |
|
if not os.path.isfile(path) and os.access(path, os.R_OK): |
|
raise argparse.ArgumentTypeError('Path {} does not exist'.format(path)) |
|
|
|
return path |
|
|
|
|
|
def main(): |
|
parser = argparse.ArgumentParser( |
|
description='Generate VTU from an Octo-tiger Silo file') |
|
parser.add_argument( |
|
'hdf5_file', type=existing_file, help='Octo-tiger .silo file to read') |
|
parser.add_argument( |
|
'vtu_file', type=str, help='Path to write the .vtu file') |
|
|
|
args = parser.parse_args() |
|
|
|
silo_to_vtu(args.hdf5_file, args.vtu_file) |
|
|
|
|
|
if __name__ == '__main__': |
|
main() |
