Last active
January 3, 2020 22:05
-
-
Save pshriwise/da30da3daf08594dddf1a58b9f10dcc8 to your computer and use it in GitHub Desktop.
Convert OpenMC Meshtally to MOAB/VTK mesh
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#!/usr/env/python3 | |
""" | |
Script for converting OpenMC mesh in an OpenMC | |
statepoint file to other mesh formats for visualization | |
""" | |
import argparse | |
import sys | |
import math | |
import openmc | |
import numpy as np | |
def write_moab(xs, ys, zs, tally_label, tally_data, error_data, outfile): | |
# attempt to import pymoab | |
try: | |
from pymoab import core | |
from pymoab.hcoord import HomCoord | |
from pymoab.scd import ScdInterface | |
from pymoab import types | |
except (ImportError, ModuleNotFoundError): | |
msg = "Conversion to MOAB .h5m file requested," \ | |
"but PyMOAB is not installed" | |
raise ImportError(msg) | |
mb = core.Core() | |
scd = ScdInterface(mb) | |
coords = [] | |
for k in zs: | |
for j in ys: | |
for i in xs: | |
coords += [i,j,k] | |
low = HomCoord([0, 0, 0, 0]) | |
high = HomCoord([len(xs) - 1, len(ys) - 1, len(zs) - 1, 0]) | |
scdbox = scd.construct_box(low, high, coords) | |
hexes = mb.get_entities_by_type(0, types.MBHEX) | |
tally_tag = mb.tag_get_handle(tally_label,1,types.MB_TYPE_DOUBLE,types.MB_TAG_DENSE,True) | |
error_tag = mb.tag_get_handle("error_tag",1,types.MB_TYPE_DOUBLE,types.MB_TAG_DENSE,True) | |
mb.tag_set_data(tally_tag,hexes,tally_data) | |
mb.tag_set_data(error_tag,hexes,error_data) | |
print('Writing %s' % outfile) | |
mb.write_file(outfile) | |
def write_vtk(xs, ys, zs, tally_label, tally_data, error_data, outfile): | |
try: | |
import vtk | |
except (ImportError, ModuleNotFoundError) as e: | |
msg = "Conversion to VTK requested," \ | |
"but the Python VTK module is not installed." | |
raise ImportError(msg) | |
vtk_box = vtk.vtkRectilinearGrid() | |
vtk_box.SetDimensions(len(xs), len(ys), len(zs)) | |
vtk_x_array = vtk.vtkDoubleArray() | |
vtk_x_array.SetName('x-coords') | |
vtk_x_array.SetArray(xs, len(xs), True) | |
print(vtk_x_array) | |
vtk_box.SetXCoordinates(vtk_x_array) | |
vtk_y_array = vtk.vtkDoubleArray() | |
vtk_y_array.SetName('y-coords') | |
vtk_y_array.SetArray(ys, len(ys), True) | |
vtk_box.SetYCoordinates(vtk_y_array) | |
vtk_z_array = vtk.vtkDoubleArray() | |
vtk_z_array.SetName('z-coords') | |
vtk_z_array.SetArray(zs, len(zs), True) | |
vtk_box.SetZCoordinates(vtk_z_array) | |
tally = np.array(tally_data) | |
tally_data = vtk.vtkDoubleArray() | |
tally_data.SetName(tally_label) | |
tally_data.SetArray(tally, tally.size, True) | |
error = np.array(error_data) | |
error_data = vtk.vtkDoubleArray() | |
error_data.SetName("error_tag") | |
error_data.SetArray(error, error.size, True) | |
vtk_box.GetCellData().AddArray(tally_data) | |
vtk_box.GetCellData().AddArray(error_data) | |
writer = vtk.vtkRectilinearGridWriter() | |
writer.SetFileName(outfile) | |
writer.SetInputData(vtk_box) | |
print('Writing %s' % outfile) | |
writer.Write() | |
def main(): | |
ap = argparse.ArgumentParser(description=__doc__) | |
ap.add_argument('-i','--input', | |
required=True, | |
help='Path to statepoint h5 file') | |
ap.add_argument('-n','--tally-name', | |
help='Tally name to add to mesh') | |
ap.add_argument('-t','--tally-id', | |
type=int, | |
required=True, | |
help='Tally name to add to mesh') | |
ap.add_argument('-m','--mesh-id', | |
type=int, | |
required=True, | |
help='ID of the mesh to convert') | |
ap.add_argument('-o', '--output', | |
action='store', | |
default='meshtally.vtk', | |
help='Name of outputfile (.h5m for MOAB, .vtk for VTK)') | |
args = ap.parse_args() | |
print('Loading file %s' % args.input) | |
sp = openmc.StatePoint(args.input) | |
print('Loading mesh with ID of %s' % args.mesh_id) | |
mesh = sp.meshes[args.mesh_id] | |
xs = np.linspace(mesh.lower_left[0], | |
mesh.upper_right[0], | |
mesh.dimension[0] + 1) | |
ys = np.linspace(mesh.lower_left[1], | |
mesh.upper_right[1], | |
mesh.dimension[1] + 1) | |
zs = np.linspace(mesh.lower_left[2], | |
mesh.upper_right[2], | |
mesh.dimension[2] + 1) | |
tally_args = {'name' : args.tally_name, | |
'id' : args.tally_id} | |
msg = "Loading retrieving tally with \n" | |
if args.tally_name is not None: | |
msg += "Name: {}\n".format(args.tally_name) | |
if args.tally_id is not None: | |
msg += "ID: {}\n".format(args.tally_id) | |
try: | |
tally = sp.get_tally(**tally_args) | |
except LookupError as e: | |
raise e | |
data = tally.mean[:,0,0] | |
error = tally.std_dev[:,0,0] | |
data = data.tolist() | |
error = error.tolist() | |
for c, i in enumerate(data): | |
if math.isnan(i): | |
data[c] = 0. | |
for c, i in enumerate(error): | |
if math.isnan(i): | |
error[c] = 0. | |
if args.tally_name is None: | |
tally_label = "tally_{}".format(args.tally_id) | |
else: | |
tally_label = args.tally_name | |
if args.output.endswith(".vtk"): | |
write_vtk(xs, ys, zs, tally_label, data, error, args.output) | |
else: | |
write_moab(xs, ys, zs, tally_label, data, error, args.output) | |
if __name__ == "__main__": | |
main() |
Hi Jon,
Sorry this email slipped through the cracks! Yeah we could certainly add
that as a flag. We could also look into some other options for limiting how
many of the HDF5 datasets are in memory at once.
I'll set a reminder to take a look some time later this week. Sound good?
Cheers,
Patrick Shriwise
608-446-8173
…On Thu, Dec 12, 2019 at 11:13 AM Jonathan Shimwell ***@***.***> wrote:
I have been stress testing this for large meshes, in the case where the
mesh tally is large then making both the tally mesh and the error mesh
pushes the computer over its RAM limits and I get a crash. I was wondering
if we can have a flag that allows the script include either the tally or
the tally and the error or just the error. This will allow us to make the
largest VTK files within the RAM limits
—
You are receiving this because you authored the thread.
Reply to this email directly, view it on GitHub
<https://gist.github.com/da30da3daf08594dddf1a58b9f10dcc8?email_source=notifications&email_token=ABC2HZMWYKSNEDYP57ODYFDQYJWL3A5CNFSM4IXNVDO2YY3PNVWWK3TUL52HS4DFVNDWS43UINXW23LFNZ2KUY3PNVWWK3TUL5UWJTQAF5ZVY#gistcomment-3109724>,
or unsubscribe
<https://github.com/notifications/unsubscribe-auth/ABC2HZK5GXBYLHQ3GZ3ZVLLQYJWL3ANCNFSM4IXNVDOQ>
.
@shimwell, were you seeing an issue with memory limits when writing a VTK or MOAB file? Maybe both?
I was using the VTK option, I removed the error tally and writing of the tally and this resulted in lower memory usage.
So I was thinking there could be some user option to allow more control of writing tally , error or both (default)
Another argument for the ap, perhaps something like this
ap.add_argument('-e', '--error',
choices=['error', 'tally'],
nargs='+',
default=['error', 'tally'],
help="Error and tally data included by default. Use this option to write just the error or just the tally data")
But this is not so important and can always be added later
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
I have been stress testing this for large meshes, in the case where the mesh tally is large then making both the tally mesh and the error mesh pushes the computer over its RAM limits and I get a crash. I was wondering if we can have a flag that allows the script include either the tally or the tally and the error or just the error. This will allow us to make the largest VTK files within the RAM limits