Allows input of OpenMC geometry

makeclean · Andrew Davis · commit 9d127e1eb749 · 2019-07-17T16:27:51.000+01:00
diff --git a/python/CellCard.py b/python/CellCard.py
@@ -27,7 +27,7 @@ def __init__(self,card_string):
         self.cell_material_number = 0
         self.cell_importance = 1 # note any importance - we assume everything else is 0
         self.cell_text_description = ""
-        self.cell_interpreted = ""
+        self.cell_interpreted = "" #this it the generalised form of the cell
         self.cell_fill = 0
         self.cell_universe = 0
         self.cell_universe_offset = 0
diff --git a/python/OpenMCCell.py b/python/OpenMCCell.py
@@ -107,9 +107,58 @@ def write_openmc_cell(cell, geometry_tree):
                       universe = str(universe))
 
 
+# take the xml attributes and populate the 
+# cell
+def cell_from_attribute(xml_attribute):
+    cell = OpenMCCell("")
+    cell.cell_id = xml_attribute["id"]
+ 
+    # todo if name based materials are used will need
+    # a helper function
+    if xml_attribute["material"] == "void":
+        cell.cell_material_number = 0
+    else:
+        cell.cell_material_number = xml_attribute["material"]
+        
+    cell.cell_text_description = xml_attribute["region"]
+    cell.cell_universe = xml_attribute["universe"]
+    cell.cell_fill = xml_attribute["fill"]
 
+    cell.generalise()
+    return cell
 
-#
+# base constructor
 class OpenMCCell(CellCard):
     def __init__(self, card_string):
         CellCard.__init__(self, card_string)
+
+    # turn the text representation of the cell
+    # into a generic description
+    def generalise(self):
+        # make an interable list of the components
+        cell_description = list(self.cell_text_description)
+        # first lets sanisise the text description  - remove
+        # double spaces trailing and leading white space
+        idx = 0
+        while True:
+            # breakout condition
+            if idx >= len(cell_description):
+                break
+            # part of the cell we're looking at
+            s = cell_description[idx]
+            if s is "|":
+                cell_description[idx] = CellCard.OperationType["UNION"]
+                idx += 1
+                continue
+            elif s is "~":
+                cell_description[idx] = CellCard.OperationType["NOT"]
+                idx += 1
+                continue
+            elif s is " ":
+                cell_description[idx] = CellCard.OperationType["AND"]
+                idx += 1
+                continue    
+            idx += 1
+        # set the generalised cell description
+        self.cell_interpreted = cell_description
+        return
diff --git a/python/OpenMCInput.py b/python/OpenMCInput.py
@@ -6,9 +6,9 @@
 import sys
 import xml.etree.ElementTree as ET
 
-from OpenMCSurface import write_openmc_surface
-from OpenMCCell import write_openmc_cell
-from OpenMCMaterial import write_openmc_material
+from OpenMCSurface import SurfaceCard,surface_from_attribute, write_openmc_surface
+from OpenMCCell import cell_from_attribute, write_openmc_cell
+from OpenMCMaterial import material_from_attribute, write_openmc_material
 '''
 copy and paste from http://effbot.org/zone/element-lib.htm#prettyprint
 it basically walks your tree and adds spaces and newlines so the tree is
@@ -35,8 +35,75 @@ class OpenMCInput(InputDeck):
     xml elements directly
     """
     # constructor
-    def __init__(self, filename = ""):
-        InputDeck.__init__(self, filename)
+    def __init__(self,geometry_file="geometry.xml",material_file="materials.xml"):
+        self.geometry = geometry_file
+        self.material = material_file
+        self.xml_geom = None
+        self.xml_material = None
+        # make a new input deck
+        InputDeck.__init__(self,filename="")
+
+    # read the openmc files
+    def read(self):
+        self.xml_geom = self.__read_geometry_xml(self.geometry)
+        self.xml_material = self.__read_material_xml(self.material)
+
+    # read the geometry file
+    def __read_geometry_xml(self, geometry_filename):
+        geom = ET.parse(geometry_filename)
+        geom_root = geom.getroot()
+        return geom_root
+
+    # read the material file
+    def __read_material_xml(self, material_filename):
+        mat = ET.parse(material_filename)
+        mat_root = mat.getroot()
+        return mat_root
+
+    # process the files
+    def process(self):
+        # do materials first 
+        self.__process_materials()
+        self.__process_geometry()
+        
+
+    # process the geometry
+    def __process_geometry(self):
+        for child in self.xml_geom:
+            if child.tag == "surface":
+                surface = surface_from_attribute(child.attrib)
+                InputDeck.surface_list.append(surface)
+            elif child.tag == "cell":
+                cell = cell_from_attribute(child.attrib)
+                InputDeck.cell_list.append(cell)
+
+        # loop over the cells and set the material 
+        # density for each cell
+        for cell in InputDeck.cell_list:
+            cell_mat = cell.cell_material_number
+            if cell_mat is not 0:
+                density = InputDeck.material_list[cell_mat].density
+                cell.cell_density = density
+
+        # identify the cells that have a surface with a vacuum
+        # boundary condition, they mark the end of the universe
+        boundary_surfs = []
+        for surface in InputDeck.surface_list:
+            if surface.boundary_condition == SurfaceCard.BoundaryCondition["VACUUM"]:
+                boundary_surfs.append(surface.surface_id)
+
+        # now loop through the cells and mark as importance 0 where
+        # appropriate
+        for cell in InputDeck.cell_list:
+            if set(boundary_surfs) & set(cell.cell_interpreted) == set(boundary_surfs):
+                cell.cell_importance = 0
+
+    # process the materials
+    def __process_materials(self):
+        for child in self.xml_material:
+            if child.tag == "material":
+                material = material_from_attribute(child.attrib, child.getchildren())
+                InputDeck.material_list[material.material_number] = material
 
     # write the collection of OpenMC surface definitions
     def __write_openmc_surfaces(self, geometry_tree):
diff --git a/python/OpenMCMaterial.py b/python/OpenMCMaterial.py
@@ -1,5 +1,6 @@
 #!/usr/env/python3
 
+import re
 from MaterialCard import MaterialCard
 import xml.etree.ElementTree as ET
 
@@ -24,6 +25,26 @@
              109:"Mt",110:"Ds",111:"Rg",112:"Cn",113:"Nh",
              114:"Fl",115:"Mc",116:"Lv",117:"Ts",118:"Og"}
 
+# return the zz given a name
+def get_zaid(name):
+    m = re.search('[A-Z]?[a-z]*',name)
+    element_name = m.group(0)
+    m = re.search('[0-9]?[0-9]?[0-9]',name)
+    nucleon_number = int(m.group(0))
+
+    for zz,el in name_zaid.items():
+        if el == element_name:
+            break
+
+    # pad the right number of zeros
+    if nucleon_number < 10:
+        nucleon_number = "00" + str(nucleon_number)
+    elif nucleon_number > 10 and nucleon_number < 100:
+        nucleon_number = "0" + str(nucleon_number)  
+        
+
+    return str(zz)+nucleon_number
+
 # convert zaid to a name for openmc
 def zaid_to_name(zaid_string):
     if len(zaid_string) <= 4:
@@ -41,6 +62,43 @@ def zaid_to_name(zaid_string):
     name = name_zaid[zz]
     return name+str(aa)
 
+################## input functions #########################
+# 
+def material_from_attribute(xml_element, children):
+    # 
+    material = MaterialCard()
+    material.material_number = xml_element["id"]
+    # loop over the constituents
+    nucs = {}
+    atom_fraction = False
+
+    for child in children:
+        # set the density
+        if "units" in child.attrib:
+            units = child.attrib["units"]
+            density = float(child.attrib["value"])
+        # set the material parameters
+        if "name" in child.attrib:
+            nucid = child.attrib["name"]
+            zaid = get_zaid(nucid)
+            if "wo" in child.attrib:     
+                # mass fractions are -ve 
+                nucs[zaid] = -1*float(child.attrib["wo"])
+            if "ao" in child.attrib:   
+                atom_fraction = True   
+                nucs[zaid] = child.attrib["ao"]      
+
+    if not atom_fraction: density = density*-1.0
+
+    material.density = density
+    material.density_units = units
+    material.composition_dictionary = nucs
+
+    return material
+
+
+################## output functions #########################
+
 # write the atomic fraction entry
 def __write_atomic_fraction(material, nuclide, mass_frac):       
     ET.SubElement(material, "nuclide", name = nuclide, ao = str(abs(mass_frac)))
diff --git a/python/OpenMCSurface.py b/python/OpenMCSurface.py
@@ -5,6 +5,15 @@
 
 import warnings
 
+def boundarystring_to_type(boundary_string):
+    if boundary_string == "transmission":
+        boundary_type = SurfaceCard.BoundaryCondition["TRANSMISSION"]
+    if boundary_string == "vacuum":
+        boundary_type = SurfaceCard.BoundaryCondition["VACUUM"]
+    if boundary_string == "reflecting":
+        boundary_type = SurfaceCard.BoundaryCondition["REFLECTING"]
+    return boundary_type
+
 def boundary_condition(boundaryCondition):
     if boundaryCondition == SurfaceCard.BoundaryCondition["TRANSMISSION"]:
         boundary = "transmission"
@@ -18,6 +27,34 @@ def boundary_condition(boundaryCondition):
 
     return boundary
 
+def type_to_generictype(surface_string):
+    if surface_string == "sphere":
+        surface = SurfaceCard.SurfaceType["SPHERE_GENERAL"]
+    if surface_string == "x-plane":
+        surface = SurfaceCard.SurfaceType["PLANE_X"]
+    if surface_string == "y-plane":
+        surface = SurfaceCard.SurfaceType["PLANE_Y"]
+    if surface_string == "z-plane":
+        surface = SurfaceCard.SurfaceType["PLANE_Z"]
+    if surface_string == "plane":
+        surface = SurfaceCard.SurfaceType["PLANE_GENERAL"]
+    if surface_string == "x-cylinder":
+        surface = SurfaceCard.SurfaceType["CYLINDER_X"]
+    if surface_string == "y-cylinder":
+        surface = SurfaceCard.SurfaceType["CYLINDER_Y"]
+    if surface_string == "z-cylinder":
+        surface = SurfaceCard.SurfaceType["CYLINDER_Z"]
+    if surface_string == "x-cone":
+        surface = SurfaceCard.SurfaceType["CONE_X"]
+    if surface_string == "y-cone":
+        surface = SurfaceCard.SurfaceType["CONE_Y"]
+    if surface_string == "z-cone":
+        surface = SurfaceCard.SurfaceType["CONE_Z"]
+    if surface_string == "quadric":
+        surface = SurfaceCard.SurfaceType["GENERAL_QUADRATIC"]
+
+    return surface
+
 def openmc_surface_info(SurfaceCard):
     if SurfaceCard.surface_type == SurfaceCard.SurfaceType["PLANE_GENERAL"]:
         type_string = "plane"        
@@ -77,7 +114,16 @@ def write_openmc_surface(SurfaceCard, geometry_tree):
     ET.SubElement(geometry_tree, "surface", id = str(id), type = str(type),
                   coeffs = str(coeffs), 
                   boundary = boundary_condition(SurfaceCard.boundary_condition))
-    
+
+# import openmc surface functions
+def surface_from_attribute(xml_attribute):
+    surface = SurfaceCard("")
+    # loop over the surface attributes and build the generic description
+    surface.boundary_condition = boundarystring_to_type(xml_attribute['boundary']) 
+    surface.surface_coefficients = xml_attribute['coeffs'].split() 
+    surface.surface_id = xml_attribute['id'] 
+    surface.surface_type = type_to_generictype(xml_attribute['type'])
+    return surface
     
 class OpenMCSurfaceCard(SurfaceCard):
     """ Class to handle the creation and translation of
