Commit 5e433f1a authored by Phathanapirom U B's avatar Phathanapirom U B
Browse files

freeze after fetch willfork

parent 420dc34f
Loading
Loading
Loading
Loading
+297 −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 numpy as np
import math
import copy
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
    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
    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]
    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,veracsHDF5=True)

    coreHeight = np.sum(np.array(dz))

    # Power as function of axial position
    axPow = lambda z: np.cos((z-coreHeight/2)*np.pi/2/coreHeight)

    # set 3d power distribution
    if power_3d_flag == 1:
        radPow = np.asarray(radPow)
        axialMesh = builder.z_centers
        pow3D = np.ones(
            (radPow.shape[0], radPow.shape[1], len(axialMesh)))
        for z_i in range(len(axialMesh)):
            pow3D[:, :, z_i] = pow3D[:, :, z_i] * radPow
        axPow = axPow(np.array(axialMesh))
        powerDist3d = pow3D * axPow
        builder.setCorePowerShape3D(powerDist3d)
        model.generateModel('asm3dpow_deck.inp')
    elif power_3d_flag == 2:
        radPow = np.asarray(radPow)
        axialMesh = builder.z_centers
        pow3D = np.ones(
            (radPow.shape[0], radPow.shape[1], 4, len(axialMesh)))
        for z_i in range(len(axialMesh)):
            for node_j in range(4):
                pow3D[:, :, node_j, z_i] = pow3D[:, :, node_j, z_i] * radPow
        axPow = axPow(np.array(axialMesh))
        powerDist3d = pow3D * axPow
        builder.setCorePowerShape3D(powerDist3d)
        model.generateModel('nodal3dpow_deck.inp')
    else:
        builder.setCorePowerShape(radialPower=radPow, axialPower=axPow)
        model.generateModel('orig_deck.inp')


if __name__ == '__main__':
    # build 3 decks which should produce nearly the same CTF TH solution output
    main(power_3d_flag=0)
    main(power_3d_flag=1)
    main(power_3d_flag=2)