WIP: MCNP additions

csg2csg/CellCard.py

@@ -35,6 +35,10 @@ def __init__(self,card_string):
         self.cell_universe_transformation_id = "0" # if there is a cell_universe tr number it should be purged
                                         # and converted into an offset and rotation
         self.cell_surface_list = set() # list of surfaces used in the cell definition
+        self.cell_lattice = None
+        self.cell_lattice_type = None
+        self.cell_volume = None
+        self.cell_temperature = None
 
         Card.__init__(self,card_string)
 

csg2csg/Input.py

@@ -17,6 +17,7 @@ def __init__(self, filename):
         # TODO maybe these should be dictionaries by index
         self.cell_list = []
         self.surface_list = []
+        self.lattice_list = []
         self.last_free_surface_index = 0
         self.importance_list = {} # dictionary of importances
         self.material_list = {}
@@ -70,6 +71,7 @@ def from_input(self,InputDeckClass):
         self.title = InputDeckClass.title
         self.cell_list = InputDeckClass.cell_list
         self.surface_list = InputDeckClass.surface_list
+        self.lattice_list = InputDeckClass.lattice_list
         self.material_list = InputDeckClass.material_list
         return
 

csg2csg/MCNPCellCard.py

@@ -10,8 +10,14 @@
 
 import logging
 
+import numpy
+import scipy.constants
+
+# define constants
+k_Mev = scipy.constants.physical_constants.get("Boltzmann constant in eV/K")[0]/1.e6
+
 # to support more keywords for cells add them here
-mcnp_cell_keywords = ["imp","u","fill","vol", "tmp"]
+mcnp_cell_keywords = ["imp","u","fill","vol", "tmp", "lat"]
 
 # if the string is a cell card or not
 def is_cell_card(line):
@@ -86,6 +92,12 @@ def write_mcnp_cell(filestream, CellCard, print_importances = True):
     if CellCard.cell_universe != 0:
         string += " u=" + CellCard.cell_universe
 
+    if CellCard.cell_volume is not None:
+        string += " vol=" + CellCard.cell_volume
+
+    if CellCard.cell_temperature is not None:
+        string += " tmp=" + str(CellCard.cell_temperature*k_Mev)
+
     if CellCard.cell_fill != 0:
         string += " fill="+CellCard.cell_fill + " "
         if CellCard.cell_universe_offset != 0 or CellCard.cell_universe_rotation != 0:
@@ -143,7 +155,10 @@ def generalise(self):
         # print(cell_description)
         # cell_description = [ s = "" for item in cell_description if item == "+"]
         cell_description = list(cell_description)
-
+
+        # need to reset cell surface list, to exclude removed macrobody, etc.
+        self.cell_surface_list = set()
+
         idx = 0
         while True:
             s = cell_description[idx]
@@ -158,7 +173,7 @@ def generalise(self):
             elif s == ("(" or  ")"):
                 idx += 1
                 continue
-            elif isinstance(s,str) and cell_description[idx-1] != "(" and cell_description[idx] != ")":
+            elif isinstance(s,str) and cell_description[idx] != "(" and cell_description[idx] != ")":
                 cell_description.insert(idx,CellCard.OperationType["AND"])
                 idx += 1
                 try:
@@ -187,7 +202,7 @@ def __sanitise(self):
     # keyword
     def __get_keyword_value(self,keyword,string):
         #regex = re.regex=re.compile("("+keyword+") ?= ?[1-9][0-9]*")
-        regex = re.regex=re.compile("("+keyword+") ?= ?(?=.)([+-]?([0-9]*)(\.([0-9]+))?)")
+        regex = re.compile("("+keyword+") ?= ?(?=.)([+-]?([0-9]*)(\.([0-9]+))?)")
         result = regex.search(string)[0]
         return result.split(" ")[2] #string[offset:end]
 
@@ -232,10 +247,11 @@ def __detect_keywords(self, keywords, string):
 
         posi = string.find("imp")
         posv = string.find("vol")
+        posl = string.find("lat")
 
         # find the posititon of the first match
-        positions = [posu, posf, posi, posv, post]
-        if posu != -1 or posf != -1 or posi != -1 or posv != -1 or post != -1:
+        positions = [posu, posf, posi, posv, post, posl]
+        if posu != -1 or posf != -1 or posi != -1 or posv != -1 or post != -1 or posl != -1:
             m = min(i for i in positions if i > 0)
         else:
             return string
@@ -252,23 +268,55 @@ def __detect_keywords(self, keywords, string):
         else:
             self.cell_universe = self.__get_keyword_value('u',end_of_string)
 
+        if posv != -1:
+            self.cell_volume = self.__get_keyword_value('vol',end_of_string)
+
+        if post != -1:
+            regex = re.compile("tmp ?= \d+\.\d+(?:[eE][+\-]?\d+)?")
+            result = regex.search(end_of_string)[0].split("=")[-1].strip()
+            self.cell_temperature = float(result)/k_Mev
+
         if posf == -1:
             self.cell_fill = 0
         else:
-            self.cell_fill = self.__get_keyword_value('fill',end_of_string).strip()
-            # if we have found fill, there may also be a rotation and translation
-            # associated with the universe of the form (0 0 0)
-            if '(' in string[posf:]:
-                rot_trans = self.__extract_string_between(string[posf:],'(',')')
-            else:
-                rot_trans = "0"
-
-            self.__set_universe_transform(rot_trans,rot_angle_degrees)
+            # if lattice, assume fully specified fill cell card
+            if posl != -1:
+               # get the lattice element bounds
+               regex = re.compile("fill ?= ([+-]?([0-9]*)):([+-]?([0-9]*)) ([+-]?([0-9]*)):([+-]?([0-9]*)) ([+-]?([0-9]*)):([+-]?([0-9]*))")
+               result = regex.search(end_of_string)[0]
+
+               self.cell_fill = result.split("=")[-1].strip()
+
+               # process the bounds
+               i, j, k  = self.cell_fill.split()
+               i_coords = list(range(int(i.split(":")[0]),int(i.split(":")[1])+1))
+               j_coords = list(range(int(j.split(":")[0]),int(j.split(":")[1])+1))
+               k_coords = list(range(int(k.split(":")[0]),int(k.split(":")[1])+1))
+               total_cells = len(i_coords)*len(j_coords)*len(k_coords)
+
+               # now read in the lattice
+               lattice = string[posf + len(result):].strip().split()[:total_cells]
+               # note this is reshaped in to z,x,y due to how numpy tiles/reshapes the data
+               self.cell_lattice = numpy.array(lattice).reshape(len(k_coords),len(i_coords),len(j_coords)).astype(int)
 
+            else:
+                self.cell_fill = self.__get_keyword_value('fill',end_of_string).strip()
+                # if we have found fill, there may also be a rotation and translation
+                # associated with the universe of the form (0 0 0)
+                if '(' in string[posf:]:
+                    rot_trans = self.__extract_string_between(string[posf:],'(',')')
+                else:
+                    rot_trans = "0"
+
+                self.__set_universe_transform(rot_trans,rot_angle_degrees)
+
         if posi == -1:
             self.cell_importance = 1.
         else:
-            self.cell_importance = float(self.__get_keyword_value('imp:n',end_of_string))
+            self.cell_importance = float(self.__get_keyword_value("imp:[npe|quvfhl+-xyo!<>g/zk%^b_~cw@dtsa*?#,]+",end_of_string))
+
+        if posl != -1:
+            self.cell_lattice_type = int(self.__get_keyword_value("lat",end_of_string))
 
         # return the string upto the posisiotn of the first detected keyword
         return string[:m]

csg2csg/MCNPInput.py

@@ -9,7 +9,7 @@
 from csg2csg.MCNPCellCard import MCNPCellCard, is_cell_card, write_mcnp_cell
 from csg2csg.MCNPSurfaceCard import MCNPSurfaceCard, is_surface_card, write_mcnp_surface
 from csg2csg.MCNPDataCard import MCNPTransformCard
-from csg2csg.MCNPMaterialCard import MCNPMaterialCard, write_mcnp_material
+from csg2csg.MCNPMaterialCard import MCNPMaterialCard, MCNPSABCard, write_mcnp_material
 
 from collections import Counter
 
@@ -23,6 +23,7 @@
 
 import warnings
 import logging
+import os
 import sys
 import re
 
@@ -39,6 +40,35 @@ def __init__(self, filename =""):
     # TODO - maybe make a function that aribitrarily extract text
     # between one keyword until another keyword is found
 
+    # extend read method to process "read" cards in MCNP input file
+    def read(self):
+        with open(self.filename, 'rU', errors="replace") as f:
+            file_content = f.read()
+
+        # look for and swap read cards
+        read_cards = re.findall("^read .*?\n", file_content, re.DOTALL|re.MULTILINE)
+        for read_card in read_cards:
+            filepath = read_card.split("=")[1].strip().split()[0]
+            # if file exists, read it in and replace
+            if os.path.isfile(filepath):
+                with open(filepath, 'rU', errors="replace") as i:
+                    read_input = i.read()
+                i.close()
+
+                file_content = re.sub(read_card, read_input, file_content)
+
+            # else, print error and comment out
+            else:
+                print("WARNING: read file {} does not exist.".format(filepath))
+                file_content = re.sub(read_card, "c " + read_card, file_content)
+
+        self.file_lines = file_content.splitlines(keepends=True)
+
+        # split into lines and count total lines
+        self.file_lines = [x.lower() for x in self.file_lines]
+
+        self.total_num_lines = len(self.file_lines)
+
     def __set_title(self):
         # set the title card
         if "message" in self.file_lines[0]:
@@ -189,12 +219,24 @@ def __get_material_card(self, start_line):
                 idx += 1
             break
 
-        material = MCNPMaterialCard(mat_num, material_string)
-        # set the colour based on the number of colours
-        # but only if its really used rather than a tally
-        # multiplier material
-        material.material_colour = get_material_colour(len(self.material_list))
-        self.material_list[material.material_number] = material
+        # process s(alpha, beta) 
+        # assumes that all s(alpha, beta) are input after corresponding material cards in the input file
+        if "t" in mat_num:
+            mat_num = mat_num.replace("t","")
+            material = MCNPSABCard(mat_num, material_string)
+
+            # update the previous material card to include thermal scattering
+            self.material_list[material.material_number].thermal_scattering = material.thermal_scattering
+
+            return 
+
+        else:
+            material = MCNPMaterialCard(mat_num, material_string)
+            # set the colour based on the number of colours
+            # but only if its really used rather than a tally
+            # multiplier material
+            material.material_colour = get_material_colour(len(self.material_list))
+            self.material_list[material.material_number] = material
 
         return
 
@@ -209,7 +251,7 @@ def __get_material_cards(self, start_line):
 
             # this crazy makes sure that we find an "m" in the line but that we dont
             # find another keyword with an m in it like prdmp
-            if re.match(" *m[0-9]/*",self.file_lines[idx]):
+            if re.match(" *m(t)?[0-9]/*",self.file_lines[idx]):
 #            if "m" in self.file_lines[idx] and not any(x in self.file_lines[idx] for x in mcnp_keywords):
                 logging.debug("%s", "material found on line " + str(idx))
                 self.__get_material_card(idx)
@@ -670,6 +712,42 @@ def __flatten_macrobodies(self):
                     text_string = ' '.join(cell.cell_text_description)
                     self.cell_list[jdx].update(text_string)
 
+                # update the lattice definition (cell in mcnp) - loop over all cells
+                for jdx, cell in enumerate(self.lattice_list):
+                    while True:
+                        # cell text description is contually updated
+                        cell_text_description = cell.cell_text_description
+
+                        # if we find the surface id of the macrobdy in the text description
+                        sub = str(surf.surface_id)
+                        regex = re.compile("^-?("+str(surf.surface_id)+")(\.+[1-9])?$")
+                        matches = [m.group(0) for l in cell_text_description for m in [regex.search(l)] if m]
+                        #if str(surf.surface_id) in cell_text_description or str(surf.surface_id)+"." in cell_text_description:
+
+                        if matches:                       
+                            # loop over each component and find the macrobody
+                            for idx, surface in enumerate(cell.cell_text_description):
+                                # if it matches we have the simmple form
+                                if str(surf.surface_id) == surface:
+                                    # replace it
+                                    cell.cell_text_description[idx] = cell_description[1]
+                                elif "-"+str(surf.surface_id) == surface:
+                                    cell.cell_text_description[idx] = cell_description[0]
+
+                                # else we have the facet form
+                                if str(surf.surface_id)+"." in surface:                                    
+                                    surface_index = int(surface.split(".")[1]) # get just the mcnp surface index
+                                    new_surface_id = new_surfaces[surface_index-1].surface_id # mcnp numbers them 1->n
+                                    if "-" in surface: # need to take care of the -sign
+                                        cell.cell_text_description[idx] = "-"+str(new_surface_id)
+                                    else:
+                                        cell.cell_text_description[idx] = str(new_surface_id)
+                        else:
+                            break
+                    # update the text description
+                    text_string = ' '.join(cell.cell_text_description)
+                    self.lattice_list[jdx].update(text_string)
+
         # clear up removed surfaces
         logging.debug("%s", "Deleting macrobody surfaces")
         for surf in to_remove:
@@ -825,11 +903,19 @@ def __get_cell_cards(self):
                 # mcnp continue line is indicated by 5 spaces
                 if cell_line.startswith("     ") and not cell_line.isspace():
                     card_line += cell_line
+                # if & continuation 
+                elif self.file_lines[jdx-1].rstrip().endswith("&"):
+                    card_line += cell_line
+                    # strip any &
+                    card_line = re.sub("&", "", card_line)
                 else: # else we have found a new cell card
                     logging.debug("%s\n", "Found new cell card " + card_line)
                     cellcard = MCNPCellCard(card_line)
-                    # we should set the comment here
-                    self.cell_list.append(cellcard)
+                    if cellcard.cell_lattice_type is not None:
+                        self.lattice_list.append(cellcard)
+                    else:
+                        # we should set the comment here
+                        self.cell_list.append(cellcard)
                     break
                 jdx += 1
             idx = jdx
@@ -853,6 +939,7 @@ def __apply_boundary_conditions(self):
         for cell in self.cell_list:
             if cell.cell_importance == 0:
                 for surf in cell.cell_surface_list:
+                    if self.get_surface_with_id(surf) is None: continue
                     self.get_surface_with_id(surf).boundary_condition = SurfaceCard.BoundaryCondition["VACUUM"]
         return
 
@@ -884,6 +971,44 @@ def __get_surface_cards(self,idx):
             idx = jdx
         return idx
 
+    # return specific surface card
+    def __get_surface_card(self, surface_id):
+        for surface in self.surface_list:
+            if surface.surface_id == surface_id:
+                return surface
+
+    # update lattices
+    def __update_lattices(self):
+        # todo - note, this is not very robust
+        # "orientation" is only used for hex lattice
+        for lattice in self.lattice_list:
+            # initialize data
+            heights = []
+
+            # determine orientation, flat-to-flat, and height
+            for surface_id in lattice.cell_surface_list:
+                surface = self.__get_surface_card(surface_id)
+
+                # if z
+                if surface.surface_type == SurfaceCard.SurfaceType["PLANE_Z"]:
+                    heights.append( surface.surface_coefficients[3] )
+                elif (surface.surface_type == SurfaceCard.SurfaceType["PLANE_GENERAL"] and 
+                      surface.surface_coefficients[0] == 0 and
+                      surface.surface_coefficients[1] == 0 and
+                      surface.surface_coefficients[2] == -1):
+                    heights.append( -surface.surface_coefficients[3] )
+                # if px, py
+                elif surface.surface_type == SurfaceCard.SurfaceType["PLANE_X"]:
+                    orientation = "x"
+                    pitch = abs(2.*surface.surface_coefficients[3])
+                elif surface.surface_type == SurfaceCard.SurfaceType["PLANE_Y"]:
+                    orientation = "y"
+                    pitch = abs(2.*surface.surface_coefficients[3])
+
+            height = abs(heights[0] - heights[1])
+
+            lattice.pitch = [pitch, height]
+            lattice.orientation = orientation
     # process the mcnp input deck and read into a generic datastructure
     # that we can translate to other formats
     def process(self):
@@ -948,6 +1073,9 @@ def process(self):
         # update the bounding coordinates of surfaces that need it
         # cones for example
         self.__update_surfaces()
+
+        # process lattices, to get geometric information, now that surfaces are defined
+        self.__update_lattices()
 
         self.split_unions()