Commit 0fb7a3a2 authored by Gurecky, William's avatar Gurecky, William
Browse files

adding method to LWR square lattice deck builder to set 3d power distribution

parent 52bd3eb8
Loading
Loading
Loading
Loading
+3 −0
Original line number Diff line number Diff line
@@ -83,6 +83,9 @@ class FuelRod(SolidGeo):
      me.num_rings = num_rings
      me.percentTheoreticalDensity = percentTheoreticalDensity

class HeatedFuelRod(FuelRod):
   def __init__(me, d_outer, d_inner, d_pellet, pellet_material=None, clad_material=None, num_rings=10, percentTheoreticalDensity=1.0, powID=None):
      super(HeatedFuelRod, me).__init__(d_outer, d_inner, d_pellet, pellet_material, clad_material, num_rings, percentTheoreticalDensity)

class _Material:
   """ Defines a solid material.
+41 −0
Original line number Diff line number Diff line
@@ -893,3 +893,44 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core):
      for i, obj in enumerate(me.model.solidObjects.values()):
         obj.power = radP[i]

   @property
   def z_centers(me):
      z_grid = np.cumsum(me.dz) - me.dz[0]
      z_grid[-1] = np.cumsum(me.dz)[-1] + 0.0001
      return z_grid

   def setCorePowerShape3D(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): 3 dim array. Contains power distribution with indixing
           (asm_i, asm_j, asm_axial)
      """
      assert len(powerDist.shape) == 3
      assert len(me.z_centers.shape) == 1
      assert powerDist.shape[0] == me.chMapObj.getAssemDimLen()
      assert powerDist.shape[1] == me.chMapObj.getAssemDimLen()
      assert me.z_centers.size == powerDist.shape[2]

      # expand powerDist to NodalMesh

      # normalize 3d powers
      if wgts is None:
          wgts = np.ones(powerDist.shape)
      powerDist /= np.average(powerDist, weights=wgts)

      radPowers = []
      weights = []
      for obj in me.model.solidObjects.values():
         assert(obj.symtype==1)
         for i in range(me.chMapObj.getAssemDimLen()):
            for j in range(me.chMapObj.getAssemDimLen()):
               assemID = me.chMapObj.getAssemIndexFromLoc(i, j)
               if assemID:
                  rods = me.rodsInAssem[assemID]
                  for rod in rods:
                     powID = "pwr_" + str(rod)
                     me.model.addAxialPowerShape(powID, me.z_centers, powerDist[i, j, :])
                     obj.powID = powID
+261 −0
Original line number Diff line number Diff line
#!/usr/bin/env python
# Author: Bob Salko
# Date: 1/23/19
# Description:
# Makes a model of CASL problem 7 using 4 subchannels per assembly
#
import sys
sys.path.insert(0, '../../')
import SubKit.build as cdb
import SubKit.utils.UnitConversions as units
from   SubKit.build.SquareLatticeLWR_Nodal import SquareLatticeLWR_Nodal
import copy
import math
import numpy as np


def main():

   dFuelOuter = 0.475e-2*2 # m
   dFuelInner = 0.418e-2*2 # m
   dFuelPellet = 0.4096e-2*2 # m
   dGuideTubeOuter = 0.602e-2*2 # m
   dGuideTubeInner = 0.561e-2*2 # m
   assemPitch = 21.5e-2 # m
   pitch = 1.26e-2 # m
   baffleGap = 0.19e-2 # m
   inletMassFlux = 0.3483429E+04 # kg/m**2/s
   inletTemp = 291.85 # C
   outletPressure = 155.132039 # bar
   linearHeatRate = 16.71977 # kW/m
   directHeat = 0.026
   gridLocations = list(np.array([13.884,  75.2, 127.4, 179.6, 231.8, 284.0, 336.2])*units.t_cm_m)
   gridHeights =   list(np.array([3.866 , 3.810, 3.810, 3.810, 3.810, 3.810, 3.810])*units.t_cm_m)
   axial_edit_bounds = list(np.array([11.951, 15.817, 24.028, 32.239, 40.45, 48.662, 56.873, 65.084, 73.295, 77.105 , 85.17,
      93.235, 101.3, 109.365, 117.43, 125.495, 129.305, 137.37, 145.435, 153.5, 161.565, 169.63, 177.695,
      181.505, 189.57, 197.635, 205.7, 213.765, 221.83, 229.895, 233.705, 241.77, 249.835, 257.9, 265.965, 274.03,
      282.095, 285.905, 293.97, 302.035, 310.1, 318.165, 326.23, 334.295, 338.105, 346.0262, 353.9474, 361.8686,
      369.7898, 377.711])*units.t_cm_m)
   formLoss = 0.9070

   # Setup axial mesh
   gridLocations = list(np.array(gridLocations)-axial_edit_bounds[0])
   axial_edit_bounds = list(np.array(axial_edit_bounds)-axial_edit_bounds[0])
   dz = []
   for j in range(1, len(axial_edit_bounds)):
      dz.append(axial_edit_bounds[j]-axial_edit_bounds[j-1])
   numLevels = len(dz)
   gridBottoms = list(np.array(gridLocations)-np.array(gridHeights)/2.0)
   gridLevels = []
   for grid in gridBottoms:
      found = False
      for level, z in enumerate(axial_edit_bounds):
         if np.isclose(grid, z):
            gridLevels.append(level+1)
            found = True
      assert(found)
   assert(len(gridLevels)==len(gridLocations))


   # The map of the core.  Each 1 represents an assembly
   coreMap = [[ 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0],
              [ 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0],
              [ 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0],
              [ 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0],
              [ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
              [ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
              [ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
              [ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
              [ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
              [ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
              [ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
              [ 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0],
              [ 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0],
              [ 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0],
              [ 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0]]

   # Fuel rod is type 1, but mark it with an "f" to make it more visible
   # Guide tube is type 2
   f=1
   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, 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, 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, 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, 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, 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]]
   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],
                 [ 0, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 0],
                 [ 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3],
                 [ 1, 1, 1, 1, 7, 7, 7, 2, 2, 2, 3, 3, 3, 3, 3],
                 [ 1, 1, 1, 1, 7, 7, 7, 7, 7, 7, 3, 3, 3, 3, 3],
                 [ 1, 1, 1, 1, 7, 7, 7, 7, 7, 7, 7, 6, 6, 6, 6],
                 [ 4, 4, 4, 4, 4, 7, 7, 7, 7, 7, 7, 6, 6, 6, 6],
                 [ 4, 4, 4, 4, 4, 7, 7, 7, 7, 7, 7, 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]]
   #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]]

   print np.unique(np.array(rodDomains), return_counts=True)

   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]]

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

   # 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 down over time
   time = [1.0, 4.0, 4.3, 4.6, 4.9, 5.2, 5.5, 5.8, 6.1, 6.4, 6.7, 7.0, 11.0]
   power = [1.0, 1.0, 0.99, 0.95, 0.90, 0.80, 0.68, 0.44, 0.27, 0.21, 0.19, 0.18, 0.16]
   model.addTotalPowerRamp(time, power)

   # Set the core boundary conditions
   # Apply a flow ramp to represent the loss of flow
   time = [ 1.0,  1.1,  1.2,  1.4,  1.6,  1.9,  2.2,  2.5,  2.8,  3.1,  3.4,  3.8,  4.0,  4.6,  5.4,  6.5,  7.9,  9.2, 10.0]
   flow = [1.00, 1.00, 0.90, 0.83, 0.78, 0.75, 0.71, 0.70, 0.67, 0.66, 0.64, 0.62, 0.61, 0.58, 0.55, 0.51, 0.46, 0.43, 0.40]
   builder.setCoreInletBC(inletMassFlux, inletTemp, flowRamp=[time, flow])
   builder.setCoreOutletBC(outletPressure, inletTemp)

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

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

   coreHeight = np.sum(np.array(dz))
   def axPow(z):
      return np.cos((z-coreHeight/2)*np.pi/2/coreHeight)

   # set 3d power distribution
   axialMesh = builder.z_centers
   radPow_arr = np.array(radPow)
   radPow_arr_full = np.ones((radPow_arr.shape[0], radPow_arr.shape[1], len(axialMesh)))
   for z_i in range(len(axialMesh)):
       radPow_arr_full[:, :, z_i] = radPow_arr_full[:, :, z_i] * radPow_arr
   axPow_arr = axPow(np.array(axialMesh))
   powerDist3d = radPow_arr_full * axPow_arr
   builder.setCorePowerShape3D(powerDist3d)

   model.generateModel()

if __name__=='__main__':
   main()