Commit a42adac9 authored by Gurecky, William's avatar Gurecky, William
Browse files

wip: fix power norm

parent 57bc3559

SubKit/build/SquareLatticeLWR_Nodal.py

+20 −4
Original line number Diff line number Diff line
@@ -1026,7 +1026,7 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core): 
          wgts = np.ones(powerDist.shape)

      assert wgts.shape == powerDist.shape

      powerDist *= wgts

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

      # powerDist /= np.average(powerDist[np.nonzero(powerDist)], weights=wgts)



      radPowers = []

      weights = []
@@ -1084,7 +1084,7 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core): 
          wgts = np.ones(powerDist.shape)

      assert wgts.shape == powerDist.shape

      powerDist *= wgts

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

      # powerDist /= np.average(powerDist[np.nonzero(powerDist)], weights=wgts)



      radPowers = []

      weights = []
@@ -1125,12 +1125,28 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core): 


                  # set rod power

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

                  node_axial_power_profile = powerDist[asm_i, asm_j, nodeID, :]

                  avg_node_axial_power_profile = \

                      np.max((1e-12, np.average(node_axial_power_profile)))

                  # scale axial power

                  # import pdb; pdb.set_trace()

                  # node_axial_power_profile *= 1.165803

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

                          powerDist[asm_i, asm_j, nodeID, :])

                          node_axial_power_profile / avg_node_axial_power_profile)

                  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_node_axial_power_profile



      # Renormalize the radial powers

      radP = []

      mults = []

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

         radP.append(obj.power)

         mults.append(obj.mult)

      radP = np.array(radP)/np.average(np.array(radP), weights=mults)

      for i, obj in enumerate(me.model.solidObjects.values()):

         obj.power = radP[i]

tests/nodal_p7_pow3d/buildDeck.py

+5 −19
Original line number Diff line number Diff line 
#!/usr/bin/env python

# Author: Bob Salko

# Date: 1/23/19

# Author: Bob Salko, W. Gurecky

# Date: 1/23/19, 1/29/20

# Description:

# Makes a model of CASL problem 7 using 4 subchannels per assembly

#

import numpy as np

import math

import copy

import sys

sys.path.insert(0, '../../')

from SubKit.build.SquareLatticeLWR_Nodal import SquareLatticeLWR_Nodal

import SubKit.utils.UnitConversions as units

import SubKit.build as cdb

import sys

sys.path.insert(0, '../../')





def main(power_3d_flag=1):



    dFuelOuter = 0.475e-2*2  # m

    dFuelInner = 0.418e-2*2  # m

    dFuelPellet = 0.4096e-2*2  # m
@@ -149,20 +148,7 @@ def main(power_3d_flag=1): 
                                     baffleGap=baffleGap, assemLosses=formLoss, assemLossLevels=gridLevels,

                                     rodDomains=rodDomains, modelGuideTubes=False)



    # Run to 1 second to let the solution reach steady state

    # model.addTransientGroup(tEnd=1.0)

    # Run a 10 second transient

    # 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)



    # Set the core boundary conditions

    #time = [1.0,  1.1]

    #flow = [1.00, 1.00]



    # Set boundary conditions

    builder.setCoreInletBC(inletMassFlux, inletTemp)

    builder.setCoreOutletBC(outletPressure, inletTemp)