Commit f4c409cf authored by Phathanapirom U B's avatar Phathanapirom U B
Browse files

updates to nodal power calc for pow3d

parent 5e433f1a

SubKit/build/calc_nodal_powers.py

+24 −15
Original line number Diff line number Diff line
@@ -33,6 +33,8 @@ class Reactor: 
        self.__get_nodal_powers__()

        self.__perform_exit_checks__()



        self.state.create_dataset('pin_powers_nodal', data=self.nodal_powers)





    def __get_inserted_pin_powers__(self):

        """ Insert duplicated pin powers along node edges """
@@ -97,17 +99,22 @@ class Reactor: 
        Sum and group powers using nodal mask

        """



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

        self.nodal_powers = np.zeros( (self.core['core_map'].shape[0],self.core['core_map'].shape[1],4,self.n_z) )

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

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



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

        self.core_map = self.core['core_map'][()]

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

        asm_i,asm_j = 0,0



        for core_asm in range(0,np.max(self.core_map)): # all assembly numbers

            for asm in range(0,len(np.argwhere(self.asm_map==core_asm))): # number of assembly matching type core_asm

                asm_i,asm_j = np.argwhere(self.asm_map==core_asm)[asm,0],np.argwhere(self.asm_map==core_asm)[asm,1] # coordinate of assembly matching type

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

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

                    weighted_pin_powers = np.multiply(self.inserted_pin_powers[:,:,z,core_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))

                        self.nodal_powers[asm_i,asm_j,n,z] = np.sum( np.multiply(weighted_pin_powers,mask) )





    def __perform_exit_checks__(self):
@@ -119,22 +126,24 @@ class Reactor: 


        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)



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

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



        for core_asm in range(np.max(self.core_map)):

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

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

                    * core_avg_linear_heatrate * mesh[z]

                for asm in range(0,len(np.argwhere(self.asm_map==core_asm))):

                    asm_i,asm_j = np.argwhere(self.asm_map==core_asm)[asm,0],np.argwhere(self.asm_map==core_asm)[asm,1]

                    local_power[asm_i,asm_j,:,z] = self.nodal_powers[asm_i,asm_j,:,z] * core_avg_linear_heatrate * mesh[z]



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

        counter = 0

        for asm_i in range(0,self.core_map.shape[0]):

            for asm_j in range(0,self.core_map.shape[1]):

                est_total_power += np.sum( local_power[asm_i,asm_j,:,:] )

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



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