diff --git a/pyne/ace.pyx b/pyne/ace.pyx index 3241a032a1..90ff6f5b48 100644 --- a/pyne/ace.pyx +++ b/pyne/ace.pyx @@ -25,6 +25,11 @@ from collections import OrderedDict cimport numpy as np import numpy as np +import math +from numpy.random import rand +from scipy.constants import physical_constants +# Complementary error function used in doppler broadening +from scipy.special import erfc from bisect import bisect_right from pyne cimport nucname @@ -37,6 +42,10 @@ cdef bint NP_LE_V15 = int(np.__version__.split('.')[1]) <= 5 and np.__version__. warn(__name__ + " is not yet QA compliant.", QAWarning) +# Constants +sqrt_pi_inv = 1. / math.sqrt(math.pi) +kb = physical_constants['Boltzmann constant in eV/K'][0] * 1e-6 + def ascii_to_binary(ascii_file, binary_file): """Convert an ACE file in ASCII format (type 1) to binary format (type 2). @@ -110,7 +119,113 @@ def ascii_to_binary(ascii_file, binary_file): # Close binary file binary.close() +def _interpolation_tab1(x_in, x, y, interp_NBT = None, interp_INT = None): + """ + INTERPOLATE_TAB1 interpolates a function between two points based on + particular interpolation scheme. The data needs to be organized as a ENDF + TAB1 type function containing the interpolation regions, break points, and + tabulated x's and y's. + """ + + cdef int i, interp + if x_in <= x[0]: + return y[0] + elif x_in >= x[-1]: + return y[-1] + else: + index = np.searchsorted(x, x_in) - 1 + + # Determine interpolation scheme + if interp_NBT is None: + #Linear-linear + interp = 2 + else: + if len(interp_NBT) == 1: + interp = int(interp_INT[0]) + elif len(interp_NBT) > 1: + for i in range(len(interp_NBT)): + if index < interp_NBT[i]: + interp = int(interp_INT[i]) + break + + # Handle special case of histogram interpolation + if interp == 1: + return y[index] + + # Determine bounding values + x0 = x[index] + x1 = x[index + 1] + y0 = y[index] + y1 = y[index + 1] + + # Determine interpolation factor and interpolated value + if interp == 2: # linear-linear + r = (x_in - x0) / (x1 - x0) + return y0 + r * (y1 - y0) + elif interp == 3: # linear _log + r = math.log(x_in / x0) / math.log(x1 / x0) + return y0 + r * (y1 - y0) + elif interp == 4: # log-linear + r = (x_in - x0) / (x1 - x0) + return y0 * math.exp(r * math.log(y1 / y0)) + elif interp == 5: # log-log + r = math.log(x_in / x0) / math.log(x1 / x0) + return y0 * math.exp(r * math.log(y1 / y0)) + +def _tabular_sample(inter_flag, out, pdf, cdf): + """ + Sample variable based on inter_flag, which is a common format + in secondary angle and energy sampling. + """ + + cdef int i + r = rand() + + # Sample from cdf + i = np.searchsorted(cdf, r) - 1 + + # Check to make sure j is <= len(cdf) - 2 + i = min(i, len(cdf) - 2) + + if inter_flag == 1: + # Histogram + if pdf[i] > 0.0: + ans = out[i] + (r - cdf[i]) / pdf[i] + else: + ans = out[i] + elif inter_flag == 2: + # Linear-linear + frac = (pdf[i + 1] - pdf[i]) / (out[i + 1] - out[i]) + if frac == 0.0: + ans = out[i] + (r - cdf[i]) / pdf[i] + else: + tmp = max(0, pdf[i] ** 2 + 2.0 * frac * (r - cdf[i])) ** 0.5 + ans = out[i] + (tmp - pdf[i]) / frac + return ans +def _find_index(value, array): + """ + Find bin of array containing value and calculate interpolation factor. + If the value is outside the range of array, choose the first or last bin. + """ + + if (value <= array[0]): + i = 0 + frac = 0.0 + elif (value >= array[-1]): + i = -2 + frac = 1.0 + else: + i = np.searchsorted(array, value) - 1 + frac = (value - array[i]) / (array[i + 1] - array[i]) + return i, frac + +def set_seed(seed): + """ + Set numpy random seed for deterministic testing. + """ + np.random.seed(seed) + class Library(object): """ A Library objects represents an ACE-formatted file which may contain @@ -517,10 +632,12 @@ class NeutronTable(AceTable): self.reactions[2] = elastic_scatter # Create all other reactions with MT values - mts = np.asarray(self.xss[self.jxs[3]:self.jxs[3] + n_reactions], dtype=int) - qvalues = np.asarray(self.xss[self.jxs[4]:self.jxs[4] + - n_reactions], dtype=float) - tys = np.asarray(self.xss[self.jxs[5]:self.jxs[5] + n_reactions], dtype=int) + mts = np.asarray(self.xss[self.jxs[3]:self.jxs[3] + n_reactions], + dtype=int) + qvalues = np.asarray(self.xss[self.jxs[4]:self.jxs[4] + n_reactions], + dtype=float) + tys = np.asarray(self.xss[self.jxs[5]:self.jxs[5] + n_reactions], + dtype=int) # Create all reactions other than elastic scatter reactions = [(mt, Reaction(mt, self)) for mt in mts] @@ -548,7 +665,8 @@ class NeutronTable(AceTable): self.jxs[7] + loc + 1 + n_energies] def _read_nu(self): - """Read the NU block -- this contains information on the prompt + """ + Read the NU block -- this contains information on the prompt and delayed neutron precursor yields, decay constants, etc """ cdef int ind, i, jxs2, KNU, LNU, NR, NE, NC @@ -567,19 +685,19 @@ class NeutronTable(AceTable): # Polynomial function form of nu if LNU == 1: self.nu_t_type = "polynomial" - NC = int(self.xss[KNU+1]) - coeffs = self.xss[KNU+2 : KNU+2+NC] + nc = int(self.xss[KNU+1]) + self.nu_t_coeffs = self.xss[KNU+2:KNU+2+nc] # Tabular data form of nu elif LNU == 2: self.nu_t_type = "tabular" NR = int(self.xss[KNU+1]) if NR > 0: - interp_NBT = self.xss[KNU+2 : KNU+2+NR ] - interp_INT = self.xss[KNU+2+NR : KNU+2+2*NR] + self.nu_t_interp_nbt = self.xss[KNU+2:KNU+2+NR] + self.nu_t_interp_int = self.xss[KNU+2+NR:KNU+2+2*NR] NE = int(self.xss[KNU+2+2*NR]) - self.nu_t_energy = self.xss[KNU+3+2*NR : KNU+3+2*NR+NE ] - self.nu_t_value = self.xss[KNU+3+2*NR+NE : KNU+3+2*NR+2*NE] + self.nu_t_energy = self.xss[KNU+3+2*NR:KNU+3+2*NR+NE] + self.nu_t_value = self.xss[KNU+3+2*NR+NE:KNU+3+2*NR+2*NE] # Both prompt nu and total nu elif self.xss[jxs2] < 0: KNU = jxs2 + 1 @@ -588,19 +706,19 @@ class NeutronTable(AceTable): # Polynomial function form of nu if LNU == 1: self.nu_p_type = "polynomial" - NC = int(self.xss[KNU+1]) - coeffs = self.xss[KNU+2 : KNU+2+NC] + nc = int(self.xss[KNU+1]) + self.nu_p_coeffs = self.xss[KNU+2:nc] # Tabular data form of nu elif LNU == 2: self.nu_p_type = "tabular" NR = int(self.xss[KNU+1]) if NR > 0: - interp_NBT = self.xss[KNU+2 : KNU+2+NR ] - interp_INT = self.xss[KNU+2+NR : KNU+2+2*NR] + self.nu_p_interp_nbt = self.xss[KNU+2:KNU+2+NR] + self.nu_p_interp_int = self.xss[KNU+2+NR:KNU+2+2*NR] NE = int(self.xss[KNU+2+2*NR]) - self.nu_p_energy = self.xss[KNU+3+2*NR : KNU+3+2*NR+NE ] - self.nu_p_value = self.xss[KNU+3+2*NR+NE : KNU+3+2*NR+2*NE] + self.nu_p_energy = self.xss[KNU+3+2*NR:KNU+3+2*NR+NE] + self.nu_p_value = self.xss[KNU+3+2*NR+NE:KNU+3+2*NR+2*NE] KNU = jxs2 + int(abs(self.xss[jxs2])) + 1 LNU = int(self.xss[KNU]) @@ -608,30 +726,30 @@ class NeutronTable(AceTable): # Polynomial function form of nu if LNU == 1: self.nu_t_type = "polynomial" - NC = int(self.xss[KNU+1]) - coeffs = self.xss[KNU+2 : KNU+2+NC] + nc = int(self.xss[KNU+1]) + self.nu_t_coeffs = self.xss[KNU+2:KNU+2+nc] # Tabular data form of nu elif LNU == 2: self.nu_t_type = "tabular" NR = int(self.xss[KNU+1]) if NR > 0: - interp_NBT = self.xss[KNU+2 : KNU+2+NR ] - interp_INT = self.xss[KNU+2+NR : KNU+2+2*NR] + self.nu_t_interp_nbt = self.xss[KNU+2:KNU+2+NR] + self.nu_t_interp_int = self.xss[KNU+2+NR:KNU+2+2*NR] NE = int(self.xss[KNU+2+2*NR]) - self.nu_t_energy = self.xss[KNU+3+2*NR : KNU+3+2*NR+NE ] - self.nu_t_value = self.xss[KNU+3+2*NR+NE : KNU+3+2*NR+2*NE] + self.nu_t_energy = self.xss[KNU+3+2*NR:KNU+3+2*NR+NE] + self.nu_t_value = self.xss[KNU+3+2*NR+NE:KNU+3+2*NR+2*NE] # Check for delayed nu data if self.jxs[24] > 0: KNU = self.jxs[24] NR = int(self.xss[KNU+1]) if NR > 0: - interp_NBT = self.xss[KNU+2 : KNU+2+NR ] - interp_INT = self.xss[KNU+2+NR : KNU+2+2*NR] + self.nu_d_interp_nbt = self.xss[KNU+2:KNU+2+NR] + self.nu_d_interp_int = self.xss[KNU+2+NR:KNU+2+2*NR] NE = int(self.xss[KNU+2+2*NR]) - self.nu_d_energy = self.xss[KNU+3+2*NR : KNU+3+2*NR+NE ] - self.nu_d_value = self.xss[KNU+3+2*NR+NE : KNU+3+2*NR+2*NE] + self.nu_d_energy = self.xss[KNU+3+2*NR:KNU+3+2*NR+NE] + self.nu_d_value = self.xss[KNU+3+2*NR+NE:KNU+3+2*NR+2*NE] # Delayed neutron precursor distribution self.nu_d_precursor_const = {} @@ -643,11 +761,11 @@ class NeutronTable(AceTable): self.nu_d_precursor_const[group] = self.xss[i] NR = int(self.xss[i+1]) if NR > 0: - interp_NBT = self.xss[i+2 : i+2+NR] - interp_INT = self.xss[i+2+NR : i+2+2*NR] + interp_NBT = self.xss[i+2:i+2+NR] + interp_INT = self.xss[i+2+NR:i+2+2*NR] NE = int(self.xss[i+2+2*NR]) - self.nu_d_precursor_energy[group] = self.xss[i+3+2*NR : i+3+2*NR+NE ] - self.nu_d_precursor_prob[group] = self.xss[i+3+2*NR+NE : i+3+2*NR+2*NE] + self.nu_d_precursor_energy[group] = self.xss[i+3+2*NR:i+3+2*NR+NE ] + self.nu_d_precursor_prob[group] = self.xss[i+3+2*NR+NE:i+3+2*NR+2*NE] i = i+3+2*NR+2*NE # Energy distribution for delayed fission neutrons @@ -659,12 +777,62 @@ class NeutronTable(AceTable): location_start, delayed_n=True) self.nu_d_energy_dist.append(energy_dist) - + def sample_nu(self, e): + """ + Sample nu based on incident neutron energy e, + return a tuple of length-3, format of [nu_total, nu_prompt, nu_delay]. + If no nu data is presented for the corresponding nu, return None instead + This implementation was inspired by the nu sampler in OpenMC v0.7.1. + """ + + cdef int i + # No nu data + if self.jxs[2] == 0: + return (None, None, None) + + # total nu always exists + if self.nu_t_type == "polynomial": + nu_t = 0.0 + for i in range(len(self.nu_t_coeffs)): + nu_t += self.nu_t_coeffs[i] * e ** i + elif self.nu_t_type == "tabular": + if hasattr(self, 'nu_t_interp_nbt'): + nu_t = _interpolation_tab1(e, self.nu_t_energy, self.nu_t_value, \ + self.nu_t_interp_nbt, self.nu_t_interp_int) + else: + nu_t = _interpolation_tab1(e, self.nu_t_energy, self.nu_t_value) + + # Prompt nu + if hasattr(self, 'nu_p_type'): + if self.nu_p_type == "polynomial": + nu_p = 0.0 + for i in range(len(self.nu_p_coeffs)): + nu_p += self.nu_p_coeffs[i] * e ** i + elif self.nu_p_type == 'tabular': + if hasattr(self, 'nu_p_interp_nbt'): + nu_p = _interpolation_tab1(e, self.nu_p_energy, self.nu_p_value, \ + self.nu_p_interp_nbt, self.nu_p_interp_int) + else: + nu_p = _interpolation_tab1(e, self.nu_p_energy, self.nu_p_value) + else: + nu_p = None + + # Delay nu + if hasattr(self, 'nu_d_energy'): + if hasattr(self, 'nu_d_interp_nbt'): + nu_d = _interpolation_tab1(e, self.nu_d_energy, self.nu_d_value, \ + self.nu_d_interp_nbt, self.nu_d_interp_int) + else: + nu_d = _interpolation_tab1(e, self.nu_d_energy, self.nu_d_value) + else: + nu_d = None + return (nu_t, nu_p, nu_d) + def _read_angular_distributions(self): """Find the angular distribution for each reaction MT """ cdef int ind, i, j, n_reactions, n_energies, n_bins - cdef dict ang_cos, ang_pdf, ang_cdf + cdef dict ang_cos, ang_pdf, ang_cdf, jj # Number of reactions with secondary neutrons (including elastic # scattering) @@ -673,18 +841,17 @@ class NeutronTable(AceTable): # Angular distribution for all reactions with secondary neutrons for i, reaction in enumerate(list(self.reactions.values())[:n_reactions]): loc = int(self.xss[self.jxs[8] + i]) - # Check if angular distribution data exist if loc == -1: - # Angular distribution data are specified through LAWi - # = 44 in the DLW block + # Angular distribution data are specified through ACE law, + # e.g., law 44 and law 61, in the DLW block continue elif loc == 0: - # No angular distribution data are given for this - # reaction, isotropic scattering is asssumed (in CM if - # TY < 0 and in LAB if TY > 0) + # No angular distribution data are given for this reaction, + # isotropic scattering is asssumed (in CM if TY < 0 and in LAB + # if TY > 0) + reaction.aflag = 'iso' continue - ind = self.jxs[9] + loc # Number of energies at which angular distributions are tabulated @@ -696,20 +863,21 @@ class NeutronTable(AceTable): # Read locations for angular distributions locations = np.asarray(self.xss[ind:ind + n_energies], dtype=int) + reaction.ang_locations = locations ind += n_energies ang_cos = {} ang_pdf = {} ang_cdf = {} + jj = {} for j, location in enumerate(locations): if location > 0: # Equiprobable 32 bin distribution - # print([reaction,'equiprobable']) ang_cos[j] = self.xss[ind:ind + 33] ind += 33 elif location < 0: # Tabular angular distribution - JJ = int(self.xss[ind]) + jj[j] = int(self.xss[ind]) n_bins = int(self.xss[ind + 1]) ind += 2 ang_dat = self.xss[ind:ind + 3*n_bins] @@ -725,7 +893,8 @@ class NeutronTable(AceTable): reaction.ang_cos = ang_cos reaction.ang_pdf = ang_pdf reaction.ang_cdf = ang_cdf - + reaction.jj = jj + def _read_energy_distributions(self): """Determine the energy distribution for secondary neutrons for each reaction MT @@ -749,7 +918,8 @@ class NeutronTable(AceTable): location_start. """ - cdef int ind, i, n_reactions, NE, n_regions, location_next_law, law, location_data, NPE, NPA + cdef int ind, i, n_reactions, NE, n_regions, location_next_law, law, \ + location_data, NPE, NPA # Create EnergyDistribution object edist = EnergyDistribution() @@ -795,6 +965,9 @@ class NeutronTable(AceTable): dat.shape = (2, n_regions) edist.NBT, edist.INT = dat ind += 2 * n_regions + raise NotImplementedError('Multiple interpolation regions ' + \ + 'not yet supported for tabular ' + \ + 'equiprobable energy distributions.') # Number of outgoing energies in each E_out table NE = int(self.xss[ind]) @@ -803,11 +976,13 @@ class NeutronTable(AceTable): # Read E_out tables NET = int(self.xss[ind]) - dat = self.xss[ind+1:ind+1+3*NET] - dat.shape = (3, NET) - self.e_dist_energy_out1, self.e_dist_energy_out2, \ - self.e_dist_energy_outNE = dat - ind += 1 + 3 * NET + # Each incident energy has a corresponding outgoing energy array, + # not just 3. The index in the instruction is from incoming energy + # bin 1, 2, ..., NE, not just 1, 2, NE + dat = self.xss[ind + 1 : ind + 1 + NE * NET] + dat.shape = (NE, NET) + edist.energy_out = dat + ind += 1 + NE * NET elif law == 2: # Discrete photon energy self.e_dist_LP = int(self.xss[ind]) @@ -824,7 +999,7 @@ class NeutronTable(AceTable): if n_regions > 0: dat = np.asarray(self.xss[ind:ind+2*n_regions], dtype=int) dat.shape = (2, n_regions) - edist.NBT, edist.INT = dat + edist.nbt, edist.int = dat ind += 2 * n_regions # Number of outgoing energies in each E_out table @@ -835,6 +1010,7 @@ class NeutronTable(AceTable): nps = [] edist.intt = [] # Interpolation scheme (1=hist, 2=lin-lin) + edist.nd = [] # Number of discrete lines edist.energy_out = [] # Outgoing E grid for each incoming E edist.pdf = [] # Probability dist for " " " edist.cdf = [] # Cumulative dist for " " " @@ -847,9 +1023,8 @@ class NeutronTable(AceTable): INTT = INTTp ND = 0 edist.intt.append(INTT) - #if ND > 0: - # print [reaction, ND, INTT] - + edist.nd.append(ND) + NP = int(self.xss[ind+1]) nps.append(NP) dat = self.xss[ind+2:ind+2+3*NP] @@ -1014,7 +1189,7 @@ class NeutronTable(AceTable): if n_regions > 0: dat = np.asarray(self.xss[ind:ind+2*n_regions], dtype=int) dat.shape = (2, n_regions) - edist.NBT, edist.INT = dat + edist.nbt, edist.int = dat ind += 2 * n_regions # Number of outgoing energies in each E_out table @@ -1025,6 +1200,7 @@ class NeutronTable(AceTable): nps = [] edist.intt = [] # Interpolation scheme (1=hist, 2=lin-lin) + edist.nd = [] # Number of discrete lines edist.energy_out = [] # Outgoing E grid for each incoming E edist.pdf = [] # Probability dist for " " " edist.cdf = [] # Cumulative dist for " " " @@ -1032,12 +1208,11 @@ class NeutronTable(AceTable): edist.ang = [] # Angular distribution slope for " " " for i in range(NE): INTTp = int(self.xss[ind]) - if INTTp > 10: - INTT = INTTp % 10 - ND = (INTTp - INTT)/10 - else: - INTT = INTTp - edist.intt.append(INTT) + # No matter INTTp > 10 or not, we can get the data + # At sample stage, use nd to determine if there is + # discrete lines error + edist.intt.append(INTTp % 10) + edist.nd.append(INTTp // 10) NP = int(self.xss[ind+1]) nps.append(NP) @@ -1050,7 +1225,7 @@ class NeutronTable(AceTable): edist.cdf.append(dat[2]) edist.frac.append(dat[3]) edist.ang.append(dat[4]) - ind += 5*NP + ind += 5 * NP # convert to arrays if possible edist.intt = np.array(edist.intt) @@ -1067,20 +1242,22 @@ class NeutronTable(AceTable): if n_regions > 0: dat = np.asarray(self.xss[ind:ind+2*n_regions], dtype=int) dat.shape = (2, n_regions) - edist.NBT, edist.INT = dat + edist.nbt, edist.int = dat ind += 2 * n_regions - # Number of outgoing energies in each E_out table + # Number of incoming energies NE = int(self.xss[ind]) edist.energy_in = self.xss[ind+1:ind+1+NE] L = np.asarray(self.xss[ind+1+NE:ind+1+2*NE], dtype=int) ind += 1 + 2*NE - + npes = [] edist.intt = [] # Interpolation scheme (1=hist, 2=lin-lin) + edist.nd = [] # Number of discrete lines edist.energy_out = [] # Outgoing E grid for each incoming E edist.pdf = [] # Probability dist for " " " edist.cdf = [] # Cumulative dist for " " " + edist.lc = [] # Use lc to determine ang. data exist or not npas = [] edist.a_dist_intt = [] @@ -1089,22 +1266,21 @@ class NeutronTable(AceTable): edist.a_dist_cdf = [] for i in range(NE): INTTp = int(self.xss[ind]) - if INTTp > 10: - INTT = INTTp % 10 - ND = (INTTp - INTT)/10 - else: - INTT = INTTp - edist.intt.append(INTT) - + # At sample stage, use nd to determine if there is + # discrete lines present + edist.intt.append(INTTp % 10) + edist.nd.append(INTTp // 10) + # Secondary energy distribution NPE = int(self.xss[ind+1]) - npes.append(NPE) + npes.append(NPE) dat = self.xss[ind+2:ind+2+4*NPE] dat.shape = (4, NPE) edist.energy_out.append(dat[0]) edist.pdf.append(dat[1]) edist.cdf.append(dat[2]) - LC = np.asarray(dat[3], dtype=int) + # For case of lc[e_out] = 0 + edist.lc.append(np.asarray(dat[3], dtype=int)) ind += 2 + 4*NPE # Secondary angular distribution @@ -1113,15 +1289,23 @@ class NeutronTable(AceTable): edist.a_dist_pdf.append([]) edist.a_dist_cdf.append([]) for j in range(NPE): - edist.a_dist_intt[-1].append(int(self.xss[ind])) - NPA = int(self.xss[ind+1]) - npas.append(NPA) - dat = self.xss[ind+2:ind+2+3*NPA] - dat.shape = (3, NPA) - edist.a_dist_mu_out[-1].append(dat[0]) - edist.a_dist_pdf[-1].append(dat[1]) - edist.a_dist_cdf[-1].append(dat[2]) - ind += 2 + 3*NPA + if edist.lc[i][j] > 0: + edist.a_dist_intt[-1].append(int(self.xss[ind])) + NPA = int(self.xss[ind+1]) + npas.append(NPA) + dat = self.xss[ind+2:ind+2+3*NPA] + dat.shape = (3, NPA) + edist.a_dist_mu_out[-1].append(dat[0]) + edist.a_dist_pdf[-1].append(dat[1]) + edist.a_dist_cdf[-1].append(dat[2]) + ind += 2 + 3*NPA + elif edist.lc[i][j] == 0: + # Insert None for correct access + edist.a_dist_intt[-1].append(None) + edist.a_dist_mu_out[-1].append(None) + edist.a.dist_pdf[-1].append(None) + edist.a_dist_cdf[-1].append(None) + # convert to arrays if possible edist.intt = np.array(edist.intt) @@ -1137,6 +1321,7 @@ class NeutronTable(AceTable): edist.a_dist_mu_out = np.array(edist.a_dist_mu_out) edist.a_dist_pdf = np.array(edist.a_dist_pdf) edist.a_dist_cdf = np.array(edist.a_dist_cdf) + elif law == 66: # N-body phase space distribution (ENDF File 6 Law 6) edist.nbodies = int(self.xss[ind]) @@ -1163,7 +1348,6 @@ class NeutronTable(AceTable): return edist - def _read_gpd(self): """Read total photon production cross section. """ @@ -1556,9 +1740,121 @@ class Reaction(object): self.IE = 0 # Energy grid index self.sigma = [] # Cross section values - def broaden(self, T_high): - pass - + def broaden(self, t_high): + """Doppler broaden cross section using SIGMA1 or piecewise-linear exact + integration method (see "Exact Doppler Broadening of Tabulated Cross + Sections," Nucl. Sci. Eng. 60, 199-229 (1976) && "Comparison of + algorithms for Doppler broadening pointwise tabulated cross sections", + Annals of Nuclear Energy 75 (2015) 358-364). The code implementation is + inspired by the doppler module in OpenMC v0.7.1. + + Parameter + ---------- + t_high : float + Broadened-to temperature in K + + Result + ---------- + sigmaNew : np array + Broadened cross sections + """ + cdef int i, k + cdef double y, y_sq, y_inv, y_inv_sq, a, b, sigma, ak, bk, slope + # Broadened cross sections + sigmaNew = np.zeros(len(self.sigma)) + # temperature difference in K, note the temp in neutron table is + # in MeV, using boltzmann constant to convert to K + t = t_high - self.table.temp / kb + alpha = self.table.awr / (kb * t) + xs = self.sigma + energy = self.table.energy[self.IE : self.IE + len(xs)] + x = np.sqrt(alpha * energy) + fa = np.zeros(5) + fb = np.zeros(5) + + for i in range(len(xs)): + sigma = 0.0 + y = x[i] + y_sq = y * y + y_inv = 1. / y + y_inv_sq = y_inv / y + + # Evaluate sigma*(y, T) from x[k] - y = 0 to -4 + k = i + a = 0.0 + self._calculate_f(fa, a) + while (a >= -4.0 and k > 0): + # Move to next point + fb[:] = fa + k -= 1 + a = x[k] - y + self._calculate_f(fa, a) + h = fa - fb + # Calculate a[k], b[k], and slope terms + ak = y_inv_sq * h[2] + 2 * y_inv * h[1] + h[0] + bk = y_inv_sq * h[4] + 4 * y_inv * h[3] + 6 * h[2] + \ + 4 * y * h[1] + y_sq * h[0] + slope = (xs[k + 1] - xs[k]) / (x[k + 1] ** 2 - x[k] ** 2) + # Add contribution to broadened cross section + sigma += ak * (xs[k] - slope * x[k] ** 2) + slope * bk + # Extend cross section to 0 assuming 1/v shape + if (k == 0 and a >= -4.0): + fb[:] = fa + a = -y + self._calculate_f(fa, a) + h = fa - fb + sigma += xs[k] * x[k] * (y_inv_sq * h[1] + y_inv * h[0]) + # Evaluate sigma*(y, T) from x[k] - y = 0 to 4 + k = i + b = 0.0 + self._calculate_f(fb, b) + while (b <= 4.0 and k < len(xs) - 1): + # Move to next point + fa[:] = fb + k += 1 + b = x[k] - y + # Calculate f and h functions + self._calculate_f(fb, b) + h = fa - fb + ak = y_inv_sq * h[2] + 2 * y_inv * h[1] + h[0] + bk = y_inv_sq * h[4] + 4 * y_inv * h[3] + 6 * h[2] + \ + 4 * y * h[1] + y_sq * h[0] + slope = (xs[k] - xs[k - 1]) / (x[k] ** 2 - x[k - 1] ** 2) + sigma += ak * (xs[k] - slope * x[k] ** 2) + slope * bk + # Extend cross section to infinity assuming constant shape + if (k == len(xs) - 1 and b <= 4.0): + a = x[k] - y + self._calculate_f(fa, a) + sigma += xs[k] * (y_inv_sq * fa[2] + 2 * y_inv * fa[1] + fa[0]) + + # Evaluate second term from x[k] + y = 0 to 4 + if (y <= 4.0): + # swip signs on y + y = -y + y_inv = -y_inv + k = 0 + # Calculate a and b based on 0 and x[0] + a = -y + b = x[k] - y + self._calculate_f(fa, a) + self._calculate_f(fb, b) + h = fa - fb + sigma = sigma - xs[k] * x[k] * (y_inv_sq * h[1] + y_inv * h[0]) + while (b <= 4.0): + fa[:] = fb + k += 1 + b = x[k] - y + self._calculate_f(fb, b) + h = fa - fb + ak = y_inv_sq * h[2] + 2 * y_inv * h[1] + h[0] + bk = y_inv_sq * h[4] + 4 * y_inv * h[3] + 6 * h[2] + \ + 4 * y * h[1] + y_sq * h[0] + slope = (xs[k] - xs[k - 1]) / (x[k] ** 2 - x[k - 1] ** 2) + sigma = sigma - ak * (xs[k] - slope * x[k] ** 2) - \ + slope * bk + sigmaNew[i] = sigma + return sigmaNew + def threshold(self): """threshold() @@ -1573,8 +1869,457 @@ class Reaction(object): else: rep = "".format(self.MT) return rep - - + + def sample(self, e): + """Sample outgoing angle and energy based on incident neutron energy. + This implementation (together with the subfunctions used) were inspired + by the second-neutron-distribution sampler in OpenMC v0.7.1. + + Parameters + ---------- + e : float + incident neutron energy + + Results + ---------- + (mu, e_out) : tuple + secondary neutron's angle, mu, and energy, e_out + """ + if self.MT == 2: + # Isotropic scattering + mu = self._sample_mu(e) + return (mu, e) + + edist = self.energy_dist + if edist.law == 1: + # Tabular equiprobable energy sample + i, frac = _find_index(e, edist.energy_in) + net = len(edist.energy_out[0]) + # Sample outgoing energy bin + k = int(net * rand()) + # Determine min and max of outgoing energy + ei1 = edist.energy_out[i, 0] + eik = edist.energy_out[i, -1] + ei11 = edist.energy_out[i + 1, 0] + ei1k = edist.energy_out[i + 1, -1] + e1 = ei1 + frac * (ei11 - ei1) + ek = eik + frac * (ei1k - eik) + # Select incident energy bin + if (rand() < frac): + l = i + 1 + else: + l = i + # Determine elk and elk1 + elk = edist.energy_out[l, k] + elk1 = edist.energy_out[l, k + 1] + # Compute unbounded outgoing energy + e_out = elk + rand() * (elk1 - elk) + # Interpolate between incident energy bins i and i + 1 + if l == i: + e_out = e1 + (e_out - ei1) * (ek - e1) / (eik - ei1) + else: + e_out = e1 + (e_out - ei11) * (ek - e1) / (ei1k - ei11) + mu = self._sample_mu(e) + return (mu, e_out) + + elif edist.law == 3: + # Inelastic level scattering + # Note the sampled outgoing energy is in CM system + e_out = edist.data[1] * (e - edist.data[0]) + mu = self._sample_mu(e) + return(mu, e_out) + + elif edist.law == 4: + # Continuous Tabular Distribution + if hasattr(edist, 'int'): + if len(edist.int) > 1: + raise NotImplementedError('Multiple interpolation regions not yet supported' + ' for continuous tabular energy distributions.') + histogram_interp = (edist.int[0] == 1) + else: + histogram_interp = False + + # Find energy bin and calculate interpolation factor -- if the energy is + # outside the range of the tabulated energies, choose the first or last bins + i, f = _find_index(e, edist.energy_in) + + # Sample between the ith and (i+1)th bin + if (histogram_interp): + if e >= edist.energy_in[-1]: + l = i + 1 + else: + l = i + else: + if (f > rand()): + l = i + 1 + else: + l = i + + # check for discrete lines present + if edist.nd[l] != 0: + raise NotImplementedError('Discrete lines in continuous tabular ' + 'distribution not yet supported.') + + # Interpolation for energy E1 and EK + E_i_1 = edist.energy_out[i][0] + E_i_K = edist.energy_out[i][-1] + + E_i1_1 = edist.energy_out[i+1][0] + E_i1_K = edist.energy_out[i+1][-1] + + E_1 = E_i_1 + f*(E_i1_1 - E_i_1) + E_K = E_i_K + f*(E_i1_K - E_i_K) + + E_out = _tabular_sample(edist.intt[l], edist.energy_out[l], \ + edist.pdf[l], edist.cdf[l]) + + # Interpolate between incident energy bins i and i + 1 + if not histogram_interp: + if (l == i): + E_out = E_1 + (E_out - E_i_1)*(E_K - E_1)/(E_i_K - E_i_1) + else: + E_out = E_1 + (E_out - E_i1_1)*(E_K - E_1)/(E_i1_K - E_i1_1) + + # Sample mu + mu = self._sample_mu(e) + return (mu, E_out) + + elif edist.law == 7: + # Simple Maxwell fission spectrum + if hasattr(edist, 'NBT'): + t = _interpolation_tab1(e, edist.energy_in, edist.T, + edist.NBT, edist.INT) + else: + t = _interpolation_tab1(e, edist.energy_in, edist.T) + e_out = self._maxwell_spectrum(t) + while e_out > e - edist.U: + e_out = self._maxwell_spectrum(t) + mu = self._sample_mu(e) + return (mu, e_out) + + elif edist.law == 9: + # Evaporation spectrum + if hasattr(edist, 'NBT'): + t = _interpolation_tab1(e, edist.energy_in, edist.T, + edist.NBT, edist.INT) + else: + t = _interpolation_tab1(e, edist.energy_in, edist.T) + w = (e - edist.U) / t + g = 1 - math.exp(-w) + e_out = -math.log((1 - g * rand()) * (1 - g * rand())) + while e_out > w: + e_out = -math.log((1 - g * rand()) * (1 - g * rand())) + e_out *= t + mu = self._sample_mu(e) + return (mu, e_out) + + elif edist.law == 11: + # Energy dependent Watt spectrum + if hasattr(edist, 'NBTa'): + a = _interpolation_tab1(e, edist.energya_in, edist.a, + edist.NBTa, edist.INTa) + else: + a = _interpolation_tab1(e, edist.energya_in, edist.a) + if hasattr(edist, 'NBTb'): + b = _interpolation_tab1(e, edist.energyb_in, edist.b, + edist.NBTb, edist.INTb) + else: + b = _interpolation_tab1(e, edist.energyb_in, edist.b) + e_out = self._watt_spectrum(a, b) + while e_out > e - edist.U: + e_out = self._watt_spectrum(a, b) + mu = self._sample_mu(e) + return (mu, e_out) + + elif edist.law == 66: + # N-body phase space distribution + e_max = (edist.massratio - 1) / edist.massratio * \ + (self.table.awr * e / (self.table.awr + 1) + self.Q) + x = self._maxwell_spectrum(1.) + if edist.nbodies == 3: + y = self._maxwell_spectrum(1.) + elif edist.nbodies == 4: + r1 = rand() + r2 = rand() + r3 = rand() + y = -math.log(r1 * r2 * r3) + elif edist.nbodies == 5: + r1 = rand() + r2 = rand() + r3 = rand() + r4 = rand() + r5 = rand() + r6 = rand() + y = -math.log(r1 * r2 * r3 * r4) - math.log(r5) * \ + math.cos(.5 * math.pi * r6) ** 2 + v = x / (x + y) + e_out = e_max * v + mu = self._sample_mu(e) + return (mu, e_out) + + elif edist.law == 44: + # Kalbach-87 Formalism + return self._sample_law44(e) + + elif edist.law == 61: + # Correlated energy and angle sample + return self._sample_law61(e) + + def _sample_mu(self, e): + """Sample the uncorrected outgoing angle + Parameters + ---------- + e: incident neutron energy + """ + + if hasattr(self, 'aflag'): + mu = 2.0 * rand() - 1.0 + return mu + else: + # Compute index and interpolation frac + i, r = _find_index(e, self.ang_energy_in) + + if (r > rand()): + i = i + 1 + + loc = self.ang_locations[i] + if loc == 0: + # Isotropic + return 2.0 * rand() - 1.0 + elif loc > 0: + # 32 equip bin + r1 = rand() + ii = 1 + int(32 * r1) + mui = self.ang_cos[ii - 1] + mui1 = self.ang_cos[ii] + mu = mui + (32 * r1 - ii) * (mui1 - mui) + + # Make sure mu is in range [-1,1] + if (abs(mu) > 1): + mu = np.sign(mu) + return mu + else: + # tabular + cos = self.ang_cos[i] + pdf = self.ang_pdf[i] + cdf = self.ang_cdf[i] + jj = self.jj[i] + mu = _tabular_sample(jj, cos, pdf, cdf) + + # Make sure mu is in range [-1,1] + if (abs(mu) > 1): + mu = np.sign(mu) + return mu + + def _sample_law44(self, e): + # Kalbach-87 Formalism (ENDF File 6 Law 1, LANG=2) + edist = self.energy_dist + # Interpolation scheme + if hasattr(edist, 'nbt'): + raise NotImplementedError('Multiple interpolation regions not yet supported ' + 'for Kalbach-Mann energy distributions') + + # Find energy bin and calculate interpolation factor -- if the energy is + # outside the range of the tabulated energies, choose the first or last bins + i, f = _find_index(e, edist.energy_in) + + # Sample between the ith and (i+1)th bin + if (f > rand()): + l = i + 1 + else: + l = i + + if edist.nd[l] != 0: + raise NotImplementedError('Discrete lines in Kalbach-Mann ' + 'distribution not yet supported.') + + # Interpolation for energy E1 and EK + E_i_1 = edist.energy_out[i][0] + E_i_K = edist.energy_out[i][-1] + + E_i1_1 = edist.energy_out[i+1][0] + E_i1_K = edist.energy_out[i+1][-1] + + E_1 = E_i_1 + f*(E_i1_1 - E_i_1) + E_K = E_i_K + f*(E_i1_K - E_i_K) + + # determine outgoing energy bin + r1 = rand() + cdf = edist.cdf[l] + k = np.searchsorted(cdf, r1) - 1 + c_k = cdf[k] + + # Check to make sure k <= len(cdf) - 2 + k = min(k, len(cdf) - 2) + + E_l_k = edist.energy_out[l][k] + p_l_k = edist.pdf[l][k] + if (edist.intt[l] == 1): + # Histogram interpolation + if (p_l_k > 0.0): + E_out = E_l_k + (r1 - c_k) / p_l_k + else: + E_out = E_l_k + + km_r = edist.frac[l][k] + km_a = edist.ang[l][k] + + elif (edist.intt[l] == 2): + # Linear-linear interpolation + E_l_k1 = edist.energy_out[l][k+1] + p_l_k1 = edist.pdf[l][k+1] + + frac = (p_l_k1 - p_l_k) / (E_l_k1 - E_l_k) + if (frac == 0.0): + E_out = E_l_k + (r1 - c_k) / p_l_k + else: + E_out = E_l_k + (max(0.0, p_l_k ** 2 + \ + 2.0 * frac * (r1-c_k)) ** 0.5 - p_l_k) / frac + + # Determine Kalbach-Mann parameters + km_r = edist.frac[l][k] + (E_out - E_l_k) / (E_l_k1 - E_l_k) * \ + (edist.frac[l][k+1] - edist.frac[l][k]) + km_a = edist.ang[l][k] + (E_out - E_l_k) / (E_l_k1 - E_l_k) * \ + (edist.ang[l][k+1] - edist.ang[l][k]) + + # Now interpolate between incident energy bins i and i + 1 + if (l == i): + E_out = E_1 + (E_out - E_i_1)*(E_K - E_1)/(E_i_K - E_i_1) + else: + E_out = E_1 + (E_out - E_i1_1)*(E_K - E_1)/(E_i1_K - E_i1_1) + + # Sampled correlated angle from Kalbach-Mann parameters + if (rand() > km_r): + t = (2.0 * rand() - 1.0) * math.sinh(km_a) + mu = math.log(t + (t * t + 1.0)**0.5) / km_a + else: + r1 = rand() + mu = math.log(r1 * math.exp(km_a) + (1.0 - r1) * math.exp(-km_a))/km_a + return (mu, E_out) + + def _sample_law61(self, e): + # Like 44, but tabular distribution instead of Kalbach-87 + + edist = self.energy_dist + + # Interpolation scheme + if hasattr(edist, 'nbt'): + raise NotImplementedError('Multiple interpolation regions not yet supported' + ' for correlated angle-energy distributions.') + + # find energy bin and calculate interpolation factor -- if the energy is + # outside the range of the tabulated energies, choose the first or last bins + i, r = _find_index(e, edist.energy_in) + + # Sample between the ith and (i+1)th bin + if (r > rand()): + l = i+1 + else: + l = i + + # check for discrete lines present + if edist.nd[l] != 0: + raise NotImplementedError('Discrete lines in correlated angle-energy' + ' distribution not yet supported.') + + # interpolation for energy E1 and EK + E_i_1 = edist.energy_out[i][0] + E_i_K = edist.energy_out[i][-1] + + E_i1_1 = edist.energy_out[i+1][0] + E_i1_K = edist.energy_out[i+1][-1] + + E_1 = E_i_1 + r*(E_i1_1 - E_i_1) + E_K = E_i_K + r*(E_i1_K - E_i_K) + + # Determine outgoing energy bin + r1 = rand() + cdf = edist.cdf[l] + k = np.searchsorted(cdf, r1) - 1 + c_k = cdf[k] + + # Make sure k <= len(cdf) - 2 + k = min(k, len(cdf) - 2) + + E_l_k = edist.energy_out[l][k] + p_l_k = edist.pdf[l][k] + if (edist.intt[l][k] == 1): + # Histogram interpolation + if(p_l_k > 0.0): + E_out = E_l_k + (r1 - c_k) / p_l_k + else: + E_out = E_l_k + elif(edist.intt[l][k] == 2): + # Linear-Linear interpolation + E_l_k1 = edist.energy_out[l][k+1] + p_l_k1 = edist.pdf[l][k+1] + + frac = (p_l_k1 - p_l_k)/(E_l_k1 - E_l_k) + if (frac == 0.0): + E_out = E_l_k + (r1 - c_k)/p_l_k + else: + E_out = E_l_k + (max(0.0, p_l_k**2 + 2.0*frac*(r1-c_k))**0.5 -\ + p_l_k)/frac + + # Now interpolate between incident energy bins i and i + 1 + if (l == i): + E_out = E_1 + (E_out - E_i_1)*(E_K - E_1)/(E_i_K - E_i_1) + else: + E_out = E_1 + (E_out - E_i1_1)*(E_K - E_1)/(E_i1_K - E_i1_1) + + # Find correlated angular distribution for closest outgoing energy bin + if(r1 - c_k < cdf[k+1] - r1): + if edist.lc[l][k] == 0: + mu = 2.0 * rand() - 1 + else: + jj = edist.a_dist_intt[l][k] + cos = edist.a_dist_mu_out[l][k] + pdf = edist.a.dist_pdf[l][k] + cdf = edist.a.dist.cdf[l][k] + mu = _tabular_sample(jj, cos, pdf, cdf) + + # Make sure mu is in range [-1,1] + if (abs(mu) > 1): + mu = np.sign(mu) + return (mu, E_out) + else: + if edist.lc[l][k+1] == 0: + mu = 2.0 * rand() - 1 + else: + jj = edist.a_dist_intt[l][k+1] + cos = edist.a_dist_mu_out[l][k+1] + pdf = edist.a.dist_pdf[l][k+1] + cdf = edist.a.dist.cdf[l][k+1] + mu = _tabular_sample(jj, cos, pdf, cdf) + + # Make sure mu is in range [-1,1] + if (abs(mu) > 1): + mu = np.sign(mu) + return (mu, E_out) + + def _maxwell_spectrum(self, t): + r1 = rand() + r2 = rand() + r3 = rand() + c = math.cos(0.5 * math.pi * r3) + e_out = -t * (math.log(r1) + math.log(r2) * c * c) + return e_out + + def _watt_spectrum(self, a, b): + w = self._maxwell_spectrum(a) + e_out = w + .25 * a ** 2 * b + (2 * rand() - 1) * \ + math.sqrt(a ** 2 * b * w) + return e_out + + def _calculate_f(self, f, a): + """Working hourse function used in Doppler broadening + """ + f[0] = .5 * erfc(a) + f[1] = .5 * sqrt_pi_inv * math.exp(-a * a) + f[2] = .5 * f[0] + a * f[1] + f[3] = f[1] * (1 + a * a) + f[4] = .75 * f[0] + f[1] * a * (1.5 + a * a) + class DosimetryTable(AceTable): def __init__(self, name, awr, temp): diff --git a/tests/test_ace.py b/tests/test_ace.py index fc9bd1abb6..52260c14d9 100644 --- a/tests/test_ace.py +++ b/tests/test_ace.py @@ -24,6 +24,10 @@ def setup(): with open('C012-n-2p0.ace','w') as f: f.writelines(lines) + + if not os.path.isfile('U235-n.ace'): + urllib.urlretrieve('ftp://ftp.nrg.eu/pub/www/talys/tendl2013/neutron_file/U/235/lib/endf/U235-n.ace', + 'U235-n.ace') def test_convert_c12(): @@ -94,7 +98,22 @@ def test_read_c12_binary(): assert_equal(table.reactions[2].sigma[0], 78.04874) assert_equal(table.reactions[2].sigma[-1], 1.00772) + +def test_sample_functions(): + u235 = pyne.ace.Library('U235-n.ace') + u235.read() + table = u235.tables['92235.00c'] + fission = table.reactions[18] + + pyne.ace.set_seed(1) + assert_almost_equal(table.sample_nu(1)[0], 2.525704) + (a, e) = fission.sample(1) + assert_almost_equal(a, -0.9997712503653102) + assert_almost_equal(e, 2.6309118671014353) def teardown(): if os.path.exists('C12-binary.ace'): os.remove('C12-binary.ace') + + if os.path.exists('U235-n.ace'): + os.remove('U235-n.ace')