Commit 7827c571 authored by Gurecky, William's avatar Gurecky, William
Browse files

wip: adds 3d nodal power shape setter to square lattice lwr nodal builder

parent b0175602

SubKit/build/SquareLatticeLWR_Nodal.py

+90 −15
Original line number Diff line number Diff line
@@ -718,8 +718,10 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core): 
            if isinstance(solidGeoDict[solidTypeID], Solid.FuelRod) or modelGuideTubes:

               me.rodsInAssem[assemID].append(solidID)

               [i, j]=chMapObj.getLoc(chIdx)

               model.addSolid(solidID, solidTypeID, Solid.HeatedSolid(mult=mult/4.0, power=power, powID=powID,

                  x=nodalMesh.getXLoc(i, j), y=nodalMesh.getYLoc(i, j)))

               x, y = nodalMesh.getXLoc(i, j), nodalMesh.getYLoc(i, j)

               model.addSolid(solidID, solidTypeID,

                       Solid.HeatedSolid(mult=mult/4.0, power=power, powID=powID,

                       x=x, y=y))

               # Connect it ot the model

               model.addSolidOuterChan(solidID, chIdx, 1.0)
@@ -958,9 +960,82 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core): 
               powID = "pwr_" + str(assemID)

               me.model.addAxialPowerShape(powID, me.z_power_grid, powerDist[i, j, :])

               rods = me.rodsInAssem[assemID]



               for rodID in rods:

                  obj_rod = me.model.solidObjects[rodID]

                     objType_rod = me.model.solidTypeIDs[rodID]

                  obj_rod_type_id = me.model.solidTypeIDs[rodID]

                  # obj_rod_type_id == 1 is heated rod

                  # obj_rod_type_id == 2 is guide tube (only present if modelGuideTubes=True)

                  # obj_rod is type SubKit.build.Solid.HeatedSolid

                  obj_rod.powID = powID

                  # radial power factors are implicit in the 3d power distribution

                  obj_rod.power = 1.0



   def setCorePowerShapeNodal3D(me, powerDist, wgts=None, coreStart=0.0):

      """

      Set a 3D power distribution for the core.  Each assembly can have it's own

      axial power shape.



      Args:

         powerDist (numpy.ndarray): 4 dim array. Contains power distribution with indexing

           (asm_i, asm_j, node_id, asm_axial)

         wgts (numpy.ndarray): Same shape as powerDist. Weights of powers (optional)

         coreStart (float): The axial location of the start of the core.  In CTF, the bottom of the

            model is zero.  If the core does not start at the bottom of the model, this must be provided

            so the power tables are correctly defined in CTF. [m]

      """

      assert len(powerDist.shape) == 4

      assert len(me.z_power_grid.shape) == 1

      # Nodal mesh is 4x larger than core map dims

      assert powerDist.shape[0] == me.chMapObj.getAssemDimLen()

      assert powerDist.shape[1] == me.chMapObj.getAssemDimLen()

      assert powerDist.shape[2] == 4

      assert me.z_power_grid.size == powerDist.shape[2]



      # normalize 3d powers

      if wgts is None:

          wgts = np.ones(powerDist.shape)

      assert wgts.shape == powerDist.shape

      powerDist *= wgts

      powerDist /= np.average(powerDist, weights=wgts)



      radPowers = []

      weights = []

      for obj in me.model.solidObjects.values():

         assert(obj.symtype==1)

      for asm_i in range(me.chMapObj.getAssemDimLen()):

         for asm_j in range(me.chMapObj.getAssemDimLen()):

            assemID = me.chMapObj.getAssemIndexFromLoc(asm_i, asm_j)

            if assemID:

               rods = me.rodsInAssem[assemID]



               # determine quadrent of rod in assem

               rod_xy = {}

               for rodID in rods:

                  obj_rod = me.model.solidObjects[rodID]

                  rod_xy[rodID] = (obj_rod.x, obj_rod.y)

               unique_xy = np.unique(np.asarray(rod_xy.values()), axis=0)

               assert unique_xy.size == 4 * 2

               centroid_xy = np.average(unique_xy, axis=0)

               for rodID in rods:

                  obj_rod = me.model.solidObjects[rodID]

                  obj_rod_type_id = me.model.solidTypeIDs[rodID]

                  if obj_rod.x > centroid_xy[0] and obj_rod.y > centroid_xy[1]:

                      # north east

                      nodeID = 0

                  elif obj_rod.x < centroid_xy[0] and obj_rod.y > centroid_xy[1]:

                      # north west

                      nodeID = 1

                  elif obj_rod.x < centroid_xy[0] and obj_rod.y < centroid_xy[1]:

                      # south west

                      nodeID = 2

                  else:

                      # south east

                      nodeID = 3



                  # set rod power

                  powID = "pwr_asm_" + str(assemID) + "_node_" + str(nodeID)

                  me.model.addAxialPowerShape(powID, me.z_power_grid,

                          powerDist[asm_i, asm_j, nodeID, :])

                  obj_rod.powID = powID

                  obj_rod.power = 1.0