Commit 2559f6ef authored by Gurecky, William's avatar Gurecky, William
Browse files

fix for axial power norm in setCorePowerShapeNodal3D

parent 3b7877dd

SubKit/build/CoreBuilder.py

+4 −6
Original line number Diff line number Diff line
@@ -131,15 +131,13 @@ class Core(object): 
   @property

   def z_power_grid(me):

      """ Adjust z center endpoints to span entire axial core height """

      #z_pwr_grid = me.z_centers

      #z_pwr_grid[0] = 0.0

      #z_pwr_grid[-1] = np.cumsum(me.dz)[-1] + 0.0001



      z_pwr_grid = me.z_bounds_pairs

      # slightly adjust the upper bounds of each axial edit zone

      z_pwr_grid[:, 1] -= 1e-8

      # flatten and return

      return z_pwr_grid.flatten()

      z_pwr_grid = z_pwr_grid.flatten()

      z_pwr_grid[-1] = np.cumsum(me.dz)[-1] + 1e-8

      return z_pwr_grid



   @property

   def z_bounds(me):
@@ -157,7 +155,7 @@ class Core(object): 
        [z_1, z_2],  # axial edit i=1

        [z_2, z_3], ...]

      """

      z_pairs = np.asarray((me.z_bounds[1:], me.z_bounds[:-1])).T

      z_pairs = np.asarray((me.z_bounds[:-1], me.z_bounds[1:])).T

      return z_pairs



   # === Delegate to me.model ===

SubKit/build/SquareLatticeLWR_Nodal.py

+25 −11
Original line number Diff line number Diff line
@@ -1047,7 +1047,7 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core): 
      assert len(me.z_power_grid.shape) == 1

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

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

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

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



      # normalize 3d powers

      if wgts is None:
@@ -1064,12 +1064,22 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core): 
            assemID = me.chMapObj.getAssemIndexFromLoc(i, j)

            if assemID:

               powID = "pwr_" + str(assemID)

               asm_axial_power_profile_avg = \

                   np.max((1e-12, np.average(powerDist[i, j, :])))

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

                       asm_axial_power_profile_avg)

               rods = me.rodsInAssem[assemID]



               asm_axial_power_profile = powerDist[i, j, :]



               # the asm_axial_power_profile is specified at z_centers

               assert me.z_centers.size == asm_axial_power_profile.size

               asm_axial_power_interpolant = \

                       interp1d(me.z_centers, asm_axial_power_profile,

                                kind='linear', fill_value='extrapolate')



               # interpolate to the z_bounds

               asm_axial_power_ctf = asm_axial_power_interpolant(me.z_power_grid)

               avg_axpwr = np.trapz(asm_axial_power_ctf, x=me.z_power_grid) / \

                     (np.max(me.z_power_grid) - np.min(me.z_power_grid))



               me.model.addAxialPowerShape(powID, me.z_power_grid, asm_axial_power_ctf / avg_axpwr)

               rods = me.rodsInAssem[assemID]

               for rodID in rods:

                  obj_rod = me.model.solidObjects[rodID]

                  obj_rod.powID = powID
@@ -1078,7 +1088,7 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core): 
                  # obj_rod_type_id == 2 is guide tube (only present if me.modelGuideTubes=True)

                  # obj_rod is type SubKit.build.Solid.HeatedSolid

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

                  obj_rod.power = asm_axial_power_profile_avg

                  obj_rod.power = avg_axpwr



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

      """
@@ -1107,7 +1117,7 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core): 
      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[3]

      assert me.z_centers.size == powerDist.shape[3]



      # normalize 3d powers

      if wgts is None:
@@ -1158,17 +1168,21 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core): 


                  # the node_axial_power_profile is specified at z_centers

                  assert me.z_centers.size == node_axial_power_profile.size

                  node_axial_power_interpolant = interp1d(me.z_centers, node_axial_power_profile)

                  node_axial_power_interpolant = \

                          interp1d(me.z_centers, node_axial_power_profile,

                                   kind='linear', fill_value='extrapolate')



                  # interpolate to the z_bounds

                  node_axial_power_ctf = node_axial_power_interpolant(me.z_power_grid)

                  avg_axpwr = np.trapz(node_axial_power_ctf, x=me.z_power_grid) / \

                        (np.max(me.z_power_grid) - np.min(me.z_power_grid))



                  # set axial power profile in the ctf model

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

                  me.model.addAxialPowerShape(powID, me.z_power_grid, node_axial_power_ctf / avg_axpwr)

                  obj_rod.powID = powID

                  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 me.modelGuideTubes=True)

                  # obj_rod is type SubKit.build.Solid.HeatedSolid

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

                  obj_rod.power = 1.0
 
                  obj_rod.power = avg_axpwr