diff --git a/openmc/data/multipole.py b/openmc/data/multipole.py index 230e944e4d..e40b0d0d0f 100644 --- a/openmc/data/multipole.py +++ b/openmc/data/multipole.py @@ -5,8 +5,12 @@ import warnings import os import h5py +import pickle import numpy as np from scipy.signal import find_peaks +import matplotlib +matplotlib.use("agg") +import matplotlib.pyplot as plt import openmc.checkvalue as cv from ..exceptions import DataError @@ -38,6 +42,7 @@ TEMPERATURE_LIMIT = 3000 # Logging control DETAILED_LOGGING = 2 + def _faddeeva(z): r"""Evaluate the complex Faddeeva function. @@ -141,16 +146,18 @@ def _broaden_wmp_polynomials(E, dopp, n): return factors + def _vectfit_xs(energy, ce_xs, mts, rtol=1e-3, atol=1e-5, orders=None, - n_vf_iter=30, log=False, path_out=None, **kwargs): - r"""Convert point-wise cross section to multipole data via Vector Fitting. + n_vf_iter=30, log=False, path_out=None): + """Convert point-wise cross section to multipole data via vector fitting. Parameters ---------- energy : np.ndarray Energy array ce_xs : np.ndarray - Point-wise cross sections to be fitted + Point-wise cross sections to be fitted, with shape (number of reactions, + number of energy points) mts : Iterable of int Reaction list rtol : float, optional @@ -162,15 +169,14 @@ def _vectfit_xs(energy, ce_xs, mts, rtol=1e-3, atol=1e-5, orders=None, n_vf_iter : int, optional Number of maximum VF iterations log : bool or int, optional - Whether to print running logs + Whether to print running logs (use int for verbosity control) path_out : str, optional - Path to save the figures - **kwargs - Additional keyword arguments + Path to save the figures to show discrepancies between the original and + fitted cross sections for different reactions Returns ------- - Tuple + tuple (poles, residues) """ @@ -184,10 +190,10 @@ def _vectfit_xs(energy, ce_xs, mts, rtol=1e-3, atol=1e-5, orders=None, raise ValueError('Inconsistent cross section data.') # construct test data: interpolate xs with finer grids - N_FINER = 10 - ne_test = (ne-1)*N_FINER + 1 + n_finer = 10 + ne_test = (ne - 1)*n_finer + 1 test_energy = np.interp(np.arange(ne_test), - np.arange(ne_test, step=N_FINER), energy) + np.arange(ne_test, step=n_finer), energy) test_energy[[0, -1]] = energy[[0, -1]] # avoid numerical issue test_xs_ref = np.zeros((nmt, ne_test)) for i in range(nmt): @@ -197,32 +203,35 @@ def _vectfit_xs(energy, ce_xs, mts, rtol=1e-3, atol=1e-5, orders=None, print("Energy: {:.3e} to {:.3e} eV ({} points)".format( energy[0], energy[-1], ne)) - # inputs - f = ce_xs * energy # sigma*E - s = np.sqrt(energy) # sqrt(E) + # transform xs (sigma) and energy (E) to f (sigma*E) and s (sqrt(E)) to be + # compatible with the multipole representation + f = ce_xs * energy + s = np.sqrt(energy) test_s = np.sqrt(test_energy) + + # inverse weighting is used for minimizing the relative deviation instead of + # absolute deviation in vector fitting weight = 1.0/f - # very small cross sections can lead to huge weights, which will harm the - # fitting accuracy - MIN_CROSS_SECTION = 1e-7 + # avoid too large weights which will harm the fitting accuracy + min_cross_section = 1e-7 for i in range(nmt): - if np.all(ce_xs[i]<=MIN_CROSS_SECTION): + if np.all(ce_xs[i] <= min_cross_section): weight[i] = 1.0 - elif np.any(ce_xs[i]<=MIN_CROSS_SECTION): - weight[i, ce_xs[i]<=MIN_CROSS_SECTION] = \ - max(weight[i, ce_xs[i]>MIN_CROSS_SECTION]) + elif np.any(ce_xs[i] <= min_cross_section): + weight[i, ce_xs[i] <= min_cross_section] = \ + max(weight[i, ce_xs[i] > min_cross_section]) # detect peaks (resonances) and determine VF order search range - peaks, _ = find_peaks(ce_xs[0]+ce_xs[1]) + peaks, _ = find_peaks(ce_xs[0] + ce_xs[1]) n_peaks = peaks.size if orders is not None: # make sure orders are even integers - orders = list(set([int(i/2)*2 for i in orders if i>=2])) + orders = list(set([int(i/2)*2 for i in orders if i >= 2])) else: lowest_order = max(2, 2*n_peaks) highest_order = max(200, 4*n_peaks) - orders = list(range(lowest_order, highest_order+1, 2)) + orders = list(range(lowest_order, highest_order + 1, 2)) if log: print("Found {} peaks".format(n_peaks)) @@ -237,20 +246,17 @@ def _vectfit_xs(energy, ce_xs, mts, rtol=1e-3, atol=1e-5, orders=None, print("Order={}({}/{})".format(order, i, len(orders))) # initial guessed poles poles = np.linspace(s[0], s[-1], order//2) - poles = poles + poles*0.01j + poles += poles*0.01j poles = np.sort(np.append(poles, np.conj(poles))) found_better = False # fitting iteration for i_vf in range(n_vf_iter): if log >= DETAILED_LOGGING: - print("VF iteration {}/{}".format(i_vf+1, n_vf_iter)) + print("VF iteration {}/{}".format(i_vf + 1, n_vf_iter)) # call vf - try: - poles, residues, cf, f_fit, rms = vf.vectfit(f, s, poles, weight) - except: - break + poles, residues, cf, f_fit, rms = vf.vectfit(f, s, poles, weight) # convert real pole to conjugate pairs n_real_poles = 0 @@ -280,9 +286,13 @@ def _vectfit_xs(energy, ce_xs, mts, rtol=1e-3, atol=1e-5, orders=None, maxre, ratio, ratio2 = 0., 1., 1. else: maxre = np.max(relerr[abserr > atol]) - ratio = np.sum((relerr piece_width: raise ValueError('Window spacing cannot be larger than piece spacing.') - # determine curve fit order - if n_cf is None: - warnings.warn("Curvefit order not specified.") - n_cf = 5 - if log: print("Windowing with # windows: {}, spacing: {}, CF order: {}".format( n_win, spacing, n_cf)) @@ -621,7 +625,7 @@ def _windowing(mp_data, rtol=1e-3, atol=1e-5, n_win=None, n_cf=None, win_data = [] for iw in range(n_win): if log >= DETAILED_LOGGING: - print("Processing window {}/{}...".format(iw+1, n_win)) + print("Processing window {}/{}...".format(iw + 1, n_win)) # inner window boundaries inbegin = sqrt(E_min) + spacing * iw @@ -631,11 +635,11 @@ def _windowing(mp_data, rtol=1e-3, atol=1e-5, n_win=None, n_cf=None, if iw == 0 or sqrt(alpha)*inbegin < 4.0: e_start = inbegin**2 else: - e_start = max(E_min, (sqrt(alpha)*inbegin-4.0)**2/alpha) + e_start = max(E_min, (sqrt(alpha)*inbegin - 4.0)**2/alpha) e_end = min(E_max, (sqrt(alpha)*inend + 4.0)**2/alpha) # locate piece and relevant poles - i_piece = min(n_pieces-1, int((inbegin - sqrt(E_min))/piece_width + 0.5)) + i_piece = min(n_pieces - 1, int((inbegin - sqrt(E_min))/piece_width + 0.5)) poles, residues = mp_poles[i_piece], mp_residues[i_piece] n_poles = poles.size @@ -648,11 +652,11 @@ def _windowing(mp_data, rtol=1e-3, atol=1e-5, n_win=None, n_cf=None, xs_ref = vf.evaluate(energy_sqrt, poles, residues*1j) / energy # curve fit matrix - matrix = np.vstack([energy**(0.5*i-1) for i in range(n_cf+1)]).T + matrix = np.vstack([energy**(0.5*i - 1) for i in range(n_cf + 1)]).T # start from 0 poles, initialize pointers to the center nearest pole center_pole_ind = np.argmin((np.fabs(poles.real - incenter))) - lp, rp = center_pole_ind, center_pole_ind + lp = rp = center_pole_ind while True: if log >= DETAILED_LOGGING: print("Trying poles {} to {}".format(lp, rp)) @@ -672,9 +676,9 @@ def _windowing(mp_data, rtol=1e-3, atol=1e-5, n_win=None, n_cf=None, abserr = np.abs(xs_fit + xs_wp - xs_ref) relerr = abserr / xs_ref if not np.any(np.isnan(abserr)): - re = relerr[abserr>atol] + re = relerr[abserr > atol] if re.size == 0 or np.all(re <= rtol) or \ - (re.max() <= 2*rtol and (re>rtol).sum() <= 0.01*relerr.size) or \ + (re.max() <= 2*rtol and (re > rtol).sum() <= 0.01*relerr.size) or \ (iw == 0 and np.all(relerr.mean(axis=1) <= rtol)): # meet tolerances if log >= DETAILED_LOGGING: @@ -688,11 +692,11 @@ def _windowing(mp_data, rtol=1e-3, atol=1e-5, n_win=None, n_cf=None, # try to include one more pole (next center nearest) if rp >= n_poles: - lp = lp - 1 - elif lp <= 0 or poles[rp] - incenter <= incenter - poles[lp-1]: - rp = rp + 1 + lp -= 1 + elif lp <= 0 or poles[rp] - incenter <= incenter - poles[lp - 1]: + rp += 1 else: - lp = lp - 1 + lp -= 1 # save data for this window win_data.append((i_piece, lp, rp, coefs)) @@ -715,8 +719,8 @@ def _windowing(mp_data, rtol=1e-3, atol=1e-5, n_win=None, n_cf=None, ip, lp, rp, coefs = win_data[iw] # adjust indices and change to 1-based for the library format n_prev_poles = sum([poles_unused[i].size for i in range(ip)]) - n_unused = sum([(poles_unused[i]==1).sum() for i in range(ip)]) + \ - (poles_unused[ip][:lp]==1).sum() + n_unused = sum([(poles_unused[i] == 1).sum() for i in range(ip)]) + \ + (poles_unused[ip][:lp] == 1).sum() lp += n_prev_poles - n_unused + 1 rp += n_prev_poles - n_unused windows.append([lp, rp]) @@ -1009,7 +1013,7 @@ class WindowedMultipole(EqualityMixin): return out @classmethod - def from_endf(cls, endf_file, vf_options={}, wmp_options={}): + def from_endf(cls, endf_file, vf_options=None, wmp_options=None): """Generate windowed multipole neutron data from an ENDF evaluation. Parameters @@ -1017,10 +1021,10 @@ class WindowedMultipole(EqualityMixin): endf_file : str Path to ENDF evaluation vf_options : dict - Dictionary of keyword arguments, e.g. {'rtol':0.001, 'log':True}, + Dictionary of keyword arguments, e.g. {'rtol': 0.001, 'log': True}, passed to :func:`openmc.data.multipole._vectfit_nuclide` wmp_options : dict - Dictionary of keyword arguments, e.g. {'search':True, 'log': True}, + Dictionary of keyword arguments, e.g. {'search': True, 'log': True}, passed to :func:`openmc.data.WindowedMultipole.from_multipole` Returns @@ -1031,6 +1035,12 @@ class WindowedMultipole(EqualityMixin): """ + if vf_options is None: + vf_options = {} + + if wmp_options is None: + wmp_options = {} + # generate multipole data from EDNF mp_data = _vectfit_nuclide(endf_file, **vf_options) @@ -1044,11 +1054,11 @@ class WindowedMultipole(EqualityMixin): Parameters ---------- mp_data : dictionary or str - Dictionary or Path to the multipole data + Dictionary or Path to the multipole data stored in a pickle file search : bool, optional Whether to search for optimal window size and curvefit order log : bool or int, optional - Whether to print running logs + Whether to print running logs (use int for verbosity control) **kwargs Keyword arguments passed to :func:`openmc.data.multipole._windowing` @@ -1062,7 +1072,6 @@ class WindowedMultipole(EqualityMixin): if isinstance(mp_data, str): # load multipole data from file - import pickle with open(mp_data, 'rb') as f: mp_data = pickle.load(f) @@ -1076,7 +1085,7 @@ class WindowedMultipole(EqualityMixin): n_poles = sum([p.size for p in mp_data["poles"]]) n_win_min = max(5, n_poles // 20) n_win_max = 2000 if n_poles < 2000 else 8000 - best_wmp, best_metric = None, None + best_wmp = best_metric = None for n_w in np.unique(np.linspace(n_win_min, n_win_max, 20, dtype=int)): for n_cf in range(10, 1, -1): if log: @@ -1163,7 +1172,7 @@ class WindowedMultipole(EqualityMixin): dopp = self.sqrtAWR / sqrtkT broadened_polynomials = _broaden_wmp_polynomials(E, dopp, self.fit_order + 1) - for i_poly in range(self.fit_order+1): + for i_poly in range(self.fit_order + 1): sig_s += (self.curvefit[i_window, i_poly, _FIT_S] * broadened_polynomials[i_poly]) sig_a += (self.curvefit[i_window, i_poly, _FIT_A] @@ -1173,7 +1182,7 @@ class WindowedMultipole(EqualityMixin): * broadened_polynomials[i_poly]) else: temp = invE - for i_poly in range(self.fit_order+1): + for i_poly in range(self.fit_order + 1): sig_s += self.curvefit[i_window, i_poly, _FIT_S] * temp sig_a += self.curvefit[i_window, i_poly, _FIT_A] * temp if self.fissionable: diff --git a/tests/unit_tests/test_data_multipole.py b/tests/unit_tests/test_data_multipole.py index 81853f820a..2a12f937a5 100644 --- a/tests/unit_tests/test_data_multipole.py +++ b/tests/unit_tests/test_data_multipole.py @@ -55,12 +55,15 @@ def test_export_to_hdf5(tmpdir, u235): u235.export_to_hdf5(filename) assert os.path.exists(filename) + @vf_only def test_from_endf(): endf_data = os.environ['OPENMC_ENDF_DATA'] endf_file = os.path.join(endf_data, 'neutrons', 'n-001_H_001.endf') return openmc.data.WindowedMultipole.from_endf( endf_file, wmp_options={"n_win": 400, "n_cf": 3}) + + @vf_only def test_from_endf_search(): endf_data = os.environ['OPENMC_ENDF_DATA']