1from collections import OrderedDict 2from collections.abc import Mapping, Callable 3from copy import deepcopy 4from io import StringIO 5from math import pi 6from numbers import Integral, Real 7import os 8 9import h5py 10import numpy as np 11import pandas as pd 12from scipy.interpolate import CubicSpline 13 14import openmc.checkvalue as cv 15from openmc.mixin import EqualityMixin 16from . import HDF5_VERSION 17from .ace import Table, get_metadata, get_table 18from .data import ATOMIC_SYMBOL, EV_PER_MEV 19from .endf import Evaluation, get_head_record, get_tab1_record, get_list_record 20from .function import Tabulated1D 21 22 23# Constants 24MASS_ELECTRON_EV = 0.5109989461e6 # Electron mass energy 25PLANCK_C = 1.2398419739062977e4 # Planck's constant times c in eV-Angstroms 26FINE_STRUCTURE = 137.035999139 # Inverse fine structure constant 27CM_PER_ANGSTROM = 1.0e-8 28# classical electron radius in cm 29R0 = CM_PER_ANGSTROM * PLANCK_C / (2.0 * pi * FINE_STRUCTURE * MASS_ELECTRON_EV) 30 31# Electron subshell labels 32_SUBSHELLS = [None, 'K', 'L1', 'L2', 'L3', 'M1', 'M2', 'M3', 'M4', 'M5', 33 'N1', 'N2', 'N3', 'N4', 'N5', 'N6', 'N7', 'O1', 'O2', 'O3', 34 'O4', 'O5', 'O6', 'O7', 'O8', 'O9', 'P1', 'P2', 'P3', 'P4', 35 'P5', 'P6', 'P7', 'P8', 'P9', 'P10', 'P11', 'Q1', 'Q2', 'Q3'] 36 37_REACTION_NAME = { 38 501: ('Total photon interaction', 'total'), 39 502: ('Photon coherent scattering', 'coherent'), 40 504: ('Photon incoherent scattering', 'incoherent'), 41 515: ('Pair production, electron field', 'pair_production_electron'), 42 516: ('Total pair production', 'pair_production_total'), 43 517: ('Pair production, nuclear field', 'pair_production_nuclear'), 44 522: ('Photoelectric absorption', 'photoelectric'), 45 525: ('Heating', 'heating'), 46 526: ('Electro-atomic scattering', 'electro_atomic_scat'), 47 527: ('Electro-atomic bremsstrahlung', 'electro_atomic_brem'), 48 528: ('Electro-atomic excitation', 'electro_atomic_excit'), 49 534: ('K (1s1/2) subshell photoelectric', 'K'), 50 535: ('L1 (2s1/2) subshell photoelectric', 'L1'), 51 536: ('L2 (2p1/2) subshell photoelectric', 'L2'), 52 537: ('L3 (2p3/2) subshell photoelectric', 'L3'), 53 538: ('M1 (3s1/2) subshell photoelectric', 'M1'), 54 539: ('M2 (3p1/2) subshell photoelectric', 'M2'), 55 540: ('M3 (3p3/2) subshell photoelectric', 'M3'), 56 541: ('M4 (3d3/2) subshell photoelectric', 'M4'), 57 542: ('M5 (3d5/2) subshell photoelectric', 'M5'), 58 543: ('N1 (4s1/2) subshell photoelectric', 'N1'), 59 544: ('N2 (4p1/2) subshell photoelectric', 'N2'), 60 545: ('N3 (4p3/2) subshell photoelectric', 'N3'), 61 546: ('N4 (4d3/2) subshell photoelectric', 'N4'), 62 547: ('N5 (4d5/2) subshell photoelectric', 'N5'), 63 548: ('N6 (4f5/2) subshell photoelectric', 'N6'), 64 549: ('N7 (4f7/2) subshell photoelectric', 'N7'), 65 550: ('O1 (5s1/2) subshell photoelectric', 'O1'), 66 551: ('O2 (5p1/2) subshell photoelectric', 'O2'), 67 552: ('O3 (5p3/2) subshell photoelectric', 'O3'), 68 553: ('O4 (5d3/2) subshell photoelectric', 'O4'), 69 554: ('O5 (5d5/2) subshell photoelectric', 'O5'), 70 555: ('O6 (5f5/2) subshell photoelectric', 'O6'), 71 556: ('O7 (5f7/2) subshell photoelectric', 'O7'), 72 557: ('O8 (5g7/2) subshell photoelectric', 'O8'), 73 558: ('O9 (5g9/2) subshell photoelectric', 'O9'), 74 559: ('P1 (6s1/2) subshell photoelectric', 'P1'), 75 560: ('P2 (6p1/2) subshell photoelectric', 'P2'), 76 561: ('P3 (6p3/2) subshell photoelectric', 'P3'), 77 562: ('P4 (6d3/2) subshell photoelectric', 'P4'), 78 563: ('P5 (6d5/2) subshell photoelectric', 'P5'), 79 564: ('P6 (6f5/2) subshell photoelectric', 'P6'), 80 565: ('P7 (6f7/2) subshell photoelectric', 'P7'), 81 566: ('P8 (6g7/2) subshell photoelectric', 'P8'), 82 567: ('P9 (6g9/2) subshell photoelectric', 'P9'), 83 568: ('P10 (6h9/2) subshell photoelectric', 'P10'), 84 569: ('P11 (6h11/2) subshell photoelectric', 'P11'), 85 570: ('Q1 (7s1/2) subshell photoelectric', 'Q1'), 86 571: ('Q2 (7p1/2) subshell photoelectric', 'Q2'), 87 572: ('Q3 (7p3/2) subshell photoelectric', 'Q3') 88} 89 90# Compton profiles are read from a pre-generated HDF5 file when they are first 91# needed. The dictionary stores an array of electron momentum values (at which 92# the profiles are tabulated) with the key 'pz' and the profile for each element 93# is a 2D array with shape (n_shells, n_momentum_values) stored on the key Z 94_COMPTON_PROFILES = {} 95 96# Scaled bremsstrahlung DCSs are read from a data file provided by Selzter and 97# Berger when they are first needed. The dictionary stores an array of n 98# incident electron kinetic energies with key 'electron_energies', an array of 99# k reduced photon energies with key 'photon_energies', and the cross sections 100# for each element are in a 2D array with shape (n, k) stored on the key 'Z'. 101# It also stores data used for calculating the density effect correction and 102# stopping power, namely, the mean excitation energy with the key 'I', number 103# of electrons per subshell with the key 'num_electrons', and binding energies 104# with the key 'ionization_energy'. 105_BREMSSTRAHLUNG = {} 106 107 108class AtomicRelaxation(EqualityMixin): 109 """Atomic relaxation data. 110 111 This class stores the binding energy, number of electrons, and electron 112 transitions possible from ioniziation for each electron subshell of an 113 atom. All of the data originates from an ENDF-6 atomic relaxation 114 sub-library (NSUB=6). Instances of this class are not normally instantiated 115 directly but rather created using the factory method 116 :math:`AtomicRelaxation.from_endf`. 117 118 Parameters 119 ---------- 120 binding_energy : dict 121 Dictionary indicating the binding energy in eV (values) for given 122 subshells (keys). The subshells should be given as strings, e.g., 'K', 123 'L1', 'L2', etc. 124 num_electrons : dict 125 Dictionary indicating the number of electrons in a subshell when neutral 126 (values) for given subshells (keys). The subshells should be given as 127 strings, e.g., 'K', 'L1', 'L2', etc. 128 transitions : pandas.DataFrame 129 Dictionary indicating allowed transitions and their probabilities 130 (values) for given subshells (keys). The subshells should be given as 131 strings, e.g., 'K', 'L1', 'L2', etc. The transitions are represented as 132 a DataFrame with columns indicating the secondary and tertiary subshell, 133 the energy of the transition in eV, and the fractional probability of 134 the transition. 135 136 Attributes 137 ---------- 138 binding_energy : dict 139 Dictionary indicating the binding energy in eV (values) for given 140 subshells (keys). The subshells should be given as strings, e.g., 'K', 141 'L1', 'L2', etc. 142 num_electrons : dict 143 Dictionary indicating the number of electrons in a subshell when neutral 144 (values) for given subshells (keys). The subshells should be given as 145 strings, e.g., 'K', 'L1', 'L2', etc. 146 transitions : pandas.DataFrame 147 Dictionary indicating allowed transitions and their probabilities 148 (values) for given subshells (keys). The subshells should be given as 149 strings, e.g., 'K', 'L1', 'L2', etc. The transitions are represented as 150 a DataFrame with columns indicating the secondary and tertiary subshell, 151 the energy of the transition in eV, and the fractional probability of 152 the transition. 153 154 See Also 155 -------- 156 IncidentPhoton 157 158 """ 159 def __init__(self, binding_energy, num_electrons, transitions): 160 self.binding_energy = binding_energy 161 self.num_electrons = num_electrons 162 self.transitions = transitions 163 self._e_fluorescence = {} 164 165 @property 166 def binding_energy(self): 167 return self._binding_energy 168 169 @property 170 def num_electrons(self): 171 return self._num_electrons 172 173 @property 174 def subshells(self): 175 return list(sorted(self.binding_energy.keys())) 176 177 @property 178 def transitions(self): 179 return self._transitions 180 181 @binding_energy.setter 182 def binding_energy(self, binding_energy): 183 cv.check_type('binding energies', binding_energy, Mapping) 184 for subshell, energy in binding_energy.items(): 185 cv.check_value('subshell', subshell, _SUBSHELLS) 186 cv.check_type('binding energy', energy, Real) 187 cv.check_greater_than('binding energy', energy, 0.0, True) 188 self._binding_energy = binding_energy 189 190 @num_electrons.setter 191 def num_electrons(self, num_electrons): 192 cv.check_type('number of electrons', num_electrons, Mapping) 193 for subshell, num in num_electrons.items(): 194 cv.check_value('subshell', subshell, _SUBSHELLS) 195 cv.check_type('number of electrons', num, Real) 196 cv.check_greater_than('number of electrons', num, 0.0, True) 197 self._num_electrons = num_electrons 198 199 @transitions.setter 200 def transitions(self, transitions): 201 cv.check_type('transitions', transitions, Mapping) 202 for subshell, df in transitions.items(): 203 cv.check_value('subshell', subshell, _SUBSHELLS) 204 cv.check_type('transitions', df, pd.DataFrame) 205 self._transitions = transitions 206 207 @classmethod 208 def from_ace(cls, ace): 209 """Generate atomic relaxation data from an ACE file 210 211 Parameters 212 ---------- 213 ace : openmc.data.ace.Table 214 ACE table to read from 215 216 Returns 217 ------- 218 openmc.data.AtomicRelaxation 219 Atomic relaxation data 220 221 """ 222 # Create data dictionaries 223 binding_energy = {} 224 num_electrons = {} 225 transitions = {} 226 227 # Get shell designators 228 n = ace.nxs[7] 229 idx = ace.jxs[11] 230 shells = [_SUBSHELLS[int(i)] for i in ace.xss[idx : idx+n]] 231 232 # Get number of electrons for each shell 233 idx = ace.jxs[12] 234 for shell, num in zip(shells, ace.xss[idx : idx+n]): 235 num_electrons[shell] = num 236 237 # Get binding energy for each shell 238 idx = ace.jxs[13] 239 for shell, e in zip(shells, ace.xss[idx : idx+n]): 240 binding_energy[shell] = e*EV_PER_MEV 241 242 # Get transition table 243 columns = ['secondary', 'tertiary', 'energy (eV)', 'probability'] 244 idx = ace.jxs[18] 245 for i, subi in enumerate(shells): 246 n_transitions = int(ace.xss[ace.jxs[15] + i]) 247 if n_transitions > 0: 248 records = [] 249 for j in range(n_transitions): 250 subj = _SUBSHELLS[int(ace.xss[idx])] 251 subk = _SUBSHELLS[int(ace.xss[idx + 1])] 252 etr = ace.xss[idx + 2]*EV_PER_MEV 253 if j == 0: 254 ftr = ace.xss[idx + 3] 255 else: 256 ftr = ace.xss[idx + 3] - ace.xss[idx - 1] 257 records.append((subj, subk, etr, ftr)) 258 idx += 4 259 260 # Create dataframe for transitions 261 transitions[subi] = pd.DataFrame.from_records( 262 records, columns=columns) 263 264 return cls(binding_energy, num_electrons, transitions) 265 266 @classmethod 267 def from_endf(cls, ev_or_filename): 268 """Generate atomic relaxation data from an ENDF evaluation 269 270 Parameters 271 ---------- 272 ev_or_filename : str or openmc.data.endf.Evaluation 273 ENDF atomic relaxation evaluation to read from. If given as a 274 string, it is assumed to be the filename for the ENDF file. 275 276 Returns 277 ------- 278 openmc.data.AtomicRelaxation 279 Atomic relaxation data 280 281 """ 282 if isinstance(ev_or_filename, Evaluation): 283 ev = ev_or_filename 284 else: 285 ev = Evaluation(ev_or_filename) 286 287 # Atomic relaxation data is always MF=28, MT=533 288 if (28, 533) not in ev.section: 289 raise IOError('{} does not appear to be an atomic relaxation ' 290 'sublibrary.'.format(ev)) 291 292 # Determine number of subshells 293 file_obj = StringIO(ev.section[28, 533]) 294 params = get_head_record(file_obj) 295 n_subshells = params[4] 296 297 # Create data dictionaries 298 binding_energy = {} 299 num_electrons = {} 300 transitions = {} 301 columns = ['secondary', 'tertiary', 'energy (eV)', 'probability'] 302 303 # Read data for each subshell 304 for i in range(n_subshells): 305 params, list_items = get_list_record(file_obj) 306 subi = _SUBSHELLS[int(params[0])] 307 n_transitions = int(params[5]) 308 binding_energy[subi] = list_items[0] 309 num_electrons[subi] = list_items[1] 310 311 if n_transitions > 0: 312 # Read transition data 313 records = [] 314 for j in range(n_transitions): 315 subj = _SUBSHELLS[int(list_items[6*(j+1)])] 316 subk = _SUBSHELLS[int(list_items[6*(j+1) + 1])] 317 etr = list_items[6*(j+1) + 2] 318 ftr = list_items[6*(j+1) + 3] 319 records.append((subj, subk, etr, ftr)) 320 321 # Create dataframe for transitions 322 transitions[subi] = pd.DataFrame.from_records( 323 records, columns=columns) 324 325 # Return instance of class 326 return cls(binding_energy, num_electrons, transitions) 327 328 @classmethod 329 def from_hdf5(cls, group): 330 """Generate atomic relaxation data from an HDF5 group 331 332 Parameters 333 ---------- 334 group : h5py.Group 335 HDF5 group to read from 336 337 Returns 338 ------- 339 openmc.data.AtomicRelaxation 340 Atomic relaxation data 341 342 """ 343 # Create data dictionaries 344 binding_energy = {} 345 num_electrons = {} 346 transitions = {} 347 348 designators = [s.decode() for s in group.attrs['designators']] 349 columns = ['secondary', 'tertiary', 'energy (eV)', 'probability'] 350 for shell in designators: 351 # Shell group 352 sub_group = group[shell] 353 354 # Read subshell binding energy and number of electrons 355 if 'binding_energy' in sub_group.attrs: 356 binding_energy[shell] = sub_group.attrs['binding_energy'] 357 if 'num_electrons' in sub_group.attrs: 358 num_electrons[shell] = sub_group.attrs['num_electrons'] 359 360 # Read transition data 361 if 'transitions' in sub_group: 362 df = pd.DataFrame(sub_group['transitions'][()], 363 columns=columns) 364 # Replace float indexes back to subshell strings 365 df[columns[:2]] = df[columns[:2]].replace( 366 np.arange(float(len(_SUBSHELLS))), _SUBSHELLS) 367 transitions[shell] = df 368 369 return cls(binding_energy, num_electrons, transitions) 370 371 def to_hdf5(self, group, shell): 372 """Write atomic relaxation data to an HDF5 group 373 374 Parameters 375 ---------- 376 group : h5py.Group 377 HDF5 group to write to 378 shell : str 379 The subshell to write data for 380 381 """ 382 383 # Write subshell binding energy and number of electrons 384 group.attrs['binding_energy'] = self.binding_energy[shell] 385 group.attrs['num_electrons'] = self.num_electrons[shell] 386 387 # Write transition data with replacements 388 if shell in self.transitions: 389 df = self.transitions[shell].replace( 390 _SUBSHELLS, range(len(_SUBSHELLS))) 391 group.create_dataset('transitions', data=df.values.astype(float)) 392 393 394class IncidentPhoton(EqualityMixin): 395 r"""Photon interaction data. 396 397 This class stores photo-atomic, photo-nuclear, atomic relaxation, 398 Compton profile, stopping power, and bremsstrahlung data assembled from 399 different sources. To create an instance, the factory method 400 :meth:`IncidentPhoton.from_endf` can be used. To add atomic relaxation or 401 Compton profile data, set the :attr:`IncidentPhoton.atomic_relaxation` and 402 :attr:`IncidentPhoton.compton_profiles` attributes directly. 403 404 Parameters 405 ---------- 406 atomic_number : int 407 Number of protons in the target nucleus 408 409 Attributes 410 ---------- 411 atomic_number : int 412 Number of protons in the target nucleus 413 atomic_relaxation : openmc.data.AtomicRelaxation or None 414 Atomic relaxation data 415 bremsstrahlung : dict 416 Dictionary of bremsstrahlung data with keys 'I' (mean excitation energy 417 in [eV]), 'num_electrons' (number of electrons in each subshell), 418 'ionization_energy' (ionization potential of each subshell), 419 'electron_energy' (incident electron kinetic energy values in [eV]), 420 'photon_energy' (ratio of the energy of the emitted photon to the 421 incident electron kinetic energy), and 'dcs' (cross section values in 422 [b]). The cross sections are in scaled form: :math:`(\beta^2/Z^2) E_k 423 (d\sigma/dE_k)`, where :math:`E_k` is the energy of the emitted photon. 424 A negative number of electrons in a subshell indicates conduction 425 electrons. 426 compton_profiles : dict 427 Dictionary of Compton profile data with keys 'num_electrons' (number of 428 electrons in each subshell), 'binding_energy' (ionization potential of 429 each subshell), and 'J' (Hartree-Fock Compton profile as a function of 430 the projection of the electron momentum on the scattering vector, 431 :math:`p_z` for each subshell). Note that subshell occupancies may not 432 match the atomic relaxation data. 433 reactions : collections.OrderedDict 434 Contains the cross sections for each photon reaction. The keys are MT 435 values and the values are instances of :class:`PhotonReaction`. 436 437 """ 438 439 def __init__(self, atomic_number): 440 self.atomic_number = atomic_number 441 self._atomic_relaxation = None 442 self.reactions = OrderedDict() 443 self.compton_profiles = {} 444 self.bremsstrahlung = {} 445 446 def __contains__(self, mt): 447 return mt in self.reactions 448 449 def __getitem__(self, mt): 450 if mt in self.reactions: 451 return self.reactions[mt] 452 else: 453 raise KeyError('No reaction with MT={}.'.format(mt)) 454 455 def __repr__(self): 456 return "<IncidentPhoton: {}>".format(self.name) 457 458 def __iter__(self): 459 return iter(self.reactions.values()) 460 461 @property 462 def atomic_number(self): 463 return self._atomic_number 464 465 @property 466 def atomic_relaxation(self): 467 return self._atomic_relaxation 468 469 @property 470 def name(self): 471 return ATOMIC_SYMBOL[self.atomic_number] 472 473 @atomic_number.setter 474 def atomic_number(self, atomic_number): 475 cv.check_type('atomic number', atomic_number, Integral) 476 cv.check_greater_than('atomic number', atomic_number, 0, True) 477 self._atomic_number = atomic_number 478 479 @atomic_relaxation.setter 480 def atomic_relaxation(self, atomic_relaxation): 481 cv.check_type('atomic relaxation data', atomic_relaxation, 482 AtomicRelaxation) 483 self._atomic_relaxation = atomic_relaxation 484 485 @classmethod 486 def from_ace(cls, ace_or_filename): 487 """Generate incident photon data from an ACE table 488 489 Parameters 490 ---------- 491 ace_or_filename : str or openmc.data.ace.Table 492 ACE table to read from. If given as a string, it is assumed to be 493 the filename for the ACE file. 494 495 Returns 496 ------- 497 openmc.data.IncidentPhoton 498 Photon interaction data 499 500 """ 501 # First obtain the data for the first provided ACE table/file 502 if isinstance(ace_or_filename, Table): 503 ace = ace_or_filename 504 else: 505 ace = get_table(ace_or_filename) 506 507 # Get atomic number based on name of ACE table 508 zaid, xs = ace.name.split('.') 509 if not xs.endswith('p'): 510 raise TypeError("{} is not a photoatomic transport ACE table.".format(ace)) 511 Z = get_metadata(int(zaid))[2] 512 513 # Read each reaction 514 data = cls(Z) 515 for mt in (502, 504, 515, 522, 525): 516 data.reactions[mt] = PhotonReaction.from_ace(ace, mt) 517 518 # Get heating cross sections [eV-barn] from factors [eV per collision] 519 # by multiplying with total xs 520 data.reactions[525].xs.y *= sum([data.reactions[mt].xs.y for mt in 521 (502, 504, 515, 522)]) 522 523 # Compton profiles 524 n_shell = ace.nxs[5] 525 if n_shell != 0: 526 # Get number of electrons in each shell 527 idx = ace.jxs[6] 528 data.compton_profiles['num_electrons'] = ace.xss[idx : idx+n_shell] 529 530 # Get binding energy for each shell 531 idx = ace.jxs[7] 532 e = ace.xss[idx : idx+n_shell]*EV_PER_MEV 533 data.compton_profiles['binding_energy'] = e 534 535 # Create Compton profile for each electron shell 536 profiles = [] 537 for k in range(n_shell): 538 # Get number of momentum values and interpolation scheme 539 loca = int(ace.xss[ace.jxs[9] + k]) 540 jj = int(ace.xss[ace.jxs[10] + loca - 1]) 541 m = int(ace.xss[ace.jxs[10] + loca]) 542 543 # Read momentum and PDF 544 idx = ace.jxs[10] + loca + 1 545 pz = ace.xss[idx : idx+m] 546 pdf = ace.xss[idx+m : idx+2*m] 547 548 # Create proflie function 549 J_k = Tabulated1D(pz, pdf, [m], [jj]) 550 profiles.append(J_k) 551 data.compton_profiles['J'] = profiles 552 553 # Subshell photoelectric xs and atomic relaxation data 554 if ace.nxs[7] > 0: 555 data.atomic_relaxation = AtomicRelaxation.from_ace(ace) 556 557 # Get subshell designators 558 n_subshells = ace.nxs[7] 559 idx = ace.jxs[11] 560 designators = [int(i) for i in ace.xss[idx : idx+n_subshells]] 561 562 # Get energy grid for subshell photoionization 563 n_energy = ace.nxs[3] 564 idx = ace.jxs[1] 565 energy = np.exp(ace.xss[idx : idx+n_energy])*EV_PER_MEV 566 567 # Get cross section for each subshell 568 idx = ace.jxs[16] 569 for d in designators: 570 # Create photon reaction 571 mt = 533 + d 572 rx = PhotonReaction(mt) 573 data.reactions[mt] = rx 574 575 # Store cross section, determining threshold 576 xs = ace.xss[idx : idx+n_energy].copy() 577 nonzero = (xs != 0.0) 578 xs[nonzero] = np.exp(xs[nonzero]) 579 threshold = np.where(xs > 0.0)[0][0] 580 rx.xs = Tabulated1D(energy[threshold:], xs[threshold:], 581 [n_energy - threshold], [5]) 582 idx += n_energy 583 584 # Copy binding energy 585 shell = _SUBSHELLS[d] 586 e = data.atomic_relaxation.binding_energy[shell] 587 rx.subshell_binding_energy = e 588 else: 589 raise ValueError("ACE table {} does not have subshell data. Only " 590 "newer ACE photoatomic libraries are supported " 591 "(e.g., eprdata14).".format(ace.name)) 592 593 # Add bremsstrahlung DCS data 594 data._add_bremsstrahlung() 595 596 return data 597 598 @classmethod 599 def from_endf(cls, photoatomic, relaxation=None): 600 """Generate incident photon data from an ENDF evaluation 601 602 Parameters 603 ---------- 604 photoatomic : str or openmc.data.endf.Evaluation 605 ENDF photoatomic data evaluation to read from. If given as a string, 606 it is assumed to be the filename for the ENDF file. 607 relaxation : str or openmc.data.endf.Evaluation, optional 608 ENDF atomic relaxation data evaluation to read from. If given as a 609 string, it is assumed to be the filename for the ENDF file. 610 611 Returns 612 ------- 613 openmc.data.IncidentPhoton 614 Photon interaction data 615 616 """ 617 if isinstance(photoatomic, Evaluation): 618 ev = photoatomic 619 else: 620 ev = Evaluation(photoatomic) 621 622 Z = ev.target['atomic_number'] 623 data = cls(Z) 624 625 # Read each reaction 626 for mf, mt, nc, mod in ev.reaction_list: 627 if mf == 23: 628 data.reactions[mt] = PhotonReaction.from_endf(ev, mt) 629 630 # Add atomic relaxation data if it hasn't been added already 631 if relaxation is not None: 632 data.atomic_relaxation = AtomicRelaxation.from_endf(relaxation) 633 634 # If Compton profile data hasn't been loaded, do so 635 if not _COMPTON_PROFILES: 636 filename = os.path.join(os.path.dirname(__file__), 'compton_profiles.h5') 637 with h5py.File(filename, 'r') as f: 638 _COMPTON_PROFILES['pz'] = f['pz'][()] 639 for i in range(1, 101): 640 group = f['{:03}'.format(i)] 641 num_electrons = group['num_electrons'][()] 642 binding_energy = group['binding_energy'][()]*EV_PER_MEV 643 J = group['J'][()] 644 _COMPTON_PROFILES[i] = {'num_electrons': num_electrons, 645 'binding_energy': binding_energy, 646 'J': J} 647 648 # Add Compton profile data 649 pz = _COMPTON_PROFILES['pz'] 650 profile = _COMPTON_PROFILES[Z] 651 data.compton_profiles['num_electrons'] = profile['num_electrons'] 652 data.compton_profiles['binding_energy'] = profile['binding_energy'] 653 data.compton_profiles['J'] = [Tabulated1D(pz, J_k) for J_k in profile['J']] 654 655 # Add bremsstrahlung DCS data 656 data._add_bremsstrahlung() 657 658 return data 659 660 @classmethod 661 def from_hdf5(cls, group_or_filename): 662 """Generate photon reaction from an HDF5 group 663 664 Parameters 665 ---------- 666 group_or_filename : h5py.Group or str 667 HDF5 group containing interaction data. If given as a string, it is 668 assumed to be the filename for the HDF5 file, and the first group is 669 used to read from. 670 671 Returns 672 ------- 673 openmc.data.IncidentPhoton 674 Photon interaction data 675 676 """ 677 if isinstance(group_or_filename, h5py.Group): 678 group = group_or_filename 679 need_to_close = False 680 else: 681 h5file = h5py.File(str(group_or_filename), 'r') 682 need_to_close = True 683 684 # Make sure version matches 685 if 'version' in h5file.attrs: 686 major, minor = h5file.attrs['version'] 687 # For now all versions of HDF5 data can be read 688 else: 689 raise IOError( 690 'HDF5 data does not indicate a version. Your installation ' 691 'of the OpenMC Python API expects version {}.x data.' 692 .format(HDF5_VERSION_MAJOR)) 693 694 group = list(h5file.values())[0] 695 696 Z = group.attrs['Z'] 697 data = cls(Z) 698 699 # Read energy grid 700 energy = group['energy'][()] 701 702 # Read cross section data 703 for mt, (name, key) in _REACTION_NAME.items(): 704 if key in group: 705 rgroup = group[key] 706 elif key in group['subshells']: 707 rgroup = group['subshells'][key] 708 else: 709 continue 710 711 data.reactions[mt] = PhotonReaction.from_hdf5(rgroup, mt, energy) 712 713 # Check for necessary reactions 714 for mt in (502, 504, 522): 715 assert mt in data, "Reaction {} not found".format(mt) 716 717 # Read atomic relaxation 718 data.atomic_relaxation = AtomicRelaxation.from_hdf5(group['subshells']) 719 720 # Read Compton profiles 721 if 'compton_profiles' in group: 722 rgroup = group['compton_profiles'] 723 profile = data.compton_profiles 724 profile['num_electrons'] = rgroup['num_electrons'][()] 725 profile['binding_energy'] = rgroup['binding_energy'][()] 726 727 # Get electron momentum values 728 pz = rgroup['pz'][()] 729 J = rgroup['J'][()] 730 if pz.size != J.shape[1]: 731 raise ValueError("'J' array shape is not consistent with the " 732 "'pz' array shape") 733 profile['J'] = [Tabulated1D(pz, Jk) for Jk in J] 734 735 # Read bremsstrahlung 736 if 'bremsstrahlung' in group: 737 rgroup = group['bremsstrahlung'] 738 data.bremsstrahlung['I'] = rgroup.attrs['I'] 739 for key in ('dcs', 'electron_energy', 'ionization_energy', 740 'num_electrons', 'photon_energy'): 741 data.bremsstrahlung[key] = rgroup[key][()] 742 743 # If HDF5 file was opened here, make sure it gets closed 744 if need_to_close: 745 h5file.close() 746 747 return data 748 749 def export_to_hdf5(self, path, mode='a', libver='earliest'): 750 """Export incident photon data to an HDF5 file. 751 752 Parameters 753 ---------- 754 path : str 755 Path to write HDF5 file to 756 mode : {'r+', 'w', 'x', 'a'} 757 Mode that is used to open the HDF5 file. This is the second argument 758 to the :class:`h5py.File` constructor. 759 libver : {'earliest', 'latest'} 760 Compatibility mode for the HDF5 file. 'latest' will produce files 761 that are less backwards compatible but have performance benefits. 762 763 """ 764 with h5py.File(str(path), mode, libver=libver) as f: 765 # Write filetype and version 766 f.attrs['filetype'] = np.string_('data_photon') 767 if 'version' not in f.attrs: 768 f.attrs['version'] = np.array(HDF5_VERSION) 769 770 group = f.create_group(self.name) 771 group.attrs['Z'] = Z = self.atomic_number 772 773 # Determine union energy grid 774 union_grid = np.array([]) 775 for rx in self: 776 union_grid = np.union1d(union_grid, rx.xs.x) 777 group.create_dataset('energy', data=union_grid) 778 779 # Write cross sections 780 shell_group = group.create_group('subshells') 781 designators = [] 782 for mt, rx in self.reactions.items(): 783 name, key = _REACTION_NAME[mt] 784 if mt in (502, 504, 515, 517, 522, 525): 785 sub_group = group.create_group(key) 786 elif mt >= 534 and mt <= 572: 787 # Subshell 788 designators.append(key) 789 sub_group = shell_group.create_group(key) 790 791 # Write atomic relaxation 792 if self.atomic_relaxation is not None: 793 if key in self.atomic_relaxation.subshells: 794 self.atomic_relaxation.to_hdf5(sub_group, key) 795 else: 796 continue 797 798 rx.to_hdf5(sub_group, union_grid, Z) 799 800 shell_group.attrs['designators'] = np.array(designators, dtype='S') 801 802 # Write Compton profiles 803 if self.compton_profiles: 804 compton_group = group.create_group('compton_profiles') 805 806 profile = self.compton_profiles 807 compton_group.create_dataset('num_electrons', 808 data=profile['num_electrons']) 809 compton_group.create_dataset('binding_energy', 810 data=profile['binding_energy']) 811 812 # Get electron momentum values 813 compton_group.create_dataset('pz', data=profile['J'][0].x) 814 815 # Create/write 2D array of profiles 816 J = np.array([Jk.y for Jk in profile['J']]) 817 compton_group.create_dataset('J', data=J) 818 819 # Write bremsstrahlung 820 if self.bremsstrahlung: 821 brem_group = group.create_group('bremsstrahlung') 822 for key, value in self.bremsstrahlung.items(): 823 if key == 'I': 824 brem_group.attrs[key] = value 825 else: 826 brem_group.create_dataset(key, data=value) 827 828 def _add_bremsstrahlung(self): 829 """Add the data used in the thick-target bremsstrahlung approximation 830 831 """ 832 # Load bremsstrahlung data if it has not yet been loaded 833 if not _BREMSSTRAHLUNG: 834 # Add data used for density effect correction 835 filename = os.path.join(os.path.dirname(__file__), 'density_effect.h5') 836 with h5py.File(filename, 'r') as f: 837 for i in range(1, 101): 838 group = f['{:03}'.format(i)] 839 _BREMSSTRAHLUNG[i] = { 840 'I': group.attrs['I'], 841 'num_electrons': group['num_electrons'][()], 842 'ionization_energy': group['ionization_energy'][()] 843 } 844 845 filename = os.path.join(os.path.dirname(__file__), 'BREMX.DAT') 846 with open(filename, 'r') as fh: 847 brem = fh.read().split() 848 849 # Incident electron kinetic energy grid in eV 850 _BREMSSTRAHLUNG['electron_energy'] = np.logspace(3, 9, 200) 851 log_energy = np.log(_BREMSSTRAHLUNG['electron_energy']) 852 853 # Get number of tabulated electron and photon energy values 854 n = int(brem[37]) 855 k = int(brem[38]) 856 857 # Index in data 858 p = 39 859 860 # Get log of incident electron kinetic energy values, used for 861 # cubic spline interpolation in log energy. Units are in MeV, so 862 # convert to eV. 863 logx = np.log(np.fromiter(brem[p:p+n], float, n)*EV_PER_MEV) 864 p += n 865 866 # Get reduced photon energy values 867 _BREMSSTRAHLUNG['photon_energy'] = np.fromiter(brem[p:p+k], float, k) 868 p += k 869 870 for i in range(1, 101): 871 dcs = np.empty([len(log_energy), k]) 872 873 # Get the scaled cross section values for each electron energy 874 # and reduced photon energy for this Z. Units are in mb, so 875 # convert to b. 876 y = np.reshape(np.fromiter(brem[p:p+n*k], float, n*k), (n, k))*1.0e-3 877 p += k*n 878 879 for j in range(k): 880 # Cubic spline interpolation in log energy and linear DCS 881 cs = CubicSpline(logx, y[:, j]) 882 883 # Get scaled DCS values (barns) on new energy grid 884 dcs[:, j] = cs(log_energy) 885 886 _BREMSSTRAHLUNG[i]['dcs'] = dcs 887 888 # Add bremsstrahlung DCS data 889 self.bremsstrahlung['electron_energy'] = _BREMSSTRAHLUNG['electron_energy'] 890 self.bremsstrahlung['photon_energy'] = _BREMSSTRAHLUNG['photon_energy'] 891 self.bremsstrahlung.update(_BREMSSTRAHLUNG[self.atomic_number]) 892 893 894class PhotonReaction(EqualityMixin): 895 """Photon-induced reaction 896 897 Parameters 898 ---------- 899 mt : int 900 The ENDF MT number for this reaction. 901 902 Attributes 903 ---------- 904 anomalous_real : openmc.data.Tabulated1D 905 Real part of the anomalous scattering factor 906 anomlaous_imag : openmc.data.Tabulated1D 907 Imaginary part of the anomalous scatttering factor 908 mt : int 909 The ENDF MT number for this reaction. 910 scattering_factor : openmc.data.Tabulated1D 911 Coherent or incoherent form factor. 912 xs : Callable 913 Cross section as a function of incident photon energy 914 915 """ 916 917 def __init__(self, mt): 918 self.mt = mt 919 self._xs = None 920 self._scattering_factor = None 921 self._anomalous_real = None 922 self._anomalous_imag = None 923 924 def __repr__(self): 925 if self.mt in _REACTION_NAME: 926 return "<Photon Reaction: MT={} {}>".format( 927 self.mt, _REACTION_NAME[self.mt][0]) 928 else: 929 return "<Photon Reaction: MT={}>".format(self.mt) 930 931 @property 932 def anomalous_real(self): 933 return self._anomalous_real 934 935 @property 936 def anomalous_imag(self): 937 return self._anomalous_imag 938 939 @property 940 def scattering_factor(self): 941 return self._scattering_factor 942 943 @property 944 def xs(self): 945 return self._xs 946 947 @anomalous_real.setter 948 def anomalous_real(self, anomalous_real): 949 cv.check_type('real part of anomalous scattering factor', 950 anomalous_real, Callable) 951 self._anomalous_real = anomalous_real 952 953 @anomalous_imag.setter 954 def anomalous_imag(self, anomalous_imag): 955 cv.check_type('imaginary part of anomalous scattering factor', 956 anomalous_imag, Callable) 957 self._anomalous_imag = anomalous_imag 958 959 @scattering_factor.setter 960 def scattering_factor(self, scattering_factor): 961 cv.check_type('scattering factor', scattering_factor, Callable) 962 self._scattering_factor = scattering_factor 963 964 @xs.setter 965 def xs(self, xs): 966 cv.check_type('reaction cross section', xs, Callable) 967 self._xs = xs 968 969 @classmethod 970 def from_ace(cls, ace, mt): 971 """Generate photon reaction from an ACE table 972 973 Parameters 974 ---------- 975 ace : openmc.data.ace.Table 976 ACE table to read from 977 mt : int 978 The MT value of the reaction to get data for 979 980 Returns 981 ------- 982 openmc.data.PhotonReaction 983 Photon reaction data 984 985 """ 986 # Create instance 987 rx = cls(mt) 988 989 # Get energy grid (stored as logarithms) 990 n = ace.nxs[3] 991 idx = ace.jxs[1] 992 energy = np.exp(ace.xss[idx : idx+n])*EV_PER_MEV 993 994 # Get index for appropriate reaction 995 if mt == 502: 996 # Coherent scattering 997 idx = ace.jxs[1] + 2*n 998 elif mt == 504: 999 # Incoherent scattering 1000 idx = ace.jxs[1] + n 1001 elif mt == 515: 1002 # Pair production 1003 idx = ace.jxs[1] + 4*n 1004 elif mt == 522: 1005 # Photoelectric 1006 idx = ace.jxs[1] + 3*n 1007 elif mt == 525: 1008 # Heating 1009 idx = ace.jxs[5] 1010 else: 1011 raise ValueError('ACE photoatomic cross sections do not have ' 1012 'data for MT={}.'.format(mt)) 1013 1014 # Store cross section 1015 xs = ace.xss[idx : idx+n].copy() 1016 if mt == 525: 1017 # Get heating factors in [eV per collision] 1018 xs *= EV_PER_MEV 1019 else: 1020 nonzero = (xs != 0.0) 1021 xs[nonzero] = np.exp(xs[nonzero]) 1022 rx.xs = Tabulated1D(energy, xs, [n], [5]) 1023 1024 # Get form factors for incoherent/coherent scattering 1025 new_format = (ace.nxs[6] > 0) 1026 if mt == 502: 1027 idx = ace.jxs[3] 1028 if new_format: 1029 n = (ace.jxs[4] - ace.jxs[3]) // 3 1030 x = ace.xss[idx : idx+n] 1031 idx += n 1032 else: 1033 x = np.array([ 1034 0.0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.08, 0.1, 0.12, 1035 0.15, 0.18, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 1036 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1037 1.7, 1.8, 1.9, 2.0, 2.2, 2.4, 2.6, 2.8, 3.0, 3.2, 3.4, 1038 3.6, 3.8, 4.0, 4.2, 4.4, 4.6, 4.8, 5.0, 5.2, 5.4, 5.6, 1039 5.8, 6.0]) 1040 n = x.size 1041 ff = ace.xss[idx+n : idx+2*n] 1042 rx.scattering_factor = Tabulated1D(x, ff) 1043 1044 elif mt == 504: 1045 idx = ace.jxs[2] 1046 if new_format: 1047 n = (ace.jxs[3] - ace.jxs[2]) // 2 1048 x = ace.xss[idx : idx+n] 1049 idx += n 1050 else: 1051 x = np.array([ 1052 0.0, 0.005, 0.01, 0.05, 0.1, 0.15, 0.2, 0.3, 0.4, 0.5, 0.6, 1053 0.7, 0.8, 0.9, 1.0, 1.5, 2.0, 3.0, 4.0, 5.0, 8.0 1054 ]) 1055 n = x.size 1056 ff = ace.xss[idx : idx+n] 1057 rx.scattering_factor = Tabulated1D(x, ff) 1058 1059 return rx 1060 1061 @classmethod 1062 def from_endf(cls, ev, mt): 1063 """Generate photon reaction from an ENDF evaluation 1064 1065 Parameters 1066 ---------- 1067 ev : openmc.data.endf.Evaluation 1068 ENDF photo-atomic interaction data evaluation 1069 mt : int 1070 The MT value of the reaction to get data for 1071 1072 Returns 1073 ------- 1074 openmc.data.PhotonReaction 1075 Photon reaction data 1076 1077 """ 1078 rx = cls(mt) 1079 1080 # Read photon cross section 1081 if (23, mt) in ev.section: 1082 file_obj = StringIO(ev.section[23, mt]) 1083 get_head_record(file_obj) 1084 params, rx.xs = get_tab1_record(file_obj) 1085 1086 # Set subshell binding energy and/or fluorescence yield 1087 if mt >= 534 and mt <= 599: 1088 rx.subshell_binding_energy = params[0] 1089 if mt >= 534 and mt <= 572: 1090 rx.fluorescence_yield = params[1] 1091 1092 # Read form factors / scattering functions 1093 if (27, mt) in ev.section: 1094 file_obj = StringIO(ev.section[27, mt]) 1095 get_head_record(file_obj) 1096 params, rx.scattering_factor = get_tab1_record(file_obj) 1097 1098 # Check for anomalous scattering factor 1099 if mt == 502: 1100 if (27, 506) in ev.section: 1101 file_obj = StringIO(ev.section[27, 506]) 1102 get_head_record(file_obj) 1103 params, rx.anomalous_real = get_tab1_record(file_obj) 1104 1105 if (27, 505) in ev.section: 1106 file_obj = StringIO(ev.section[27, 505]) 1107 get_head_record(file_obj) 1108 params, rx.anomalous_imag = get_tab1_record(file_obj) 1109 1110 return rx 1111 1112 @classmethod 1113 def from_hdf5(cls, group, mt, energy): 1114 """Generate photon reaction from an HDF5 group 1115 1116 Parameters 1117 ---------- 1118 group : h5py.Group 1119 HDF5 group to read from 1120 mt : int 1121 The MT value of the reaction to get data for 1122 energy : Iterable of float 1123 arrays of energies at which cross sections are tabulated at 1124 1125 Returns 1126 ------- 1127 openmc.data.PhotonReaction 1128 Photon reaction data 1129 1130 """ 1131 # Create instance 1132 rx = cls(mt) 1133 1134 # Cross sections 1135 xs = group['xs'][()] 1136 # Replace zero elements to small non-zero to enable log-log 1137 xs[xs == 0.0] = np.exp(-500.0) 1138 1139 # Threshold 1140 threshold_idx = 0 1141 if 'threshold_idx' in group['xs'].attrs: 1142 threshold_idx = group['xs'].attrs['threshold_idx'] 1143 1144 # Store cross section 1145 rx.xs = Tabulated1D(energy[threshold_idx:], xs, [len(xs)], [5]) 1146 1147 # Check for anomalous scattering factor 1148 if 'anomalous_real' in group: 1149 rx.anomalous_real = Tabulated1D.from_hdf5(group['anomalous_real']) 1150 if 'anomalous_imag' in group: 1151 rx.anomalous_imag = Tabulated1D.from_hdf5(group['anomalous_imag']) 1152 1153 # Check for factors / scattering functions 1154 if 'scattering_factor' in group: 1155 rx.scattering_factor = Tabulated1D.from_hdf5(group['scattering_factor']) 1156 1157 return rx 1158 1159 def to_hdf5(self, group, energy, Z): 1160 """Write photon reaction to an HDF5 group 1161 1162 Parameters 1163 ---------- 1164 group : h5py.Group 1165 HDF5 group to write to 1166 energy : Iterable of float 1167 arrays of energies at which cross sections are tabulated at 1168 Z : int 1169 atomic number 1170 1171 """ 1172 1173 # Write cross sections 1174 if self.mt >= 534 and self.mt <= 572: 1175 # Determine threshold 1176 threshold = self.xs.x[0] 1177 idx = np.searchsorted(energy, threshold, side='right') - 1 1178 1179 # Interpolate cross section onto union grid and write 1180 photoionization = self.xs(energy[idx:]) 1181 group.create_dataset('xs', data=photoionization) 1182 assert len(energy) == len(photoionization) + idx 1183 group['xs'].attrs['threshold_idx'] = idx 1184 else: 1185 group.create_dataset('xs', data=self.xs(energy)) 1186 1187 # Write scattering factor 1188 if self.scattering_factor is not None: 1189 if self.mt == 502: 1190 # Create integrated form factor 1191 ff = deepcopy(self.scattering_factor) 1192 ff.x *= ff.x 1193 ff.y *= ff.y/Z**2 1194 int_ff = Tabulated1D(ff.x, ff.integral()) 1195 int_ff.to_hdf5(group, 'integrated_scattering_factor') 1196 self.scattering_factor.to_hdf5(group, 'scattering_factor') 1197 if self.anomalous_real is not None: 1198 self.anomalous_real.to_hdf5(group, 'anomalous_real') 1199 if self.anomalous_imag is not None: 1200 self.anomalous_imag.to_hdf5(group, 'anomalous_imag') 1201