updated nodal power calc

## ------------------------------------------------------ ##

# Author: Birdy

# Date: 1/10/20

# Description:

# Converts pin-by-pin powers into nodal powers

import h5py

import numpy as np

import sys

def write_nodal_powers(file_name):

    """

    Write node-averaged power to an existing H5 file. A node is a quadrant of an assembly.

    Args:

        file_name (str):        Path to H5 file

    Assumes assemblies are symmetric across x and y dimensions (e.g., 17x17 assembly)

    """

    with h5py.File(file_name, 'r+') as data:

        nStates = get_nStates(data)

        for istate in range(nStates):

            core = data['CORE']

            state = data['STATE_{:04d}'.format(istate+1)]

            if 'pin_powers_nodal' in state.keys():

                del state['pin_powers_nodal']

            inserted_pin_powers = insert_pin_powers(state)

            mult_mask = get_mult_mask(state)

            weighted_pin_powers = get_weighted_powers(inserted_pin_powers,mult_mask)

            node_mask = get_nodal_mask(state)

            nodal_power = get_nodal_powers(weighted_pin_powers,node_mask,mult_mask)

import os

            state['pin_powers_nodal'] = nodal_power

class Reactor:

    data.close()

    def __init__(self,core,state):

        self.core = core

        self.state = state

def get_nStates(data):

        if 'pin_powers_nodal' in self.state.keys():

            del self.state['pin_powers_nodal']

    nStates = 0

    for key in data.keys():

        if 'STATE' in key: nStates+=1

    return(nStates)

def get_assembly_dims(state):

    pin_powers = np.squeeze(state['pin_powers'][:])

    n_asm_pin_x = pin_powers.shape[0]

    n_asm_pin_y = pin_powers.shape[1]

    n_asm_z = pin_powers.shape[2]

    return(n_asm_pin_x,n_asm_pin_y,n_asm_z)

def get_assembly_midpoint(state):

    n_x,n_y,_ = get_assembly_dims(state)

    if n_x % 2 == 0:

        mid_x, mid_y = n_x/2, n_y/2

        odd = False

        self.n_x, self.n_y, self.n_z, self.n_asm = self.state['pin_powers'][()].shape

        if self.n_x % 2 == 0:

            self.mid_x, self.mid_y = self.n_x/2, self.n_y/2

            self.odd = False

        else:

        mid_x, mid_y = (n_x-1)/2, (n_y-1)/2 

        odd = True

    return(mid_x,mid_y,odd)

            self.mid_x, self.mid_y = (self.n_x-1)/2, (self.n_y-1)/2

            self.odd = True

        self.pin_powers = state['pin_powers'][()]

        self.__get_inserted_pin_powers__()

        self.__get_mult_mask__()

        self.__get_nodal_mask__()

        self.__get_nodal_powers__()

        self.__perform_exit_checks__()

def get_nodal_mask(state):

    # (inline) Duplicate pin powers at the node boundaries

    def __get_inserted_pin_powers__(self):

    """

    Map of nodal regions for weighted pin power mapping

    Insert duplicated pin powers along node edges

    """

        if self.odd:

    mid_x,mid_y,odd = get_assembly_midpoint(state)

            self.inserted_pin_powers = np.zeros((self.n_x+1,self.n_y+1,self.n_z,self.n_asm))

            for asm in range(0, self.n_asm):

                for z in range(0, self.n_z):

    if odd:

        mask_dims = (mid_x+1,mid_y+1)

    else:

        mask_dims = (mid_x,mid_y)

    mask_0 = np.zeros(mask_dims, dtype=int)

    mask_1 = np.ones((mask_0.shape),dtype=int)

    mask_2 = np.full((mask_0.shape),2,dtype=int)

    mask_3 = np.full((mask_0.shape),3,dtype=int)

    mask = np.vstack(( np.hstack((mask_0,mask_1)),np.hstack((mask_2,mask_3)) ))

                    dummy = self.pin_powers[:,:,z,asm]

                    ins_row = dummy[self.mid_x,:]

                    dummy = np.insert(dummy,self.mid_x,ins_row,axis=0)

                    ins_col = dummy[:,self.mid_y]

                    dummy = np.insert(dummy,self.mid_y,ins_col,axis=1)

                    self.inserted_pin_powers[:,:,z,asm] = dummy

    return(mask)

        else: self.inserted_pin_powers = self.pin_powers

def insert_pin_powers(state):

    def __get_mult_mask__(self):

    """

    Duplicate pin powers at the node boundaries

        via inline insertion of centerline row and column

    Passed with mult_mask to get weighted pin powers across the assembly

    Define the pin power weights across the duplicated pin powers

    for a single axial slice

    Weights at node boundaries are 0.5 and at the vertex 0.25

    Regions are symmetric about x,y plane

    """

    n_x,n_y,n_z = get_assembly_dims(state)

    mid_x,mid_y,odd = get_assembly_midpoint(state)

    pin_powers = np.squeeze(state['pin_powers'][:]) 

    if odd:

        inserted_pin_powers = np.zeros((n_x+1,n_y+1,n_z))

        for z in range(0,n_z):

            dummy = pin_powers[:,:,z]

            ins_row = dummy[mid_x,:]

            dummy = np.insert(dummy,mid_x,ins_row,axis=0)

            ins_col = dummy[:,mid_y]

            dummy = np.insert(dummy,mid_y,ins_col,axis=1)

            inserted_pin_powers[:,:,z] = dummy

    else: inserted_pin_powers = pin_powers

    return(inserted_pin_powers)

def get_mult_mask(state):

    """

    Map of pin power weights across the inserted_pin_powers

        One assembly is broken into regions ((a,b),(c,d)) = (L,R)

    Weights at node boundaries are 0.5

    Weights at node vertex are 0.25

    """

    n_x,n_y,n_z  = get_assembly_dims(state)

    mid_x,mid_y,odd = get_assembly_midpoint(state)

    # core = ((a,b),(c,d)) = ((L,R))

    # pin powers at node edges weighted at 0.5

    # pin powers at node vertices weighted at 0.25

        if self.odd:

    if odd:

        mask_d = np.ones((mid_x+1,mid_y+1),dtype=float)

            mask_d = np.ones((self.mid_x+1,self.mid_y+1),dtype=float)

            mask_d[0,:],mask_d[:,0] = 0.5,0.5

            mask_d[0,0] = 0.25

            R = np.concatenate((mask_b,mask_d),axis=0)

            L = np.fliplr(R[:,:])

        mask = np.concatenate((L,R),axis=1)

    else:

        mask = np.ones((n_x,n_y))

    return(mask)

            self.mult_mask = np.concatenate((L,R),axis=1)

        else: self.mult_mask = np.ones((self.n_x,self.n_y))

def get_weighted_powers(pin_powers,mask):

    def __get_nodal_mask__(self):

    """

    Get weighted pin powers across the assembly

        pin_powers have been duplicated at node boundaries 

    Define the nodal mask

    Axial slices are mapped to ((0,1),(2,3)) = (T,B)

    """

    weighted_powers = np.zeros((pin_powers.shape))

    for z in range(0,pin_powers.shape[2]):

        weighted_powers[:,:,z] = np.multiply(pin_powers[:,:,z],mask)

    return(weighted_powers)

        if self.odd: mask_dims = (self.mid_x+1,self.mid_y+1)

        else: mask_dims = (self.mid_x,self.mid_y)

        mask_0 = np.zeros(mask_dims,dtype=int)

        mask_1 = np.ones(mask_dims,dtype=int)

        mask_2 = np.full(mask_dims,2,dtype=int)

        mask_3 = np.full(mask_dims,3,dtype=int)

        self.node_mask = np.vstack(( np.hstack((mask_0,mask_1)),np.hstack((mask_2,mask_3)) ))

def get_nodal_powers(pin_powers,node_mask,mult_mask):

    def __get_nodal_powers__(self):

    """

    **** Apply mult_mask to inserted_pin_powers ****

    renormalize power

    Multiply element-by-element inserted pin powers; pin power weights

    Sum and group powers using nodal mask

    """

    n_z = pin_powers.shape[2]

    nodal_power = np.zeros((2,2,n_z))

    for z in range(0,pin_powers.shape[2]):

        power, weight, weighted_power = np.zeros((4,1)), np.zeros((4,1)), np.zeros((4,1))

        self.nodal_powers = np.zeros((2,2,self.n_z,self.n_asm))

        weighted_pin_powers = np.zeros(self.inserted_pin_powers.shape)

        power = np.zeros((4,1))

        for asm in range(0,self.n_asm):

            for z in range(0,self.n_z):

                weighted_pin_powers = np.multiply(self.inserted_pin_powers[:,:,z,asm],self.mult_mask)

                for n in range(0,4):

                    mask = self.node_mask == n

                    power[n] = np.sum( np.multiply(weighted_pin_powers,mask) )

                self.nodal_powers[:,:,z,asm] = power.reshape((2,2))

            for i in range(pin_powers.shape[0]):

                for j in range(pin_powers.shape[1]):

                    if node_mask[i,j] == n:

                        power[n] += pin_powers[i,j,z]

                        weight[n] += mult_mask[i,j] 

        weighted_power = np.divide(power,weight)

        nodal_power[:,:,z] = weighted_power.reshape(2,2) 

    def __perform_exit_checks__(self):

    """

    Check estimated nodal powers result in consistent:

    1. total power

       across entire core

    """

    return(nodal_power)

        total_power = self.core['rated_power'] * (self.state['power'][()] / 100.)

        core_avg_linear_heatrate = self.state['core_avg_linear_heatrate'][()]

        core_map = self.core['core_map'][()].flatten()-1

        axial_mesh = self.core['axial_mesh'][:]

        mesh = axial_mesh[1:] - axial_mesh[:-1]

        ax_length = sum(mesh)

def main():

        local_power = np.zeros(self.nodal_powers.shape)

        for asm in range(0,self.n_asm):

            for z in range(0,self.n_z):

                local_power[:,:,z,asm] = self.nodal_powers[:,:,z,asm] \

                    * core_avg_linear_heatrate * mesh[z]

    file_name = sys.argv[1]

    write_nodal_powers(file_name)

        # total power

        est_total_power = 0.

        for core_asm in range(0,len(core_map)):

           est_total_power += np.sum(local_power[:,:,:,core_map[core_asm]])

        est_total_power = est_total_power * 1e-6 # convert W to MW

if __name__ == '__main__':

    main()

        if not np.isclose(est_total_power, total_power, rtol=1e-3):

            print('rated power out of tolerance 1e-3')