Commit 57bc3559 authored by Gurecky, William's avatar Gurecky, William
Browse files

allow Core instances to delegate CTF model spec tasks to contained model instance

parent 986a5e55
Loading
Loading
Loading
Loading
+32 −2
Original line number Diff line number Diff line
@@ -22,7 +22,7 @@ import Section
class Core(object):
   """ Base class for building core maps."""
   def __init__(me):
      pass
      me._model = None

   def getSection(me):
      """ Returns the Section object that contains the core. """
@@ -40,6 +40,10 @@ class Core(object):
      """ See overriding procedure doc"""
      raise NotImplementedError()

   def setEditOptions(me, *args, **kwrags):
      """ Set output edit options """
      raise NotImplementedError()

   def getCoreMaxRadialBound(me):
      """ Returns the furthest point in the model from the core origin."""
      coreRadius = 0.0
@@ -59,7 +63,22 @@ class Core(object):
      """ Returns a list of the channel IDs of the channels in the core region."""
      return me.core_ch

  ### Core Attributes ###
  ### Core Base Attributes ###
   @property
   def model(me):
      """ Model container """
      return me._model

   @model.setter
   def model(me, my_model):
      """
      Set Model container
      Args:
        my_model: instance of Model.Model
      """
      assert isinstance(my_model, Model.Model)
      me._model = my_model

   @property
   def dz(me):
      """ Mesh partition heights """
@@ -105,6 +124,17 @@ class Core(object):
      z_pairs = np.asarray((me.z_bounds[1:], me.z_bounds[:-1])).T
      return z_pairs

   # === Delegate to me.model ===
   def __getattr__(me, attr):
      """
      Delegates methods called on this Core builder which do not match any method currently
      implemented in this builder to the contained Model instance.

      Args:
          attr:  python method or attribute of the Model class
      """
      return getattr(me.model, attr)


class MSRE(Core):
   """ Core extension for building core geometry map for the MSRE.
+30 −24
Original line number Diff line number Diff line
@@ -616,7 +616,7 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core):
   """
   def __init__(me, model, coreMap, assemMaps, pitch, dz, solidGeoms, startChanIndex=1, assemPitch=None,
            baffleGap=0.0, sectionID=1, symOpt=1, assemLosses=None, assemLossLevels=None, rodDomains=None,
            modelGuideTubes=False, setCoreMaps=False):
            modelGuideTubes=False):
      assert(isinstance(model, Model.Model))
      assert(len(coreMap)>0)
      assert(utils.isArraySquare(coreMap))
@@ -804,9 +804,25 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core):
      me.chMapObj = chMapObj
      me.nodalGeom = nodalGeom
      me.dz = dz
      if setCoreMaps:

   def setEditOptions(me, chanVTK=None, rodVTK=None, ctfHDF5=None, veracsHDF5=None, chanASCII=None,
         rodEdits=None, dnbEdits=None):
      """
      Delegates setting edits to the Model instance contained in this builder.

      Args:
          chanVTK (bool): Set to True to write the channel VTK file (default True).
          rodVTK (bool): Set to True to write the rod VTK file (default False).
          ctfHDF5 (bool): Set to True to write the native CTF HDF5 file (default True).
          veracsHDF5 (bool): Set to True to write the VERA-CS HDF5 file (requires that core map be specified and will only work for square rod lattices) (default False).
          chanASCII (bool): Set to True to write the ASCII channel output file.
          rodEdits (bool): Set to True to write rod edit information to the standard CTF .out file (Note the rod information in the .out file is limited in scope and it is better to get rod data from the native HDF5 file).
          dnbEdits (bool): Set  to True to produce the CTF .dnb file, which contains information about DNB (default False).
      """
      if veracsHDF5 == True:
         rodMap, chanMap, assemMap, symOpt = me.getChannelAndRodMaps()
         me.model.setCoreMaps(rodMap, chanMap, assemMap, symOpt)
      me.model.setEditOptions(chanVTK, rodVTK, ctfHDF5, veracsHDF5, chanASCII, rodEdits, dnbEdits)

   def getChannelAndRodMaps(me):
      """
@@ -820,16 +836,18 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core):
      >>> my_lwr_nodal.model.setEditOptions(veracsHDF5=True)

      Returns:
            rodMap, chanMap, assemMap, symOpt:
            where rodMap is a 2 dimensional list
            where chanMap is a 2 dimensional list
            where assemMap is a 2 dimensional list
          rodMap, chanMap, assemMap, symOpt
          rodMap is a 2 dimensional list
          chanMap is a 2 dimensional list
          assemMap is a 2 dimensional list
          symOpt is an int: 1=full sym, 4=qtr sym
      """
      if not me.modelGuideTubes:
          rodMap = np.zeros((len(me.chMapObj.chAssemIdx),len(me.chMapObj.chAssemIdx[0])), dtype=int)
      else:
      if me.modelGuideTubes:
          raise RuntimeError("Model guide tubes is currently not supported " + \
                             "by the nodal veracsHDF5 writer. Set modelGuideTubes=False.")
          rodMap = np.zeros((len(me.chMapObj.chAssemIdx),2*len(me.chMapObj.chAssemIdx[0])), dtype=int)
      else:
          rodMap = np.zeros((len(me.chMapObj.chAssemIdx),len(me.chMapObj.chAssemIdx[0])), dtype=int)
      symOpt = 1
      for obj in me.model.solidObjects.values():
         # this only works in full sym mode for now
@@ -845,17 +863,6 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core):
               if me.modelGuideTubes: assert len(rods) == 8
               if not me.modelGuideTubes: assert len(rods) == 4

               tmp_rod_types = []
               tmp_rod_x = []
               tmp_rod_y = []
               for rodID in rods:
                  obj_rod = me.model.solidObjects[rodID]
                  obj_rod_type_id = me.model.solidTypeIDs[rodID]
                  tmp_rod_types.append(obj_rod_type_id)
                  tmp_rod_x.append(obj_rod.x)
                  # in nodal model: xy loc of guide tube == xy locs of heated rods
                  tmp_rod_y.append(obj_rod.y)

               # fill rodMap
               if me.modelGuideTubes:
                   # top nodes in asm
@@ -877,7 +884,6 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core):
                   rodMap[asm_i * 2 + 1][asm_j * 2 + 1] = int(assemID)
      # convert rodMap to list of lists to be consistent with SubKit
      rodMap = list([list(rodRow) for rodRow in rodMap])

      return rodMap, me.chMapObj.chAssemIdx, me.chMapObj.assemIdx, symOpt

   def setCoreInletBC(me, massflux, temperature, flowRamp=None):
+49 −152
Original line number Diff line number Diff line
@@ -77,87 +77,27 @@ def main(power_3d_flag=1):

    # Fuel rod is type 1, but mark it with an "f" to make it more visible
    # Guide tube is type 2
    f = 1
    f, g = 1, 2
    pwr17x17 = [
        [f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f],
        [f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f],
        [f, f, f, f, f, 2, f, f, 2, f, f, 2, f, f, f, f, f],
        [f, f, f, 2, f, f, f, f, f, f, f, f, f, 2, f, f, f],
        [f, f, f, f, f, g, f, f, g, f, f, g, f, f, f, f, f],
        [f, f, f, g, f, f, f, f, f, f, f, f, f, g, f, f, f],
        [f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f],
        [f, f, 2, f, f, 2, f, f, 2, f, f, 2, f, f, 2, f, f],
        [f, f, g, f, f, g, f, f, g, f, f, g, f, f, g, f, f],
        [f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f],
        [f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f],
        [f, f, 2, f, f, 2, f, f, 2, f, f, 2, f, f, 2, f, f],
        [f, f, g, f, f, g, f, f, g, f, f, g, f, f, g, f, f],
        [f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f],
        [f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f],
        [f, f, 2, f, f, 2, f, f, 2, f, f, 2, f, f, 2, f, f],
        [f, f, g, f, f, g, f, f, g, f, f, g, f, f, g, f, f],
        [f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f],
        [f, f, f, 2, f, f, f, f, f, f, f, f, f, 2, f, f, f],
        [f, f, f, f, f, 2, f, f, 2, f, f, 2, f, f, f, f, f],
        [f, f, f, g, f, f, f, f, f, f, f, f, f, g, f, f, f],
        [f, f, f, f, f, g, f, f, g, f, f, g, f, f, f, f, f],
        [f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f],
        [f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f]]

    # The multiple gap connects are not working with CTF parallel for some reason
    rodDomains = [[0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 0, 0, 0, 0],
                  [0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 0, 0],
                  [0, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 0],
                  [0, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 0],
                  [1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2],
                  [1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2],
                  [1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2],
                  [3, 3, 3, 3, 3, 3, 3, 4, 2, 2, 2, 2, 2, 2, 2],
                  [3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4],
                  [3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4],
                  [3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4],
                  [0, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 0],
                  [0, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 0],
                  [0, 0, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 0, 0],
                  [0, 0, 0, 0, 3, 3, 3, 4, 4, 4, 4, 0, 0, 0, 0]]
    # rodDomains = [[ 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 0, 0, 0, 0],
    #              [ 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 0, 0],
    #              [ 0, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 0],
    #              [ 0, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 0],
    #              [ 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2],
    #              [ 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2],
    #              [ 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2],
    #              [ 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2],
    #              [ 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2],
    #              [ 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2],
    #              [ 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2],
    #              [ 0, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 0],
    #              [ 0, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 0],
    #              [ 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 0, 0],
    #              [ 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 0, 0, 0, 0]]
    # rodDomains = [[ 0, 0, 0, 0, 1, 1, 2, 2, 2, 3, 3, 0, 0, 0, 0],
    #              [ 0, 0, 1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 3, 0, 0],
    #              [ 0, 1, 1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 3, 3, 0],
    #              [ 0, 1, 1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 3, 3, 0],
    #              [ 1, 1, 1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 3, 3, 3],
    #              [ 1, 1, 1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 3, 3, 3],
    #              [ 1, 1, 1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 3, 3, 3],
    #              [ 1, 1, 1, 1, 1, 1, 2, 2, 2, 6, 6, 6, 6, 6, 6],
    #              [ 4, 4, 4, 4, 4, 4, 5, 5, 5, 6, 6, 6, 6, 6, 6],
    #              [ 4, 4, 4, 4, 4, 4, 5, 5, 5, 6, 6, 6, 6, 6, 6],
    #              [ 4, 4, 4, 4, 4, 4, 5, 5, 5, 6, 6, 6, 6, 6, 6],
    #              [ 0, 4, 4, 4, 4, 4, 5, 5, 5, 6, 6, 6, 6, 6, 0],
    #              [ 0, 4, 4, 4, 4, 4, 5, 5, 5, 6, 6, 6, 6, 6, 0],
    #              [ 0, 0, 4, 4, 4, 4, 5, 5, 5, 6, 6, 6, 6, 0, 0],
    #              [ 0, 0, 0, 0, 4, 4, 5, 5, 5, 6, 6, 0, 0, 0, 0]]
    # rodDomains = [[ 0, 0, 0, 0, 1, 2, 2, 2, 2, 3, 3, 0, 0, 0, 0],
    #              [ 0, 0, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 0, 0],
    #              [ 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 0],
    #              [ 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 0],
    #              [ 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3],
    #              [ 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3],
    #              [ 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3],
    #              [ 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 6, 6, 6, 6, 6],
    #              [ 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6],
    #              [ 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6],
    #              [ 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6],
    #              [ 0, 4, 4, 4, 4, 5, 5, 5, 5, 5, 6, 6, 6, 6, 0],
    #              [ 0, 4, 4, 4, 4, 5, 5, 5, 5, 5, 6, 6, 6, 6, 0],
    #              [ 0, 0, 4, 4, 4, 5, 5, 5, 5, 5, 6, 6, 6, 0, 0],
    #              [ 0, 0, 0, 0, 4, 5, 5, 5, 5, 5, 6, 0, 0, 0, 0]]
    # Parallel decomposition map
    rodDomains = [[0, 0, 0, 0, 1, 2, 2, 2, 2, 2, 3, 0, 0, 0, 0],
                  [0, 0, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 0, 0],
                  [0, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 0],
@@ -173,111 +113,68 @@ def main(power_3d_flag=1):
                  [0, 4, 4, 4, 4, 5, 5, 5, 5, 5, 6, 6, 6, 6, 0],
                  [0, 0, 4, 4, 4, 5, 5, 5, 5, 5, 6, 6, 6, 0, 0],
                  [0, 0, 0, 0, 4, 5, 5, 5, 5, 5, 6, 0, 0, 0, 0]]
    # rodDomains = [[ 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 0, 0, 0, 0],
    #              [ 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 0, 0],
    #              [ 0, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 0],
    #              [ 0, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 0],
    #              [ 1, 1, 1, 1, 1, 1, 5, 5, 5, 2, 2, 2, 2, 2, 2],
    #              [ 1, 1, 1, 1, 1, 5, 5, 5, 5, 5, 2, 2, 2, 2, 2],
    #              [ 1, 1, 1, 1, 5, 5, 5, 5, 5, 5, 5, 2, 2, 2, 2],
    #              [ 3, 3, 3, 3, 5, 5, 5, 5, 5, 5, 5, 2, 2, 2, 2],
    #              [ 3, 3, 3, 3, 5, 5, 5, 5, 5, 5, 5, 4, 4, 4, 4],
    #              [ 3, 3, 3, 3, 3, 5, 5, 5, 5, 5, 4, 4, 4, 4, 4],
    #              [ 3, 3, 3, 3, 3, 3, 5, 5, 5, 4, 4, 4, 4, 4, 4],
    #              [ 0, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 0],
    #              [ 0, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 0],
    #              [ 0, 0, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 0, 0],
    #              [ 0, 0, 0, 0, 3, 3, 3, 4, 4, 4, 4, 0, 0, 0, 0]]

    # crate assembly map for vera HDF5 edits
    asmMap = np.zeros((15,15), dtype=int)
    n_asm = 0
    for i in range(asmMap.shape[0]):
        for j in range(asmMap.shape[1]):
            if coreMap[i][j] >= 1:
                n_asm += 1
            asmMap[i, j] = n_asm
    asmMap = [list(asmRow) for asmRow in asmMap]

    # crate rod map for vera HDF5 edits
    pass

    # crate channel map for vera HDF5 edits
    pass

    print np.unique(np.array(rodDomains), return_counts=True)
    # print number of assemblies associated with each parallel process
    print("=== Parallel Zones, N Asm in Zone ===")
    print(np.unique(np.array(rodDomains), return_counts=True))

    # initilize the CTF Model
    model = cdb.Model()

    model.setTitle("4-Loop nodal mesh loss of flow")

    solidGeoms = {1: cdb.FuelRod(d_outer=dFuelOuter, d_inner=dFuelInner, d_pellet=dFuelPellet),  # Fuel rod
                  2: cdb.Tube(d_outer=dGuideTubeOuter, d_inner=dGuideTubeInner)}                 # Guide tube

    # Radial power shape
    radPow = [[0.00, 0.00, 0.00, 0.00, 0.75, 0.76, 0.77, 0.78,
               0.77, 0.76, 0.75, 0.00, 0.00, 0.00, 0.00],
              [0.00, 0.00, 0.00, 0.75, 1.10, 0.93, 1.10, 0.91,
               1.10, 0.93, 1.10, 0.75, 0.00, 0.00, 0.00],
              [0.00, 0.00, 0.74, 1.12, 0.93, 1.17, 0.97, 1.18,
               0.97, 1.17, 0.93, 1.12, 0.74, 0.00, 0.00],
              [0.00, 0.75, 1.12, 0.93, 1.12, 0.94, 1.19, 0.93,
               1.19, 0.94, 1.12, 0.93, 1.12, 0.75, 0.00],
              [0.75, 1.10, 0.93, 1.12, 0.95, 1.19, 0.95, 1.22,
               0.95, 1.19, 0.95, 1.12, 0.93, 1.10, 0.75],
              [0.76, 0.93, 1.17, 0.94, 1.19, 0.98, 1.19, 0.94,
               1.19, 0.98, 1.19, 0.94, 1.17, 0.93, 0.76],
              [0.77, 1.10, 0.97, 1.19, 0.95, 1.19, 0.94, 1.12,
               0.94, 1.19, 0.95, 1.19, 0.97, 1.10, 0.77],
              [0.78, 0.91, 1.18, 0.93, 1.22, 0.94, 1.12, 0.90,
               1.12, 0.94, 1.22, 0.93, 1.18, 0.91, 0.78],
              [0.77, 1.10, 0.97, 1.19, 0.95, 1.19, 0.94, 1.12,
               0.94, 1.19, 0.95, 1.19, 0.97, 1.10, 0.77],
              [0.76, 0.93, 1.17, 0.94, 1.19, 0.98, 1.19, 0.94,
               1.19, 0.98, 1.19, 0.94, 1.17, 0.93, 0.76],
              [0.75, 1.10, 0.93, 1.12, 0.95, 1.19, 0.95, 1.22,
               0.95, 1.19, 0.95, 1.12, 0.93, 1.10, 0.75],
              [0.00, 0.75, 1.12, 0.93, 1.12, 0.94, 1.19, 0.93,
               1.19, 0.94, 1.12, 0.93, 1.12, 0.75, 0.00],
              [0.00, 0.00, 0.74, 1.12, 0.93, 1.17, 0.97, 1.18,
               0.97, 1.17, 0.93, 1.12, 0.74, 0.00, 0.00],
              [0.00, 0.00, 0.00, 0.75, 1.10, 0.93, 1.10, 0.91,
               1.10, 0.93, 1.10, 0.75, 0.00, 0.00, 0.00],
              [0.00, 0.00, 0.00, 0.00, 0.75, 0.76, 0.77,
               0.78, 0.77, 0.76, 0.75, 0.00, 0.00, 0.00, 0.00]]
    # Assembly by assembly radial power shape
    radPow = [[0.00, 0.00, 0.00, 0.00, 0.75, 0.76, 0.77, 0.78, 0.77, 0.76, 0.75, 0.00, 0.00, 0.00, 0.00],
              [0.00, 0.00, 0.00, 0.75, 1.10, 0.93, 1.10, 0.91, 1.10, 0.93, 1.10, 0.75, 0.00, 0.00, 0.00],
              [0.00, 0.00, 0.74, 1.12, 0.93, 1.17, 0.97, 1.18, 0.97, 1.17, 0.93, 1.12, 0.74, 0.00, 0.00],
              [0.00, 0.75, 1.12, 0.93, 1.12, 0.94, 1.19, 0.93, 1.19, 0.94, 1.12, 0.93, 1.12, 0.75, 0.00],
              [0.75, 1.10, 0.93, 1.12, 0.95, 1.19, 0.95, 1.22, 0.95, 1.19, 0.95, 1.12, 0.93, 1.10, 0.75],
              [0.76, 0.93, 1.17, 0.94, 1.19, 0.98, 1.19, 0.94, 1.19, 0.98, 1.19, 0.94, 1.17, 0.93, 0.76],
              [0.77, 1.10, 0.97, 1.19, 0.95, 1.19, 0.94, 1.12, 0.94, 1.19, 0.95, 1.19, 0.97, 1.10, 0.77],
              [0.78, 0.91, 1.18, 0.93, 1.22, 0.94, 1.12, 0.90, 1.12, 0.94, 1.22, 0.93, 1.18, 0.91, 0.78],
              [0.77, 1.10, 0.97, 1.19, 0.95, 1.19, 0.94, 1.12, 0.94, 1.19, 0.95, 1.19, 0.97, 1.10, 0.77],
              [0.76, 0.93, 1.17, 0.94, 1.19, 0.98, 1.19, 0.94, 1.19, 0.98, 1.19, 0.94, 1.17, 0.93, 0.76],
              [0.75, 1.10, 0.93, 1.12, 0.95, 1.19, 0.95, 1.22, 0.95, 1.19, 0.95, 1.12, 0.93, 1.10, 0.75],
              [0.00, 0.75, 1.12, 0.93, 1.12, 0.94, 1.19, 0.93, 1.19, 0.94, 1.12, 0.93, 1.12, 0.75, 0.00],
              [0.00, 0.00, 0.74, 1.12, 0.93, 1.17, 0.97, 1.18, 0.97, 1.17, 0.93, 1.12, 0.74, 0.00, 0.00],
              [0.00, 0.00, 0.00, 0.75, 1.10, 0.93, 1.10, 0.91, 1.10, 0.93, 1.10, 0.75, 0.00, 0.00, 0.00],
              [0.00, 0.00, 0.00, 0.00, 0.75, 0.76, 0.77, 0.78, 0.77, 0.76, 0.75, 0.00, 0.00, 0.00, 0.00]]

    # Power as function of axial position
    axPow = lambda z: np.cos((z-coreHeight/2)*np.pi/2/coreHeight)

    builder = SquareLatticeLWR_Nodal(model, coreMap, pwr17x17, pitch, dz, solidGeoms, assemPitch=assemPitch,
                                     baffleGap=baffleGap, assemLosses=formLoss, assemLossLevels=gridLevels,
                                     rodDomains=rodDomains, modelGuideTubes=False, setCoreMaps=True)
                                     rodDomains=rodDomains, modelGuideTubes=False)

    # Run to 1 second to let the solution reach steady state
    model.addTransientGroup(tEnd=1.0)
    # model.addTransientGroup(tEnd=1.0)
    # Run a 10 second transient
    model.addTransientGroup(tEnd=11.0, editInterval=1.0)
    # model.addTransientGroup(tEnd=11.0, editInterval=1.0)

    # Ramp the power
    time = [1.0, 4.0]
    power = [1.0, 1.0]
    model.addTotalPowerRamp(time, power)
    #time = [1.0, 4.0]
    #power = [1.0, 1.0]
    #model.addTotalPowerRamp(time, power)

    # Set the core boundary conditions
    time = [1.0,  1.1]
    flow = [1.00, 1.00]
    builder.setCoreInletBC(inletMassFlux, inletTemp, flowRamp=[time, flow])
    #time = [1.0,  1.1]
    #flow = [1.00, 1.00]

    builder.setCoreInletBC(inletMassFlux, inletTemp)
    builder.setCoreOutletBC(outletPressure, inletTemp)

    # Set the nominal power and direct heating to coolant
    model.setAveragePower(qprime=linearHeatRate, dhfrac=directHeat)
    builder.setAveragePower(qprime=linearHeatRate, dhfrac=directHeat)

    # Turn on rod VTK and DNB edits.  They are both off by default.
    model.setEditOptions(dnbEdits=True,veracsHDF5=True)
    builder.setEditOptions(dnbEdits=True,veracsHDF5=True)

    coreHeight = np.sum(np.array(dz))

    # Power as function of axial position
    axPow = lambda z: np.cos((z-coreHeight/2)*np.pi/2/coreHeight)

    # set 3d power distribution
    # Set 3d power distribution
    if power_3d_flag == 1:
        radPow = np.asarray(radPow)
        axialMesh = builder.z_centers
@@ -288,7 +185,7 @@ def main(power_3d_flag=1):
        axPow = axPow(np.array(axialMesh))
        powerDist3d = pow3D * axPow
        builder.setCorePowerShape3D(powerDist3d)
        model.generateModel('asm3dpow_deck.inp')
        builder.generateModel('asm3dpow_deck.inp')
    elif power_3d_flag == 2:
        radPow = np.asarray(radPow)
        axialMesh = builder.z_centers
@@ -300,10 +197,10 @@ def main(power_3d_flag=1):
        axPow = axPow(np.array(axialMesh))
        powerDist3d = pow3D * axPow
        builder.setCorePowerShape3D(powerDist3d)
        model.generateModel('nodal3dpow_deck.inp')
        builder.generateModel('nodal3dpow_deck.inp')
    else:
        builder.setCorePowerShape(radialPower=radPow, axialPower=axPow)
        model.generateModel('orig_deck.inp')
        builder.generateModel('orig_deck.inp')


if __name__ == '__main__':