adding assembly power calculation

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

# Author: Birdy

# Date: 1/10/20

# Description:

# Converts pin-by-pin powers into nodal powers

import h5py

import numpy as np

import os

class Reactor:

    def __init__(self,core,state):

        self.core = core

        self.state = state

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

            del self.state['pin_powers_assembly']

        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:

            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_assembly_powers__()

        self.__perform_exit_checks__()

        self.state.create_dataset('pin_powers_assembly', data=self.assembly_powers)

    def __get_assembly_powers__(self):

        """ Sum pin powers across assembly """

        self.assembly_powers = np.zeros((self.n_z,self.n_asm))

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

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

                self.assembly_powers[z,asm] = np.sum( self.pin_powers[:,:,z,asm] )

    def __perform_exit_checks__(self):

        """

        Check estimated nodal powers result in consistent:

        1. total power

           across entire core

        """

        total_power = self.core['rated_power'] * np.divide(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)

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

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

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

                local_power[z,asm] = self.assembly_powers[z,asm] * core_avg_linear_heatrate * mesh[z]

        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 not np.isclose(est_total_power, total_power, rtol=1e-3):

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