SHARE
TWEET

Untitled

a guest Sep 23rd, 2019 71 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. #!/usr/env/python3
  2.  
  3. """
  4. Script for converting OpenMC mesh in an OpenMC
  5. statepoint file to other mesh formats for visualization
  6. """
  7.  
  8. import argparse
  9. import sys
  10. import math
  11. import openmc
  12. import numpy as np
  13.  
  14. def write_moab(xs, ys, zs, tally_label, tally_data, error_data, outfile):
  15.     # attempt to import pymoab
  16.     try:
  17.         from pymoab import core
  18.         from pymoab.hcoord import HomCoord
  19.         from pymoab.scd import ScdInterface
  20.         from pymoab import types
  21.     except (ImportError, ModuleNotFoundError):
  22.         msg = "Conversion to MOAB .h5m file requested," \
  23.               "but PyMOAB is not installed"
  24.         raise ImportError(msg)
  25.  
  26.     mb = core.Core()
  27.  
  28.     scd = ScdInterface(mb)
  29.  
  30.     coords = []
  31.     for k in zs:
  32.         for j in ys:
  33.             for i in xs:
  34.                 coords += [i,j,k]
  35.  
  36.     low = HomCoord([0, 0, 0, 0])
  37.     high = HomCoord([len(xs) - 1, len(ys) - 1, len(zs) - 1, 0])
  38.  
  39.     scdbox = scd.construct_box(low, high, coords)
  40.  
  41.     hexes = mb.get_entities_by_type(0, types.MBHEX)
  42.  
  43.     tally_tag = mb.tag_get_handle(tally_label,1,types.MB_TYPE_DOUBLE,types.MB_TAG_DENSE,True)
  44.     error_tag = mb.tag_get_handle("error_tag",1,types.MB_TYPE_DOUBLE,types.MB_TAG_DENSE,True)
  45.  
  46.     mb.tag_set_data(tally_tag,hexes,tally_data)
  47.     mb.tag_set_data(error_tag,hexes,error_data)
  48.  
  49.     print('Writing %s' % outfile)
  50.  
  51.     mb.write_file(outfile)
  52.  
  53. def write_vtk(xs, ys, zs, tally_label, tally_data, error_data, outfile):
  54.     try:
  55.         import vtk
  56.     except (ImportError, ModuleNotFoundError) as e:
  57.         msg = "Conversion to VTK requested," \
  58.               "but the Python VTK module is not installed."
  59.         raise ImportError(msg)
  60.  
  61.     vtk_box = vtk.vtkRectilinearGrid()
  62.  
  63.     vtk_box.SetDimensions(len(xs), len(ys), len(zs))
  64.  
  65.     vtk_x_array = vtk.vtkDoubleArray()
  66.     vtk_x_array.SetName('x-coords')
  67.     vtk_x_array.SetArray(xs, len(xs), True)
  68.     print(vtk_x_array)
  69.     vtk_box.SetXCoordinates(vtk_x_array)
  70.  
  71.     vtk_y_array = vtk.vtkDoubleArray()
  72.     vtk_y_array.SetName('y-coords')
  73.     vtk_y_array.SetArray(ys, len(ys), True)
  74.     vtk_box.SetYCoordinates(vtk_y_array)
  75.  
  76.     vtk_z_array = vtk.vtkDoubleArray()
  77.     vtk_z_array.SetName('z-coords')
  78.     vtk_z_array.SetArray(zs, len(zs), True)
  79.     vtk_box.SetZCoordinates(vtk_z_array)
  80.  
  81.     tally = np.array(tally_data)
  82.     tally_data = vtk.vtkDoubleArray()
  83.     tally_data.SetName(tally_label)
  84.     tally_data.SetArray(tally, tally.size, True)
  85.  
  86.     error = np.array(error_data)
  87.     error_data = vtk.vtkDoubleArray()
  88.     error_data.SetName("error_tag")
  89.     error_data.SetArray(error, error.size, True)
  90.  
  91.     vtk_box.GetCellData().AddArray(tally_data)
  92.     vtk_box.GetCellData().AddArray(error_data)
  93.  
  94.     writer = vtk.vtkRectilinearGridWriter()
  95.  
  96.     writer.SetFileName(outfile)
  97.  
  98.     writer.SetInputData(vtk_box)
  99.  
  100.     print('Writing %s' % outfile)
  101.  
  102.     writer.Write()
  103.  
  104. def main():
  105.  
  106.     ap = argparse.ArgumentParser(description=__doc__)
  107.  
  108.     ap.add_argument('-i','--input',
  109.                     required=True,
  110.                     help='Path to statepoint h5 file')
  111.  
  112.     ap.add_argument('-n','--tally-name',
  113.                     help='Tally name to add to mesh')
  114.  
  115.     ap.add_argument('-t','--tally-id',
  116.                     type=int,
  117.                     required=True,
  118.                     help='Tally name to add to mesh')
  119.  
  120.     ap.add_argument('-m','--mesh-id',
  121.                     type=int,
  122.                     required=True,
  123.                     help='ID of the mesh to convert')
  124.  
  125.     ap.add_argument('-o', '--output',
  126.                     action='store',
  127.                     default='meshtally.vtk',
  128.                     help='Name of outputfile (.h5m for MOAB, .vtk for VTK)')
  129.  
  130.     args = ap.parse_args()
  131.  
  132.     print('Loading file %s' % args.input)
  133.  
  134.     sp = openmc.StatePoint(args.input)
  135.  
  136.     print('Loading mesh with ID of %s' % args.mesh_id)
  137.     mesh = sp.meshes[args.mesh_id]
  138.  
  139.     xs = np.linspace(mesh.lower_left[0],
  140.                      mesh.upper_right[0],
  141.                      mesh.dimension[0] + 1)
  142.     ys = np.linspace(mesh.lower_left[1],
  143.                      mesh.upper_right[1],
  144.                      mesh.dimension[1] + 1)
  145.     zs = np.linspace(mesh.lower_left[2],
  146.                      mesh.upper_right[2],
  147.                      mesh.dimension[2] + 1)
  148.  
  149.     tally_args = {'name' : args.tally_name,
  150.                   'id'   : args.tally_id}
  151.  
  152.     msg = "Loading retrieving tally with \n"
  153.     if args.tally_name is not None:
  154.         msg += "Name: {}\n".format(args.tally_name)
  155.     if args.tally_id is not None:
  156.         msg += "ID: {}\n".format(args.tally_id)
  157.     try:
  158.         tally = sp.get_tally(**tally_args)
  159.     except LookupError as e:
  160.         raise e
  161.  
  162.     data = tally.mean[:,0,0]
  163.     error = tally.std_dev[:,0,0]
  164.  
  165.     data = data.tolist()
  166.     error = error.tolist()
  167.  
  168.     for c, i in enumerate(data):
  169.         if math.isnan(i):
  170.             data[c] = 0.
  171.  
  172.     for c, i in enumerate(error):
  173.         if math.isnan(i):
  174.             error[c] = 0.
  175.  
  176.     if args.tally_name is None:
  177.         tally_label = "tally_{}".format(args.tally_id)
  178.  
  179.     if args.output.endswith(".vtk"):
  180.         write_vtk(xs, ys, zs, tally_label, data, error, args.output)
  181.     else:
  182.         write_moab(xs, ys, zs, tally_label, data, error, args.output)
  183.  
  184. if __name__ == "__main__":
  185.     main()
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top