crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-21 14:35:27 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5
OpenMC/openmc/data/multipole.py
azim-givron 360ec24b41
Implement vector fitting to replace external vectfit package (#3493) 
Co-authored-by: azim_givron <a.givron@naarea.fr>
Co-authored-by: Paul Romano <paul.k.romano@gmail.com>
Co-authored-by: GuySten <guyste@post.bgu.ac.il>
Co-authored-by: GuySten <62616591+GuySten@users.noreply.github.com>
2026-02-11 10:00:18 -06:00

1301 lines
46 KiB
Python
Raw Permalink Blame History

from numbers import Real
from math import exp, erf, pi, sqrt
from copy import deepcopy
import os
import h5py
import pickle
import numpy as np
from scipy.signal import find_peaks
import openmc.checkvalue as cv
from ..exceptions import DataError
from ..mixin import EqualityMixin
from . import WMP_VERSION, WMP_VERSION_MAJOR
from .data import K_BOLTZMANN
from .neutron import IncidentNeutron
from .resonance import ResonanceRange
from .vectfit import vectfit, evaluate
# Constants that determine which value to access
_MP_EA = 0 # Pole
# Residue indices
_MP_RS = 1 # Residue scattering
_MP_RA = 2 # Residue absorption
_MP_RF = 3 # Residue fission
# Polynomial fit indices
_FIT_S = 0 # Scattering
_FIT_A = 1 # Absorption
_FIT_F = 2 # Fission
# Upper temperature limit (K)
TEMPERATURE_LIMIT = 3000
# Logging control
DETAILED_LOGGING = 2
def _faddeeva(z):
r"""Evaluate the complex Faddeeva function.
Technically, the value we want is given by the equation:
.. math::
w(z) = \frac{i}{\pi} \int_{-\infty}^{\infty} \frac{1}{z - t}
\exp(-t^2) \text{d}t
as shown in Equation 63 from Hwang, R. N. "A rigorous pole
representation of multilevel cross sections and its practical
applications." Nuclear Science and Engineering 96.3 (1987): 192-209.
The :func:`scipy.special.wofz` function evaluates
:math:`w(z) = \exp(-z^2) \text{erfc}(-iz)`. These two forms of the Faddeeva
function are related by a transformation.
If we call the integral form :math:`w_\text{int}`, and the function form
:math:`w_\text{fun}`:
.. math::
w_\text{int}(z) =
\begin{cases}
w_\text{fun}(z) & \text{for } \text{Im}(z) > 0\\
-w_\text{fun}(z^*)^* & \text{for } \text{Im}(z) < 0
\end{cases}
Parameters
----------
z : complex
Argument to the Faddeeva function.
Returns
-------
complex
:math:`\frac{i}{\pi} \int_{-\infty}^{\infty} \frac{1}{z - t} \exp(-t^2)
\text{d}t`
"""
from scipy.special import wofz
if np.angle(z) > 0:
return wofz(z)
else:
return -np.conj(wofz(z.conjugate()))
def _broaden_wmp_polynomials(E, dopp, n):
r"""Evaluate Doppler-broadened windowed multipole curvefit.
The curvefit is a polynomial of the form :math:`\frac{a}{E}
+ \frac{b}{\sqrt{E}} + c + d \sqrt{E} + \ldots`
Parameters
----------
E : float
Energy to evaluate at.
dopp : float
sqrt(atomic weight ratio / kT) in units of eV.
n : int
Number of components to the polynomial.
Returns
-------
np.ndarray
The value of each Doppler-broadened curvefit polynomial term.
"""
sqrtE = sqrt(E)
beta = sqrtE * dopp
half_inv_dopp2 = 0.5 / dopp**2
quarter_inv_dopp4 = half_inv_dopp2**2
if beta > 6.0:
# Save time, ERF(6) is 1 to machine precision.
# beta/sqrtpi*exp(-beta**2) is also approximately 1 machine epsilon.
erf_beta = 1.0
exp_m_beta2 = 0.0
else:
erf_beta = erf(beta)
exp_m_beta2 = exp(-beta**2)
# Assume that, for sure, we'll use a second order (1/E, 1/V, const)
# fit, and no less.
factors = np.zeros(n)
factors[0] = erf_beta / E
factors[1] = 1.0 / sqrtE
factors[2] = (factors[0] * (half_inv_dopp2 + E)
+ exp_m_beta2 / (beta * sqrt(pi)))
# Perform recursive broadening of high order components. range(1, n-2)
# replaces a do i = 1, n-3. All indices are reduced by one due to the
# 1-based vs. 0-based indexing.
for i in range(1, n-2):
if i != 1:
factors[i+2] = (-factors[i-2] * (i - 1.0) * i * quarter_inv_dopp4
+ factors[i] * (E + (1.0 + 2.0 * i) * half_inv_dopp2))
else:
factors[i+2] = factors[i]*(E + (1.0 + 2.0 * i) * half_inv_dopp2)
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):
"""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, with shape (number of reactions,
number of energy points)
mts : Iterable of int
Reaction list
rtol : float, optional
Relative error tolerance
atol : float, optional
Absolute error tolerance
orders : Iterable of int, optional
A list of orders (number of poles) to be searched
n_vf_iter : int, optional
Number of maximum VF iterations
log : bool or int, optional
Whether to print running logs (use int for verbosity control)
path_out : str, optional
Path to save the figures to show discrepancies between the original and
fitted cross sections for different reactions
Returns
-------
tuple
(poles, residues)
"""
ne = energy.size
nmt = len(mts)
if ce_xs.shape != (nmt, ne):
raise ValueError('Inconsistent cross section data.')
# construct test data: interpolate xs with finer grids
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)
test_energy[[0, -1]] = energy[[0, -1]] # avoid numerical issue
test_xs_ref = np.zeros((nmt, ne_test))
for i in range(nmt):
test_xs_ref[i] = np.interp(test_energy, energy, ce_xs[i])
if log:
print(f"\tenergy: {energy[0]:.3e} to {energy[-1]:.3e} eV ({ne} points)")
print(f"\terror tolerance: rtol={rtol}, atol={atol}")
# 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
with np.errstate(divide='ignore'):
weight = 1.0/f
# 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):
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])
# detect peaks (resonances) and determine VF order search range
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]))
else:
lowest_order = max(2, 2*n_peaks)
highest_order = max(200, 4*n_peaks)
orders = list(range(lowest_order, highest_order + 1, 2))
if log:
print(f"Found {n_peaks} peaks")
print(f"Fitting orders from {orders[0]} to {orders[-1]}")
# perform VF with increasing orders
found_ideal = False
n_discarded = 0 # for accelation, number of discarded searches
best_quality = best_ratio = -np.inf
for i, order in enumerate(orders):
if log:
print(f"Order={order}({i}/{len(orders)})")
# initial guessed poles
poles_r = np.linspace(s[0], s[-1], order//2)
poles = poles_r + poles_r*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(f"VF iteration {i_vf + 1}/{n_vf_iter}")
# call vf
poles, residues, *_ = vectfit(f, s, poles, weight)
# convert real pole to conjugate pairs
n_real_poles = 0
new_poles = []
for p in poles:
p_r, p_i = np.real(p), np.imag(p)
if (s[0] <= p_r <= s[-1]) and p_i == 0.:
new_poles += [p_r+p_r*0.01j, p_r-p_r*0.01j]
n_real_poles += 1
else:
new_poles += [p]
new_poles = np.array(new_poles)
# re-calculate residues if poles changed
if n_real_poles > 0:
if log >= DETAILED_LOGGING:
print(f" # real poles: {n_real_poles}")
new_poles, residues, *_ = \
vectfit(f, s, new_poles, weight, skip_pole_update=True)
# assess the result on test grid
test_xs = evaluate(test_s, new_poles, residues) / test_energy
abserr = np.abs(test_xs - test_xs_ref)
with np.errstate(invalid='ignore', divide='ignore'):
relerr = abserr / test_xs_ref
if np.any(np.isnan(abserr)):
maxre, ratio, ratio2 = np.inf, -np.inf, -np.inf
elif np.all(abserr <= atol):
maxre, ratio, ratio2 = 0., 1., 1.
else:
maxre = np.max(relerr[abserr > atol])
ratio = np.sum((relerr < rtol) | (abserr < atol)) / relerr.size
ratio2 = np.sum((relerr < 10*rtol) | (abserr < atol)) / relerr.size
# define a metric for choosing the best fitting results
# basically, it is preferred to have more points within accuracy
# tolerance, smaller maximum deviation and fewer poles
#TODO: improve the metric with clearer basis
quality = ratio + ratio2 - min(0.1*maxre, 1) - 0.001*new_poles.size
if np.any(test_xs < -atol):
quality = -np.inf
if log >= DETAILED_LOGGING:
print(f" # poles: {new_poles.size}")
print(f" Max relative error: {maxre * 100:.3f}%")
print(f" Satisfaction: {ratio * 100:.1f}%, {ratio2 * 100:.1f}%")
print(f" Quality: {quality:.2f}")
if quality > best_quality:
if log >= DETAILED_LOGGING:
print(" Best so far!")
found_better = True
best_quality, best_ratio = quality, ratio
best_poles, best_residues = new_poles, residues
best_test_xs, best_relerr = test_xs, relerr
if best_ratio >= 1.0:
if log:
print("Found ideal results. Stop!")
found_ideal = True
break
else:
if log >= DETAILED_LOGGING:
print(" Discarded!")
if found_ideal:
break
# acceleration
if found_better:
n_discarded = 0
else:
if order > max(2*n_peaks, 50) and best_ratio > 0.7:
n_discarded += 1
if n_discarded >= 10 or (n_discarded >= 5 and best_ratio > 0.9):
if log >= DETAILED_LOGGING:
print("Couldn't get better results. Stop!")
break
# merge conjugate poles
real_idx = []
conj_idx = []
found_conj = False
for i, p in enumerate(best_poles):
if found_conj:
found_conj = False
continue
if np.imag(p) == 0.:
real_idx.append(i)
else:
if i < best_poles.size and np.conj(p) == best_poles[i + 1]:
found_conj = True
conj_idx.append(i)
else:
raise RuntimeError("Complex poles are not conjugate!")
if log:
print("Found {} real poles and {} conjugate complex pairs.".format(
len(real_idx), len(conj_idx)))
mp_poles = best_poles[real_idx + conj_idx]
mp_residues = np.concatenate((best_residues[:, real_idx],
best_residues[:, conj_idx]*2), axis=1)/1j
if log:
print(f"Final number of poles: {mp_poles.size}")
if path_out:
if not os.path.exists(path_out):
os.makedirs(path_out)
for i, mt in enumerate(mts):
if not test_xs_ref[i].any():
continue
import matplotlib.pyplot as plt
fig, ax1 = plt.subplots()
lns1 = ax1.loglog(test_energy, test_xs_ref[i], 'g', label="ACE xs")
lns2 = ax1.loglog(test_energy, best_test_xs[i], 'b', label="VF xs")
ax2 = ax1.twinx()
lns3 = ax2.loglog(test_energy, best_relerr[i], 'r',
label="Relative error", alpha=0.5)
lns = lns1 + lns2 + lns3
labels = [l.get_label() for l in lns]
ax1.legend(lns, labels, loc='best')
ax1.set_xlabel('energy (eV)')
ax1.set_ylabel('cross section (b)', color='b')
ax1.tick_params('y', colors='b')
ax2.set_ylabel('relative error', color='r')
ax2.tick_params('y', colors='r')
plt.title(f"MT {mt} vector fitted with {mp_poles.size} poles")
fig.tight_layout()
fig_file = os.path.join(path_out, "{:.0f}-{:.0f}_MT{}.png".format(
energy[0], energy[-1], mt))
plt.savefig(fig_file)
plt.close()
if log:
print(f"Saved figure: {fig_file}")
return (mp_poles, mp_residues)
def vectfit_nuclide(endf_file, njoy_error=5e-4, vf_pieces=None,
log=False, path_out=None, mp_filename=None,
**kwargs):
r"""Generate multipole data for a nuclide from ENDF.
Parameters
----------
endf_file : str
Path to ENDF evaluation
njoy_error : float, optional
Fractional error tolerance for processing point-wise data with NJOY
vf_pieces : integer, optional
Number of equal-in-momentum spaced energy pieces for data fitting
log : bool or int, optional
Whether to print running logs (use int for verbosity control)
path_out : str, optional
Path to write out mutipole data file and vector fitting figures
mp_filename : str, optional
File name to write out multipole data
**kwargs
Keyword arguments passed to :func:`openmc.data.multipole._vectfit_xs`
Returns
-------
mp_data
Dictionary containing necessary multipole data of the nuclide
"""
# ======================================================================
# PREPARE POINT-WISE XS
# make 0K ACE data using njoy
if log:
print(f"Running NJOY to get 0K point-wise data (error={njoy_error})...")
nuc_ce = IncidentNeutron.from_njoy(endf_file, temperatures=[0.0],
error=njoy_error, broadr=False, heatr=False, purr=False)
if log:
print("Parsing cross sections within resolved resonance range...")
# Determine upper energy: the lower of RRR upper bound and first threshold
endf_res = IncidentNeutron.from_endf(endf_file).resonances
if hasattr(endf_res, 'resolved') and \
hasattr(endf_res.resolved, 'energy_max') and \
type(endf_res.resolved) is not ResonanceRange:
E_max = endf_res.resolved.energy_max
elif hasattr(endf_res, 'unresolved') and \
hasattr(endf_res.unresolved, 'energy_min'):
E_max = endf_res.unresolved.energy_min
else:
E_max = nuc_ce.energy['0K'][-1]
E_max_idx = np.searchsorted(nuc_ce.energy['0K'], E_max, side='right') - 1
for mt in nuc_ce.reactions:
if hasattr(nuc_ce.reactions[mt].xs['0K'], '_threshold_idx'):
threshold_idx = nuc_ce.reactions[mt].xs['0K']._threshold_idx
if 0 < threshold_idx < E_max_idx:
E_max_idx = threshold_idx
# parse energy and cross sections
energy = nuc_ce.energy['0K'][:E_max_idx + 1]
E_min, E_max = energy[0], energy[-1]
n_points = energy.size
total_xs = nuc_ce[1].xs['0K'](energy)
elastic_xs = nuc_ce[2].xs['0K'](energy)
try:
absorption_xs = nuc_ce[27].xs['0K'](energy)
except KeyError:
absorption_xs = np.zeros_like(total_xs)
fissionable = False
try:
fission_xs = nuc_ce[18].xs['0K'](energy)
fissionable = True
except KeyError:
pass
# make vectors
if fissionable:
ce_xs = np.vstack((elastic_xs, absorption_xs, fission_xs))
mts = [2, 27, 18]
else:
ce_xs = np.vstack((elastic_xs, absorption_xs))
mts = [2, 27]
if log:
print(f" MTs: {mts}")
print(f" Energy range: {E_min:.3e} to {E_max:.3e} eV ({n_points} points)")
# ======================================================================
# PERFORM VECTOR FITTING
if vf_pieces is None:
# divide into pieces for complex nuclides
peaks, _ = find_peaks(total_xs)
n_peaks = peaks.size
if n_peaks > 200 or n_points > 30000 or n_peaks * n_points > 100*10000:
vf_pieces = max(5, n_peaks // 50, n_points // 2000)
else:
vf_pieces = 1
piece_width = (sqrt(E_max) - sqrt(E_min)) / vf_pieces
alpha = nuc_ce.atomic_weight_ratio/(K_BOLTZMANN*TEMPERATURE_LIMIT)
poles, residues = [], []
# VF piece by piece
for i_piece in range(vf_pieces):
if log:
print(f"Vector fitting piece {i_piece + 1}/{vf_pieces}...")
# start E of this piece
e_bound = (sqrt(E_min) + piece_width*(i_piece-0.5))**2
if i_piece == 0 or sqrt(alpha*e_bound) < 4.0:
e_start = E_min
e_start_idx = 0
else:
e_start = max(E_min, (sqrt(alpha*e_bound) - 4.0)**2/alpha)
e_start_idx = np.searchsorted(energy, e_start, side='right') - 1
# end E of this piece
e_bound = (sqrt(E_min) + piece_width*(i_piece + 1))**2
e_end = min(E_max, (sqrt(alpha*e_bound) + 4.0)**2/alpha)
e_end_idx = np.searchsorted(energy, e_end, side='left') + 1
e_idx = range(e_start_idx, min(e_end_idx + 1, n_points))
p, r = _vectfit_xs(energy[e_idx], ce_xs[:, e_idx], mts, log=log,
path_out=path_out, **kwargs)
poles.append(p)
residues.append(r)
# collect multipole data into a dictionary
mp_data = {"name": nuc_ce.name,
"AWR": nuc_ce.atomic_weight_ratio,
"E_min": E_min,
"E_max": E_max,
"poles": poles,
"residues": residues}
# dump multipole data to file
if path_out:
if not os.path.exists(path_out):
os.makedirs(path_out)
if not mp_filename:
mp_filename = f"{nuc_ce.name}_mp.pickle"
mp_filename = os.path.join(path_out, mp_filename)
with open(mp_filename, 'wb') as f:
pickle.dump(mp_data, f)
if log:
print(f"Dumped multipole data to file: {mp_filename}")
return mp_data
def _windowing(mp_data, n_cf, rtol=1e-3, atol=1e-5, n_win=None, spacing=None,
log=False):
"""Generate windowed multipole library from multipole data with specific
settings of window size, curve fit order, etc.
Parameters
----------
mp_data : dict
Multipole data
n_cf : int
Curve fitting order
rtol : float, optional
Maximum relative error tolerance
atol : float, optional
Minimum absolute error tolerance
n_win : int, optional
Number of equal-in-mementum spaced energy windows
spacing : float, optional
Inner window spacing (sqrt energy space)
log : bool or int, optional
Whether to print running logs (use int for verbosity control)
Returns
-------
openmc.data.WindowedMultipole
Resonant cross sections represented in the windowed multipole
format.
"""
# unpack multipole data
name = mp_data["name"]
awr = mp_data["AWR"]
E_min = mp_data["E_min"]
E_max = mp_data["E_max"]
mp_poles = mp_data["poles"]
mp_residues = mp_data["residues"]
n_pieces = len(mp_poles)
piece_width = (sqrt(E_max) - sqrt(E_min)) / n_pieces
alpha = awr / (K_BOLTZMANN*TEMPERATURE_LIMIT)
# determine window size
if n_win is None:
if spacing is not None:
# ensure the windows are within the multipole energy range
n_win = int((sqrt(E_max) - sqrt(E_min)) / spacing)
E_max = (sqrt(E_min) + n_win*spacing)**2
else:
n_win = 1000
# inner window size
spacing = (sqrt(E_max) - sqrt(E_min)) / n_win
# make sure inner window size is smaller than energy piece size
if spacing > piece_width:
raise ValueError('Window spacing cannot be larger than piece spacing.')
if log:
print("Windowing:")
print(f" config: # windows={n_win}, spacing={spacing}, CF order={n_cf}")
print(f" error tolerance: rtol={rtol}, atol={atol}")
# sort poles (and residues) by the real component of the pole
for ip in range(n_pieces):
indices = mp_poles[ip].argsort()
mp_poles[ip] = mp_poles[ip][indices]
mp_residues[ip] = mp_residues[ip][:, indices]
# initialize an array to record whether each pole is used or not
poles_unused = [np.ones_like(p, dtype=int) for p in mp_poles]
# optimize the windows: the goal is to find the least set of significant
# consecutive poles and curve fit coefficients to reproduce cross section
win_data = []
for iw in range(n_win):
if log >= DETAILED_LOGGING:
print(f"Processing window {iw + 1}/{n_win}...")
# inner window boundaries
inbegin = sqrt(E_min) + spacing * iw
inend = inbegin + spacing
incenter = (inbegin + inend) / 2.0
# extend window energy range for Doppler broadening
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_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))
poles, residues = mp_poles[i_piece], mp_residues[i_piece]
n_poles = poles.size
# generate energy points for fitting: equally spaced in momentum
n_points = min(max(100, int((e_end - e_start)*4)), 10000)
energy_sqrt = np.linspace(np.sqrt(e_start), np.sqrt(e_end), n_points)
energy = energy_sqrt**2
# reference xs from multipole form, note the residue terms in the
# multipole and vector fitting representations differ by a 1j
xs_ref = 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
# 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
while True:
if log >= DETAILED_LOGGING:
print(f"Trying poles {lp} to {rp}")
# calculate the cross sections contributed by the windowed poles
if rp > lp:
xs_wp = evaluate(energy_sqrt, poles[lp:rp],
residues[:, lp:rp]*1j) / energy
else:
xs_wp = np.zeros_like(xs_ref)
# do least square curve fit on the remains
coefs = np.linalg.lstsq(matrix, (xs_ref - xs_wp).T, rcond=None)[0]
xs_fit = (matrix @ coefs).T
# assess the result
abserr = np.abs(xs_fit + xs_wp - xs_ref)
with np.errstate(invalid='ignore', divide='ignore'):
relerr = abserr / xs_ref
if not np.any(np.isnan(abserr)):
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 \
(iw == 0 and np.all(relerr.mean(axis=1) <= rtol)):
# meet tolerances
if log >= DETAILED_LOGGING:
print("Accuracy satisfied.")
break
# we expect pure curvefit will succeed for the first window
# TODO: find the energy boundary below which no poles are allowed
if iw == 0:
raise RuntimeError('Pure curvefit failed for the first window!')
# try to include one more pole (next center nearest)
if rp >= n_poles:
lp -= 1
elif lp <= 0 or poles[rp] - incenter <= incenter - poles[lp - 1]:
rp += 1
else:
lp -= 1
# save data for this window
win_data.append((i_piece, lp, rp, coefs))
# mark the windowed poles as used poles
poles_unused[i_piece][lp:rp] = 0
# flatten and shrink by removing unused poles
data = [] # used poles and residues
for ip in range(n_pieces):
used = (poles_unused[ip] == 0)
# stack poles and residues for library format
data.append(np.vstack([mp_poles[ip][used], mp_residues[ip][:, used]]).T)
# stack poles/residues in sequence vertically
data = np.vstack(data)
# new start/end pole indices
windows = []
curvefit = []
for iw in range(n_win):
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()
lp += n_prev_poles - n_unused + 1
rp += n_prev_poles - n_unused
windows.append([lp, rp])
curvefit.append(coefs)
# construct the WindowedMultipole object
wmp = WindowedMultipole(name)
wmp.spacing = spacing
wmp.sqrtAWR = sqrt(awr)
wmp.E_min = E_min
wmp.E_max = E_max
wmp.data = data
wmp.windows = np.asarray(windows)
wmp.curvefit = np.asarray(curvefit)
# TODO: check if Doppler brodening of the polynomial curvefit is negligible
wmp.broaden_poly = np.ones((n_win,), dtype=bool)
return wmp
class WindowedMultipole(EqualityMixin):
"""Resonant cross sections represented in the windowed multipole format.
Parameters
----------
name : str
Name of the nuclide using the GNDS naming convention
Attributes
----------
name : str
Name of the nuclide using the GNDS naming convention
spacing : float
The width of each window in sqrt(E)-space. For example, the frst window
will end at (sqrt(E_min) + spacing)**2 and the second window at
(sqrt(E_min) + 2*spacing)**2.
sqrtAWR : float
Square root of the atomic weight ratio of the target nuclide.
E_min : float
Lowest energy in eV the library is valid for.
E_max : float
Highest energy in eV the library is valid for.
data : np.ndarray
A 2D array of complex poles and residues. data[i, 0] gives the energy
at which pole i is located. data[i, 1:] gives the residues associated
with the i-th pole. There are 3 residues, one each for the scattering,
absorption, and fission channels.
windows : np.ndarray
A 2D array of Integral values. windows[i, 0] - 1 is the index of the
first pole in window i. windows[i, 1] - 1 is the index of the last pole
in window i.
broaden_poly : np.ndarray
A 1D array of boolean values indicating whether or not the polynomial
curvefit in that window should be Doppler broadened.
curvefit : np.ndarray
A 3D array of Real curvefit polynomial coefficients. curvefit[i, 0, :]
gives coefficients for the scattering cross section in window i.
curvefit[i, 1, :] gives absorption coefficients and curvefit[i, 2, :]
gives fission coefficients. The polynomial terms are increasing powers
of sqrt(E) starting with 1/E e.g:
a/E + b/sqrt(E) + c + d sqrt(E) + ...
"""
def __init__(self, name):
self.name = name
self.spacing = None
self.sqrtAWR = None
self.E_min = None
self.E_max = None
self.data = None
self.windows = None
self.broaden_poly = None
self.curvefit = None
@property
def name(self):
return self._name
@name.setter
def name(self, name):
cv.check_type('name', name, str)
self._name = name
@property
def fit_order(self):
return self.curvefit.shape[1] - 1
@property
def fissionable(self):
return self.data.shape[1] == 4
@property
def n_poles(self):
return self.data.shape[0]
@property
def n_windows(self):
return self.windows.shape[0]
@property
def poles_per_window(self):
return (self.windows[:, 1] - self.windows[:, 0] + 1).mean()
@property
def spacing(self):
return self._spacing
@spacing.setter
def spacing(self, spacing):
if spacing is not None:
cv.check_type('spacing', spacing, Real)
cv.check_greater_than('spacing', spacing, 0.0, equality=False)
self._spacing = spacing
@property
def sqrtAWR(self):
return self._sqrtAWR
@sqrtAWR.setter
def sqrtAWR(self, sqrtAWR):
if sqrtAWR is not None:
cv.check_type('sqrtAWR', sqrtAWR, Real)
cv.check_greater_than('sqrtAWR', sqrtAWR, 0.0, equality=False)
self._sqrtAWR = sqrtAWR
@property
def E_min(self):
return self._E_min
@E_min.setter
def E_min(self, E_min):
if E_min is not None:
cv.check_type('E_min', E_min, Real)
cv.check_greater_than('E_min', E_min, 0.0, equality=True)
self._E_min = E_min
@property
def E_max(self):
return self._E_max
@E_max.setter
def E_max(self, E_max):
if E_max is not None:
cv.check_type('E_max', E_max, Real)
cv.check_greater_than('E_max', E_max, 0.0, equality=False)
self._E_max = E_max
@property
def data(self):
return self._data
@data.setter
def data(self, data):
if data is not None:
cv.check_type('data', data, np.ndarray)
if len(data.shape) != 2:
raise ValueError('Multipole data arrays must be 2D')
if data.shape[1] not in (3, 4):
raise ValueError(
'data.shape[1] must be 3 or 4. One value for the pole.'
' One each for the scattering and absorption residues. '
'Possibly one more for a fission residue.')
if not np.issubdtype(data.dtype, np.complexfloating):
raise TypeError('Multipole data arrays must be complex dtype')
self._data = data
@property
def windows(self):
return self._windows
@windows.setter
def windows(self, windows):
if windows is not None:
cv.check_type('windows', windows, np.ndarray)
if len(windows.shape) != 2:
raise ValueError('Multipole windows arrays must be 2D')
if not np.issubdtype(windows.dtype, np.integer):
raise TypeError('Multipole windows arrays must be integer'
' dtype')
self._windows = windows
@property
def broaden_poly(self):
return self._broaden_poly
@broaden_poly.setter
def broaden_poly(self, broaden_poly):
if broaden_poly is not None:
cv.check_type('broaden_poly', broaden_poly, np.ndarray)
if len(broaden_poly.shape) != 1:
raise ValueError('Multipole broaden_poly arrays must be 1D')
if not np.issubdtype(broaden_poly.dtype, np.bool_):
raise TypeError('Multipole broaden_poly arrays must be boolean'
' dtype')
self._broaden_poly = broaden_poly
@property
def curvefit(self):
return self._curvefit
@curvefit.setter
def curvefit(self, curvefit):
if curvefit is not None:
cv.check_type('curvefit', curvefit, np.ndarray)
if len(curvefit.shape) != 3:
raise ValueError('Multipole curvefit arrays must be 3D')
if curvefit.shape[2] not in (2, 3): # sig_s, sig_a (maybe sig_f)
raise ValueError('The third dimension of multipole curvefit'
' arrays must have a length of 2 or 3')
if not np.issubdtype(curvefit.dtype, np.floating):
raise TypeError('Multipole curvefit arrays must be float dtype')
self._curvefit = curvefit
@classmethod
def from_hdf5(cls, group_or_filename):
"""Construct a WindowedMultipole object from an HDF5 group or file.
Parameters
----------
group_or_filename : h5py.Group or str
HDF5 group containing multipole data. If given as a string, it is
assumed to be the filename for the HDF5 file, and the first group is
used to read from.
Returns
-------
openmc.data.WindowedMultipole
Resonant cross sections represented in the windowed multipole
format.
"""
if isinstance(group_or_filename, h5py.Group):
group = group_or_filename
need_to_close = False
else:
h5file = h5py.File(str(group_or_filename), 'r')
need_to_close = True
# Make sure version matches
if 'version' in h5file.attrs:
major, minor = h5file.attrs['version']
if major != WMP_VERSION_MAJOR:
raise DataError(
'WMP data format uses version {}. {} whereas your '
'installation of the OpenMC Python API expects version '
'{}.x.'.format(major, minor, WMP_VERSION_MAJOR))
else:
raise DataError(
'WMP data does not indicate a version. Your installation of '
'the OpenMC Python API expects version {}.x data.'
.format(WMP_VERSION_MAJOR))
group = list(h5file.values())[0]
name = group.name[1:]
out = cls(name)
# Read scalars.
out.spacing = group['spacing'][()]
out.sqrtAWR = group['sqrtAWR'][()]
out.E_min = group['E_min'][()]
out.E_max = group['E_max'][()]
# Read arrays.
err = "WMP '{}' array shape is not consistent with the '{}' array shape"
out.data = group['data'][()]
out.windows = group['windows'][()]
out.broaden_poly = group['broaden_poly'][...].astype(bool)
if out.broaden_poly.shape[0] != out.windows.shape[0]:
raise ValueError(err.format('broaden_poly', 'windows'))
out.curvefit = group['curvefit'][()]
if out.curvefit.shape[0] != out.windows.shape[0]:
raise ValueError(err.format('curvefit', 'windows'))
# _broaden_wmp_polynomials assumes the curve fit has at least 3 terms.
if out.fit_order < 2:
raise ValueError("Windowed multipole is only supported for "
"curvefits with 3 or more terms.")
# If HDF5 file was opened here, make sure it gets closed
if need_to_close:
h5file.close()
return out
@classmethod
def from_endf(cls, endf_file, log=False, vf_options=None, wmp_options=None):
"""Generate windowed multipole neutron data from an ENDF evaluation.
.. versionadded:: 0.12.1
Parameters
----------
endf_file : str
Path to ENDF evaluation
log : bool or int, optional
Whether to print running logs (use int for verbosity control)
vf_options : dict, optional
Dictionary of keyword arguments, e.g. {'njoy_error': 0.001},
passed to :func:`openmc.data.multipole.vectfit_nuclide`
wmp_options : dict, optional
Dictionary of keyword arguments, e.g. {'search': True, 'rtol': 0.01},
passed to :func:`openmc.data.WindowedMultipole.from_multipole`
Returns
-------
openmc.data.WindowedMultipole
Resonant cross sections represented in the windowed multipole
format.
"""
if vf_options is None:
vf_options = {}
if wmp_options is None:
wmp_options = {}
if log:
vf_options.update(log=log)
wmp_options.update(log=log)
# generate multipole data from EDNF
mp_data = vectfit_nuclide(endf_file, **vf_options)
# windowing
return cls.from_multipole(mp_data, **wmp_options)
@classmethod
def from_multipole(
cls,
mp_data,
search=None,
log=False,
search_n_win=20,
search_cf_orders=None,
**kwargs,
):
"""Generate windowed multipole neutron data from multipole data.
Parameters
----------
mp_data : dictionary or str
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.
Defaults to True if no windowing parameters are specified.
log : bool or int, optional
Whether to print running logs (use int for verbosity control)
search_n_win : int, optional
Number of window sizes to consider in the search grid when
``search`` is True.
search_cf_orders : iterable of int, optional
Curve-fit orders to consider in the search grid when ``search`` is
True. Defaults to integers from 10 down to 2.
**kwargs
Keyword arguments passed to :func:`openmc.data.multipole._windowing`.
Returns
-------
openmc.data.WindowedMultipole
Resonant cross sections represented in the windowed multipole
format.
"""
if isinstance(mp_data, str):
# load multipole data from file
with open(mp_data, 'rb') as f:
mp_data = pickle.load(f)
if search is None:
if 'n_cf' in kwargs and ('n_win' in kwargs or 'spacing' in kwargs):
search = False
else:
search = True
# windowing with specific options
if not search:
# set default value for curvefit order if not specified
if 'n_cf' not in kwargs:
kwargs.update(n_cf=5)
return _windowing(mp_data, log=log, **kwargs)
# search optimal WMP from a range of window sizes and CF orders
if log:
print("Start searching ...")
if search_cf_orders is None:
search_cf_orders = range(10, 1, -1)
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
for n_w in np.unique(
np.linspace(n_win_min, n_win_max, search_n_win, dtype=int)
):
for n_cf in search_cf_orders:
if log:
print(f"Testing N_win={n_w} N_cf={n_cf}")
# update arguments dictionary
kwargs.update(n_win=n_w, n_cf=n_cf)
# windowing
try:
wmp = _windowing(mp_data, log=log, **kwargs)
except Exception as e:
if log:
print('Failed: ' + str(e))
break
# select wmp library with metric:
# - performance: average # used poles per window and CF order
# - memory: # windows
metric = -(wmp.poles_per_window * 10. + wmp.fit_order * 1. +
wmp.n_windows * 0.01)
if best_wmp is None or metric > best_metric:
if log:
print("Best library so far.")
best_wmp = deepcopy(wmp)
best_metric = metric
# return the best wmp library
if log:
print("Final library: {} poles, {} windows, {:.2g} poles per window, "
"{} CF order".format(best_wmp.n_poles, best_wmp.n_windows,
best_wmp.poles_per_window, best_wmp.fit_order))
return best_wmp
def _evaluate(self, E, T):
"""Compute scattering, absorption, and fission cross sections.
Parameters
----------
E : Real
Energy of the incident neutron in eV.
T : Real
Temperature of the target in K.
Returns
-------
3-tuple of Real
Scattering, absorption, and fission microscopic cross sections
at the given energy and temperature.
"""
if E < self.E_min: return (0, 0, 0)
if E > self.E_max: return (0, 0, 0)
# ======================================================================
# Bookkeeping
# Define some frequently used variables.
sqrtkT = sqrt(K_BOLTZMANN * T)
sqrtE = sqrt(E)
invE = 1.0 / E
# Locate us. The i_window calc omits a + 1 present from the legacy
# Fortran version of OpenMC because of the 1-based vs. 0-based
# indexing. Similarly startw needs to be decreased by 1. endw does
# not need to be decreased because range(startw, endw) does not include
# endw.
i_window = min(self.n_windows - 1,
int(np.floor((sqrtE - sqrt(self.E_min)) / self.spacing)))
startw = self.windows[i_window, 0] - 1
endw = self.windows[i_window, 1]
# Initialize the ouptut cross sections.
sig_s = 0.0
sig_a = 0.0
sig_f = 0.0
# ======================================================================
# Add the contribution from the curvefit polynomial.
if sqrtkT != 0 and self.broaden_poly[i_window]:
# Broaden the curvefit.
dopp = self.sqrtAWR / sqrtkT
broadened_polynomials = _broaden_wmp_polynomials(E, dopp,
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]
* broadened_polynomials[i_poly])
if self.fissionable:
sig_f += (self.curvefit[i_window, i_poly, _FIT_F]
* broadened_polynomials[i_poly])
else:
temp = invE
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:
sig_f += self.curvefit[i_window, i_poly, _FIT_F] * temp
temp *= sqrtE
# ======================================================================
# Add the contribution from the poles in this window.
if sqrtkT == 0.0:
# If at 0K, use asymptotic form.
for i_pole in range(startw, endw):
psi_chi = -1j / (self.data[i_pole, _MP_EA] - sqrtE)
c_temp = psi_chi / E
sig_s += (self.data[i_pole, _MP_RS] * c_temp).real
sig_a += (self.data[i_pole, _MP_RA] * c_temp).real
if self.fissionable:
sig_f += (self.data[i_pole, _MP_RF] * c_temp).real
else:
# At temperature, use Faddeeva function-based form.
dopp = self.sqrtAWR / sqrtkT
for i_pole in range(startw, endw):
Z = (sqrtE - self.data[i_pole, _MP_EA]) * dopp
w_val = _faddeeva(Z) * dopp * invE * sqrt(pi)
sig_s += (self.data[i_pole, _MP_RS] * w_val).real
sig_a += (self.data[i_pole, _MP_RA] * w_val).real
if self.fissionable:
sig_f += (self.data[i_pole, _MP_RF] * w_val).real
return sig_s, sig_a, sig_f
def __call__(self, E, T):
"""Compute scattering, absorption, and fission cross sections.
Parameters
----------
E : Real or Iterable of Real
Energy of the incident neutron in eV.
T : Real
Temperature of the target in K.
Returns
-------
3-tuple of Real or 3-tuple of numpy.ndarray
Scattering, absorption, and fission microscopic cross sections
at the given energy and temperature.
"""
fun = np.vectorize(lambda x: self._evaluate(x, T))
return fun(E)
def export_to_hdf5(self, path, mode='a', libver='earliest'):
"""Export windowed multipole data to an HDF5 file.
Parameters
----------
path : str
Path to write HDF5 file to
mode : {'r+', 'w', 'x', 'a'}
Mode that is used to open the HDF5 file. This is the second argument
to the :class:`h5py.File` constructor.
libver : {'earliest', 'latest'}
Compatibility mode for the HDF5 file. 'latest' will produce files
that are less backwards compatible but have performance benefits.
"""
# Open file and write version.
with h5py.File(str(path), mode, libver=libver) as f:
f.attrs['filetype'] = np.bytes_('data_wmp')
f.attrs['version'] = np.array(WMP_VERSION)
g = f.create_group(self.name)
# Write scalars.
g.create_dataset('spacing', data=np.array(self.spacing))
g.create_dataset('sqrtAWR', data=np.array(self.sqrtAWR))
g.create_dataset('E_min', data=np.array(self.E_min))
g.create_dataset('E_max', data=np.array(self.E_max))
# Write arrays.
g.create_dataset('data', data=self.data)
g.create_dataset('windows', data=self.windows)
g.create_dataset('broaden_poly',
data=self.broaden_poly.astype(np.int8))
g.create_dataset('curvefit', data=self.curvefit)
Reference in a new issue View git blame Copy permalink