From 16c0243c480728f02cc91c64293e34a5ffb92b73 Mon Sep 17 00:00:00 2001 From: "Ryan M. Bergmann" Date: Fri, 18 Mar 2016 21:28:48 +0100 Subject: [PATCH] fixed some nu stuff... not many error messages now --- benchmarks.cpp | 2 +- pop_fission.cu | 2 +- unionize.py | 35 ++++++++++++++++++++++++++++------- 3 files changed, 30 insertions(+), 9 deletions(-) diff --git a/benchmarks.cpp b/benchmarks.cpp index 603cf20..1606a93 100644 --- a/benchmarks.cpp +++ b/benchmarks.cpp @@ -507,7 +507,7 @@ int main(int argc, char* argv[]){ bc = 1; runtype = "criticality"; if(card0name.compare(argv[2])==0){ - N = int(1e6); + N = int(1e5); dev = 0; } else if(card1name.compare(argv[2])==0){ diff --git a/pop_fission.cu b/pop_fission.cu index d2cac84..babdb90 100644 --- a/pop_fission.cu +++ b/pop_fission.cu @@ -133,7 +133,7 @@ __global__ void pop_fission_kernel(unsigned N, cross_section_data* d_xsdata, par // get interpolated beta value, beta = nu_d / nu_t beta = interpolate_linear_energy( this_E, e0, e1, nu_d0, nu_d1 ) / interpolate_linear_energy( this_E, e0, e1, nu_t0, nu_t1 ) ; - if(this_E > e1 | this_E < e0){printf("OUTSIDE bounds in pop_fission! this_E %6.4E e0 %6.4E e1 %6.4E\n",this_E,e0,e1);} + if(this_E > e1 | this_E < e0){printf("OUTSIDE bounds in pop_fission! this_E %6.4E e0 %6.4E e1 %6.4E col %u row %u\n",this_E,e0,e1,this_col,this_row);} // write new histories for this yield number for(unsigned k=0 ; k < this_yield ; k++ ){ diff --git a/unionize.py b/unionize.py index 3d71f9c..11e5a78 100644 --- a/unionize.py +++ b/unionize.py @@ -146,7 +146,15 @@ def _unionize(self): print " --------- unionizing grid --------- " for table in self.tables: + # main xs self.MT_E_grid=numpy.union1d(self.MT_E_grid,table.energy) + # nu if present + if hasattr(table,"nu_t_energy"): + self.MT_E_grid=numpy.union1d(self.MT_E_grid,table.nu_t_energy) + if hasattr(table,"nu_d_energy"): + self.MT_E_grid=numpy.union1d(self.MT_E_grid,table.nu_d_energy) + if hasattr(table,"nu_p_energy"): + self.MT_E_grid=numpy.union1d(self.MT_E_grid,table.nu_p_energy) # unionize the scattering energies in as well! if present of course for MT in table.reactions: rxn = table.reactions[MT] @@ -381,14 +389,18 @@ def _get_scatter_data(self,row,col): if nu_t_upper_index == None: nu_t_upper_index = len(table.nu_t_energy)-1 nu_t_lower_index = len(table.nu_t_energy)-1 + above_last_t = True else: nu_t_lower_index = nu_t_upper_index - 1 + above_last_t = False if nu_d_upper_index == None: nu_d_upper_index = len(table.nu_d_energy)-1 nu_d_lower_index = len(table.nu_d_energy)-1 + above_last_d = True else: nu_d_lower_index = nu_d_upper_index - 1 + above_last_d = False # make sure above threshold if nu_t_lower_index < 0: @@ -440,8 +452,8 @@ def _get_scatter_data(self,row,col): upper_nu_t_intt = table.nu_t_interp_INT[nu_p_upper_index] lower_pre_intt = table.nu_d_energy_dist[0].intt[0] upper_pre_intt = table.nu_d_energy_dist[0].intt[1] - lower_pre_law = table.nu_d_energy_dist[0].law - upper_pre_law = table.nu_d_energy_dist[0].law + lower_pre_law = table.nu_d_energy_dist[0].law + upper_pre_law = table.nu_d_energy_dist[0].law # set values in vars lower_law = -1 @@ -452,10 +464,18 @@ def _get_scatter_data(self,row,col): upper_erg = min(upper_e_t,upper_e_d) # take narrowest interval # evaluate nu on this interval - lower_nu_t = lower_nu_t_grid + (lower_erg - lower_e_t)/(upper_e_t - lower_e_t) * (upper_nu_t_grid - lower_nu_t_grid) - lower_nu_d = lower_nu_d_grid + (lower_erg - lower_e_d)/(upper_e_d - lower_e_d) * (upper_nu_d_grid - lower_nu_d_grid) - upper_nu_t = lower_nu_t_grid + (upper_erg - lower_e_t)/(upper_e_t - lower_e_t) * (upper_nu_t_grid - lower_nu_t_grid) - upper_nu_d = lower_nu_d_grid + (upper_erg - lower_e_d)/(upper_e_d - lower_e_d) * (upper_nu_d_grid - lower_nu_d_grid) + if above_last_t: + lower_nu_t = upper_nu_t_grid + upper_nu_t = upper_nu_t_grid + else: + lower_nu_t = lower_nu_t_grid + (lower_erg - lower_e_t)/(upper_e_t - lower_e_t) * (upper_nu_t_grid - lower_nu_t_grid) + upper_nu_t = lower_nu_t_grid + (upper_erg - lower_e_t)/(upper_e_t - lower_e_t) * (upper_nu_t_grid - lower_nu_t_grid) + if above_last_d: + lower_nu_d = upper_nu_d_grid + upper_nu_d = upper_nu_d_grid + else: + lower_nu_d = lower_nu_d_grid + (lower_erg - lower_e_d)/(upper_e_d - lower_e_d) * (upper_nu_d_grid - lower_nu_d_grid) + upper_nu_d = lower_nu_d_grid + (upper_erg - lower_e_d)/(upper_e_d - lower_e_d) * (upper_nu_d_grid - lower_nu_d_grid) lower_len = numpy.array([lower_nu_t,lower_nu_d]) upper_len = numpy.array([upper_nu_t,upper_nu_d]) @@ -484,7 +504,8 @@ def _get_scatter_data(self,row,col): #print lower_cdf[ pre_position], lower_cdf[ pre_position + lower_pdf[6]] ,lower_cdf[ pre_position+ lower_pdf[6]*2] # next index - if max(nu_t_upper_index,nu_d_upper_index) == max(len(table.nu_t_energy),len(table.nu_d_energy))-1 : # above last dist energy bin + if above_last_d and above_last_t: + #if max(nu_t_upper_index,nu_d_upper_index) == max(len(table.nu_t_energy),len(table.nu_d_energy))-1 : # above last dist energy bin next_dex = len(self.MT_E_grid) else: next_dex = next((i for i, x in enumerate(upper_erg <= self.MT_E_grid) if x), len(self.MT_E_grid))