Commit ec3b7b8d authored by Gurecky, William's avatar Gurecky, William Committed by Salko Jr, Robert
Browse files

Adds alternate interpolation method to get nodal powers onto the ctf grid

parent 79d4d5ae
Loading
Loading
Loading
Loading
+48 −16
Original line number Diff line number Diff line
@@ -990,7 +990,7 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core):
      for i, obj in enumerate(me.model.solidObjects.values()):
         obj.power = radP[i]

   def setCorePowerShape3D(me, powerDist, wgts=None, coreStart=0.0):
   def setCorePowerShape3D(me, powerDist, wgts=None, coreStart=0.0, interp_type="nearest"):
      """
      Set a 3D power distribution for the core.
      Args:
@@ -1013,11 +1013,24 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core):
         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]
         interp_type (str): One of "nearest" or "linear".
            case nearest:
               a     bc   de         fg
               |--p--||-p-||----p----||--...
                  1     2       3
               Where the user specifies nodal power at p locations
               and a...g represent upper and lower bounds of each
               axial zone as indicated in z_power_grid. In this case:
               The power given at p1 is copied to 'a' and 'b'.
               The power at p2 is copied to 'c' and 'd'.
            case linear:
               Linear interpolation is used to obtain powers at the
               axial zone boundaries.
      """
      if len(powerDist.shape) == 3:
         me._setCorePowerShapeAsm3D(powerDist, wgts, coreStart)
         me._setCorePowerShapeAsm3D(powerDist, wgts, coreStart, interp_type)
      elif len(powerDist.shape) == 4:
         me._setCorePowerShapeNodal3D(powerDist, wgts, coreStart)
         me._setCorePowerShapeNodal3D(powerDist, wgts, coreStart, interp_type)
      else:
         raise RuntimeError("Invalid power distribution shape: %s" % str(powerDist.shape))

@@ -1031,7 +1044,7 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core):
      for i, obj in enumerate(me.model.solidObjects.values()):
         obj.power = radP[i]

   def _setCorePowerShapeAsm3D(me, powerDist, wgts=None, coreStart=0.0):
   def _setCorePowerShapeAsm3D(me, powerDist, wgts=None, coreStart=0.0, interp_type="nearest"):
      """
      Set a 3D power distribution for the core.  Each assembly can have its own
      axial power shape.
@@ -1043,6 +1056,7 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core):
         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]
         interp_type (str): One of "nearest" or "linear"
      """
      assert len(powerDist.shape) == 3
      assert len(me.z_power_grid.shape) == 1
@@ -1068,17 +1082,25 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core):

               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

               if interp_type == "linear":
                  # lin interpolate to the z_bounds
                  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)
               else:
                  asm_axial_power_ctf = np.zeros(me.z_power_grid.size)
                  for z_lvl in range(me.z_centers.size):
                     asm_axial_power_ctf[[2*z_lvl, 2*z_lvl+1]] = \
                             asm_axial_power_profile[z_lvl]

               # Normalize such that axially integrated asm nodal power is 1
               avg_axpwr = np.max((np.trapz(asm_axial_power_ctf, x=me.z_power_grid) / \
                     (np.max(me.z_power_grid) - np.min(me.z_power_grid)), 1e-12))

               # set axial power profile in the ctf model
               me.model.addAxialPowerShape(powID, me.z_power_grid, asm_axial_power_ctf / avg_axpwr)
               rods = me.rodsInAssem[assemID]
               for rodID in rods:
@@ -1091,7 +1113,7 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core):
                  # radial power factors are implicit in the 3d power distribution
                  obj_rod.power = avg_axpwr

   def _setCorePowerShapeNodal3D(me, powerDist, wgts=None, coreStart=0.0):
   def _setCorePowerShapeNodal3D(me, powerDist, wgts=None, coreStart=0.0, interp_type="nearest"):
      """
      Set a 3D power distribution for the core.  Each node within each assembly
      can have its own axial power shape.
@@ -1111,6 +1133,7 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core):
         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]
         interp_type (str): One of "nearest" or "linear"
      """
      assert len(powerDist.shape) == 4
      assert len(me.z_power_grid.shape) == 1
@@ -1169,12 +1192,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

                  if interp_type == "linear":
                     # lin interpolate to the z_bounds
                     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)
                  else:
                     # alternative interp to ctf grid
                     node_axial_power_ctf = np.zeros(me.z_power_grid.size)
                     for z_lvl in range(me.z_centers.size):
                        node_axial_power_ctf[[2*z_lvl, 2*z_lvl+1]] = \
                                node_axial_power_profile[z_lvl]

                  # Normalize such that axially integrated asm nodal power is 1
                  avg_axpwr = np.max((np.trapz(node_axial_power_ctf, x=me.z_power_grid) / \
                        (np.max(me.z_power_grid) - np.min(me.z_power_grid)), 1e-12))

+72520 −72520

File changed.

Preview size limit exceeded, changes collapsed.