1from itertools import chain 2from numbers import Integral, Real 3import string 4 5import numpy as np 6 7import openmc.checkvalue as cv 8import openmc.data 9 10# Supported keywords for continuous-energy cross section plotting 11PLOT_TYPES = ['total', 'scatter', 'elastic', 'inelastic', 'fission', 12 'absorption', 'capture', 'nu-fission', 'nu-scatter', 'unity', 13 'slowing-down power', 'damage'] 14 15# Supported keywords for multi-group cross section plotting 16PLOT_TYPES_MGXS = ['total', 'absorption', 'scatter', 'fission', 17 'kappa-fission', 'nu-fission', 'prompt-nu-fission', 18 'deleyed-nu-fission', 'chi', 'chi-prompt', 'chi-delayed', 19 'inverse-velocity', 'beta', 'decay-rate', 'unity'] 20# Create a dictionary which can be used to convert PLOT_TYPES_MGXS to the 21# openmc.XSdata attribute name needed to access the data 22_PLOT_MGXS_ATTR = {line: line.replace(' ', '_').replace('-', '_') 23 for line in PLOT_TYPES_MGXS} 24_PLOT_MGXS_ATTR['scatter'] = 'scatter_matrix' 25 26# Special MT values 27UNITY_MT = -1 28XI_MT = -2 29 30# MTs to combine to generate associated plot_types 31_INELASTIC = [mt for mt in openmc.data.SUM_RULES[3] if mt != 27] 32PLOT_TYPES_MT = { 33 'total': openmc.data.SUM_RULES[1], 34 'scatter': [2] + _INELASTIC, 35 'elastic': [2], 36 'inelastic': _INELASTIC, 37 'fission': [18], 38 'absorption': [27], 39 'capture': [101], 40 'nu-fission': [18], 41 'nu-scatter': [2] + _INELASTIC, 42 'unity': [UNITY_MT], 43 'slowing-down power': [2] + [XI_MT], 44 'damage': [444] 45} 46 47# Types of plots to plot linearly in y 48PLOT_TYPES_LINEAR = {'nu-fission / fission', 'nu-scatter / scatter', 49 'nu-fission / absorption', 'fission / absorption'} 50 51# Minimum and maximum energies for plotting (units of eV) 52_MIN_E = 1.e-5 53_MAX_E = 20.e6 54 55 56def plot_xs(this, types, divisor_types=None, temperature=294., data_type=None, 57 axis=None, sab_name=None, ce_cross_sections=None, 58 mg_cross_sections=None, enrichment=None, plot_CE=True, orders=None, 59 divisor_orders=None, **kwargs): 60 """Creates a figure of continuous-energy cross sections for this item. 61 62 Parameters 63 ---------- 64 this : str or openmc.Material 65 Object to source data from 66 types : Iterable of values of PLOT_TYPES 67 The type of cross sections to include in the plot. 68 divisor_types : Iterable of values of PLOT_TYPES, optional 69 Cross section types which will divide those produced by types 70 before plotting. A type of 'unity' can be used to effectively not 71 divide some types. 72 temperature : float, optional 73 Temperature in Kelvin to plot. If not specified, a default 74 temperature of 294K will be plotted. Note that the nearest 75 temperature in the library for each nuclide will be used as opposed 76 to using any interpolation. 77 data_type : {'nuclide', 'element', 'material', 'macroscopic'}, optional 78 Type of object to plot. If not specified, a guess is made based on the 79 `this` argument. 80 axis : matplotlib.axes, optional 81 A previously generated axis to use for plotting. If not specified, 82 a new axis and figure will be generated. 83 sab_name : str, optional 84 Name of S(a,b) library to apply to MT=2 data when applicable; only used 85 for items which are instances of openmc.Element or openmc.Nuclide 86 ce_cross_sections : str, optional 87 Location of cross_sections.xml file. Default is None. 88 mg_cross_sections : str, optional 89 Location of MGXS HDF5 Library file. Default is None. 90 enrichment : float, optional 91 Enrichment for U235 in weight percent. For example, input 4.95 for 92 4.95 weight percent enriched U. Default is None. This is only used for 93 items which are instances of openmc.Element 94 plot_CE : bool, optional 95 Denotes whether or not continuous-energy will be plotted. Defaults to 96 plotting the continuous-energy data. 97 orders : Iterable of Integral, optional 98 The scattering order or delayed group index to use for the 99 corresponding entry in types. Defaults to the 0th order for scattering 100 and the total delayed neutron data. This only applies to plots of 101 multi-group data. 102 divisor_orders : Iterable of Integral, optional 103 Same as orders, but for divisor_types 104 **kwargs 105 All keyword arguments are passed to 106 :func:`matplotlib.pyplot.figure`. 107 108 Returns 109 ------- 110 fig : matplotlib.figure.Figure 111 If axis is None, then a Matplotlib Figure of the generated 112 cross section will be returned. Otherwise, a value of 113 None will be returned as the figure and axes have already been 114 generated. 115 116 """ 117 import matplotlib.pyplot as plt 118 119 cv.check_type("plot_CE", plot_CE, bool) 120 121 if data_type is None: 122 if isinstance(this, openmc.Nuclide): 123 data_type = 'nuclide' 124 elif isinstance(this, openmc.Element): 125 data_type = 'element' 126 elif isinstance(this, openmc.Material): 127 data_type = 'material' 128 elif isinstance(this, openmc.Macroscopic): 129 data_type = 'macroscopic' 130 elif isinstance(this, str): 131 if this[-1] in string.digits: 132 data_type = 'nuclide' 133 else: 134 data_type = 'element' 135 else: 136 raise TypeError("Invalid type for plotting") 137 138 if plot_CE: 139 # Calculate for the CE cross sections 140 E, data = calculate_cexs(this, data_type, types, temperature, sab_name, 141 ce_cross_sections, enrichment) 142 if divisor_types: 143 cv.check_length('divisor types', divisor_types, len(types)) 144 Ediv, data_div = calculate_cexs(this, divisor_types, temperature, 145 sab_name, ce_cross_sections, 146 enrichment) 147 148 # Create a new union grid, interpolate data and data_div on to that 149 # grid, and then do the actual division 150 Enum = E[:] 151 E = np.union1d(Enum, Ediv) 152 data_new = np.zeros((len(types), len(E))) 153 154 for line in range(len(types)): 155 data_new[line, :] = \ 156 np.divide(np.interp(E, Enum, data[line, :]), 157 np.interp(E, Ediv, data_div[line, :])) 158 if divisor_types[line] != 'unity': 159 types[line] = types[line] + ' / ' + divisor_types[line] 160 data = data_new 161 else: 162 # Calculate for MG cross sections 163 E, data = calculate_mgxs(this, data_type, types, orders, temperature, 164 mg_cross_sections, ce_cross_sections, 165 enrichment) 166 167 if divisor_types: 168 cv.check_length('divisor types', divisor_types, len(types)) 169 Ediv, data_div = calculate_mgxs(this, data_type, divisor_types, 170 divisor_orders, temperature, 171 mg_cross_sections, 172 ce_cross_sections, enrichment) 173 174 # Perform the division 175 for line in range(len(types)): 176 data[line, :] /= data_div[line, :] 177 if divisor_types[line] != 'unity': 178 types[line] += ' / ' + divisor_types[line] 179 180 # Generate the plot 181 if axis is None: 182 fig, ax = plt.subplots() 183 else: 184 fig = None 185 ax = axis 186 # Set to loglog or semilogx depending on if we are plotting a data 187 # type which we expect to vary linearly 188 if set(types).issubset(PLOT_TYPES_LINEAR): 189 plot_func = ax.semilogx 190 else: 191 plot_func = ax.loglog 192 193 # Plot the data 194 for i in range(len(data)): 195 data[i, :] = np.nan_to_num(data[i, :]) 196 if np.sum(data[i, :]) > 0.: 197 plot_func(E, data[i, :], label=types[i]) 198 199 ax.set_xlabel('Energy [eV]') 200 if plot_CE: 201 ax.set_xlim(_MIN_E, _MAX_E) 202 else: 203 ax.set_xlim(E[-1], E[0]) 204 if divisor_types: 205 if data_type == 'nuclide': 206 ylabel = 'Nuclidic Microscopic Data' 207 elif data_type == 'element': 208 ylabel = 'Elemental Microscopic Data' 209 elif data_type == 'material' or data_type == 'macroscopic': 210 ylabel = 'Macroscopic Data' 211 else: 212 if data_type == 'nuclide': 213 ylabel = 'Microscopic Cross Section [b]' 214 elif data_type == 'element': 215 ylabel = 'Elemental Cross Section [b]' 216 elif data_type == 'material' or data_type == 'macroscopic': 217 ylabel = 'Macroscopic Cross Section [1/cm]' 218 ax.set_ylabel(ylabel) 219 ax.legend(loc='best') 220 name = this.name if data_type == 'material' else this 221 if len(types) > 1: 222 ax.set_title('Cross Sections for ' + name) 223 else: 224 ax.set_title('Cross Section for ' + name) 225 226 return fig 227 228 229def calculate_cexs(this, data_type, types, temperature=294., sab_name=None, 230 cross_sections=None, enrichment=None): 231 """Calculates continuous-energy cross sections of a requested type. 232 233 Parameters 234 ---------- 235 this : {str, openmc.Nuclide, openmc.Element, openmc.Material} 236 Object to source data from 237 data_type : {'nuclide', 'element', 'material'} 238 Type of object to plot 239 types : Iterable of values of PLOT_TYPES 240 The type of cross sections to calculate 241 temperature : float, optional 242 Temperature in Kelvin to plot. If not specified, a default 243 temperature of 294K will be plotted. Note that the nearest 244 temperature in the library for each nuclide will be used as opposed 245 to using any interpolation. 246 sab_name : str, optional 247 Name of S(a,b) library to apply to MT=2 data when applicable. 248 cross_sections : str, optional 249 Location of cross_sections.xml file. Default is None. 250 enrichment : float, optional 251 Enrichment for U235 in weight percent. For example, input 4.95 for 252 4.95 weight percent enriched U. Default is None 253 (natural composition). 254 255 Returns 256 ------- 257 energy_grid : numpy.ndarray 258 Energies at which cross sections are calculated, in units of eV 259 data : numpy.ndarray 260 Cross sections calculated at the energy grid described by energy_grid 261 262 """ 263 264 # Check types 265 cv.check_type('temperature', temperature, Real) 266 if sab_name: 267 cv.check_type('sab_name', sab_name, str) 268 if enrichment: 269 cv.check_type('enrichment', enrichment, Real) 270 271 if data_type == 'nuclide': 272 if isinstance(this, str): 273 nuc = openmc.Nuclide(this) 274 else: 275 nuc = this 276 energy_grid, xs = _calculate_cexs_nuclide(nuc, types, temperature, 277 sab_name, cross_sections) 278 # Convert xs (Iterable of Callable) to a grid of cross section values 279 # calculated on the points in energy_grid for consistency with the 280 # element and material functions. 281 data = np.zeros((len(types), len(energy_grid))) 282 for line in range(len(types)): 283 data[line, :] = xs[line](energy_grid) 284 elif data_type == 'element': 285 if isinstance(this, str): 286 elem = openmc.Element(this) 287 else: 288 elem = this 289 energy_grid, data = _calculate_cexs_elem_mat(elem, types, temperature, 290 cross_sections, sab_name, 291 enrichment) 292 elif data_type == 'material': 293 cv.check_type('this', this, openmc.Material) 294 energy_grid, data = _calculate_cexs_elem_mat(this, types, temperature, 295 cross_sections) 296 else: 297 raise TypeError("Invalid type") 298 299 return energy_grid, data 300 301 302def _calculate_cexs_nuclide(this, types, temperature=294., sab_name=None, 303 cross_sections=None): 304 """Calculates continuous-energy cross sections of a requested type. 305 306 Parameters 307 ---------- 308 this : openmc.Nuclide 309 Nuclide object to source data from 310 types : Iterable of str or Integral 311 The type of cross sections to calculate; values can either be those 312 in openmc.PLOT_TYPES or keys from openmc.data.REACTION_MT which 313 correspond to a reaction description e.g '(n,2n)' or integers which 314 correspond to reaction channel (MT) numbers. 315 temperature : float, optional 316 Temperature in Kelvin to plot. If not specified, a default 317 temperature of 294K will be plotted. Note that the nearest 318 temperature in the library for each nuclide will be used as opposed 319 to using any interpolation. 320 sab_name : str, optional 321 Name of S(a,b) library to apply to MT=2 data when applicable. 322 cross_sections : str, optional 323 Location of cross_sections.xml file. Default is None. 324 325 Returns 326 ------- 327 energy_grid : numpy.ndarray 328 Energies at which cross sections are calculated, in units of eV 329 data : Iterable of Callable 330 Requested cross section functions 331 332 """ 333 334 # Load the library 335 library = openmc.data.DataLibrary.from_xml(cross_sections) 336 337 # Convert temperature to format needed for access in the library 338 strT = "{}K".format(int(round(temperature))) 339 T = temperature 340 341 # Now we can create the data sets to be plotted 342 energy_grid = [] 343 xs = [] 344 lib = library.get_by_material(this) 345 if lib is not None: 346 nuc = openmc.data.IncidentNeutron.from_hdf5(lib['path']) 347 # Obtain the nearest temperature 348 if strT in nuc.temperatures: 349 nucT = strT 350 else: 351 delta_T = np.array(nuc.kTs) - T * openmc.data.K_BOLTZMANN 352 closest_index = np.argmin(np.abs(delta_T)) 353 nucT = nuc.temperatures[closest_index] 354 355 # Prep S(a,b) data if needed 356 if sab_name: 357 sab = openmc.data.ThermalScattering.from_hdf5(sab_name) 358 # Obtain the nearest temperature 359 if strT in sab.temperatures: 360 sabT = strT 361 else: 362 delta_T = np.array(sab.kTs) - T * openmc.data.K_BOLTZMANN 363 closest_index = np.argmin(np.abs(delta_T)) 364 sabT = sab.temperatures[closest_index] 365 366 # Create an energy grid composed the S(a,b) and the nuclide's grid 367 grid = nuc.energy[nucT] 368 sab_Emax = 0. 369 sab_funcs = [] 370 if sab.elastic is not None: 371 elastic = sab.elastic.xs[sabT] 372 if isinstance(elastic, openmc.data.CoherentElastic): 373 grid = np.union1d(grid, elastic.bragg_edges) 374 if elastic.bragg_edges[-1] > sab_Emax: 375 sab_Emax = elastic.bragg_edges[-1] 376 elif isinstance(elastic, openmc.data.Tabulated1D): 377 grid = np.union1d(grid, elastic.x) 378 if elastic.x[-1] > sab_Emax: 379 sab_Emax = elastic.x[-1] 380 sab_funcs.append(elastic) 381 if sab.inelastic is not None: 382 inelastic = sab.inelastic.xs[sabT] 383 grid = np.union1d(grid, inelastic.x) 384 if inelastic.x[-1] > sab_Emax: 385 sab_Emax = inelastic.x[-1] 386 sab_funcs.append(inelastic) 387 energy_grid = grid 388 else: 389 energy_grid = nuc.energy[nucT] 390 391 # Parse the types 392 mts = [] 393 ops = [] 394 yields = [] 395 for line in types: 396 if line in PLOT_TYPES: 397 tmp_mts = [mtj for mti in PLOT_TYPES_MT[line] for mtj in 398 nuc.get_reaction_components(mti)] 399 mts.append(tmp_mts) 400 if line.startswith('nu'): 401 yields.append(True) 402 else: 403 yields.append(False) 404 if XI_MT in tmp_mts: 405 ops.append((np.add,) * (len(tmp_mts) - 2) + (np.multiply,)) 406 else: 407 ops.append((np.add,) * (len(tmp_mts) - 1)) 408 elif line in openmc.data.REACTION_MT: 409 mt_number = openmc.data.REACTION_MT[line] 410 cv.check_type('MT in types', mt_number, Integral) 411 cv.check_greater_than('MT in types', mt_number, 0) 412 tmp_mts = nuc.get_reaction_components(mt_number) 413 mts.append(tmp_mts) 414 ops.append((np.add,) * (len(tmp_mts) - 1)) 415 yields.append(False) 416 elif isinstance(line, int): 417 # Not a built-in type, we have to parse it ourselves 418 cv.check_type('MT in types', line, Integral) 419 cv.check_greater_than('MT in types', line, 0) 420 tmp_mts = nuc.get_reaction_components(line) 421 mts.append(tmp_mts) 422 ops.append((np.add,) * (len(tmp_mts) - 1)) 423 yields.append(False) 424 else: 425 raise TypeError("Invalid type", line) 426 427 for i, mt_set in enumerate(mts): 428 # Get the reaction xs data from the nuclide 429 funcs = [] 430 op = ops[i] 431 for mt in mt_set: 432 if mt == 2: 433 if sab_name: 434 # Then we need to do a piece-wise function of 435 # The S(a,b) and non-thermal data 436 sab_sum = openmc.data.Sum(sab_funcs) 437 pw_funcs = openmc.data.Regions1D( 438 [sab_sum, nuc[mt].xs[nucT]], 439 [sab_Emax]) 440 funcs.append(pw_funcs) 441 else: 442 funcs.append(nuc[mt].xs[nucT]) 443 elif mt in nuc: 444 if yields[i]: 445 # Get the total yield first if available. This will be 446 # used primarily for fission. 447 for prod in chain(nuc[mt].products, 448 nuc[mt].derived_products): 449 if prod.particle == 'neutron' and \ 450 prod.emission_mode == 'total': 451 func = openmc.data.Combination( 452 [nuc[mt].xs[nucT], prod.yield_], 453 [np.multiply]) 454 funcs.append(func) 455 break 456 else: 457 # Total doesn't exist so we have to create from 458 # prompt and delayed. This is used for scatter 459 # multiplication. 460 func = None 461 for prod in chain(nuc[mt].products, 462 nuc[mt].derived_products): 463 if prod.particle == 'neutron' and \ 464 prod.emission_mode != 'total': 465 if func: 466 func = openmc.data.Combination( 467 [prod.yield_, func], [np.add]) 468 else: 469 func = prod.yield_ 470 if func: 471 funcs.append(openmc.data.Combination( 472 [func, nuc[mt].xs[nucT]], [np.multiply])) 473 else: 474 # If func is still None, then there were no 475 # products. In that case, assume the yield is 476 # one as its not provided for some summed 477 # reactions like MT=4 478 funcs.append(nuc[mt].xs[nucT]) 479 else: 480 funcs.append(nuc[mt].xs[nucT]) 481 elif mt == UNITY_MT: 482 funcs.append(lambda x: 1.) 483 elif mt == XI_MT: 484 awr = nuc.atomic_weight_ratio 485 alpha = ((awr - 1.) / (awr + 1.))**2 486 xi = 1. + alpha * np.log(alpha) / (1. - alpha) 487 funcs.append(lambda x: xi) 488 else: 489 funcs.append(lambda x: 0.) 490 funcs = funcs if funcs else [lambda x: 0.] 491 xs.append(openmc.data.Combination(funcs, op)) 492 else: 493 raise ValueError(this + " not in library") 494 495 return energy_grid, xs 496 497 498def _calculate_cexs_elem_mat(this, types, temperature=294., 499 cross_sections=None, sab_name=None, 500 enrichment=None): 501 """Calculates continuous-energy cross sections of a requested type. 502 503 Parameters 504 ---------- 505 this : openmc.Material or openmc.Element 506 Object to source data from 507 types : Iterable of values of PLOT_TYPES 508 The type of cross sections to calculate 509 temperature : float, optional 510 Temperature in Kelvin to plot. If not specified, a default 511 temperature of 294K will be plotted. Note that the nearest 512 temperature in the library for each nuclide will be used as opposed 513 to using any interpolation. 514 cross_sections : str, optional 515 Location of cross_sections.xml file. Default is None. 516 sab_name : str, optional 517 Name of S(a,b) library to apply to MT=2 data when applicable. 518 enrichment : float, optional 519 Enrichment for U235 in weight percent. For example, input 4.95 for 520 4.95 weight percent enriched U. Default is None 521 (natural composition). 522 523 Returns 524 ------- 525 energy_grid : numpy.ndarray 526 Energies at which cross sections are calculated, in units of eV 527 data : numpy.ndarray 528 Cross sections calculated at the energy grid described by energy_grid 529 530 """ 531 532 if isinstance(this, openmc.Material): 533 if this.temperature is not None: 534 T = this.temperature 535 else: 536 T = temperature 537 else: 538 T = temperature 539 540 # Load the library 541 library = openmc.data.DataLibrary.from_xml(cross_sections) 542 543 if isinstance(this, openmc.Material): 544 # Expand elements in to nuclides with atomic densities 545 nuclides = this.get_nuclide_atom_densities() 546 # For ease of processing split out the nuclide and its fraction 547 nuc_fractions = {nuclide[1][0]: nuclide[1][1] 548 for nuclide in nuclides.items()} 549 # Create a dict of [nuclide name] = nuclide object to carry forward 550 # with a common nuclides format between openmc.Material and 551 # openmc.Element objects 552 nuclides = {nuclide[1][0]: nuclide[1][0] 553 for nuclide in nuclides.items()} 554 else: 555 # Expand elements in to nuclides with atomic densities 556 nuclides = this.expand(1., 'ao', enrichment=enrichment, 557 cross_sections=cross_sections) 558 # For ease of processing split out the nuclide and its fraction 559 nuc_fractions = {nuclide[0]: nuclide[1] for nuclide in nuclides} 560 # Create a dict of [nuclide name] = nuclide object to carry forward 561 # with a common nuclides format between openmc.Material and 562 # openmc.Element objects 563 nuclides = {nuclide[0]: nuclide[0] for nuclide in nuclides} 564 565 # Identify the nuclides which have S(a,b) data 566 sabs = {} 567 for nuclide in nuclides.items(): 568 sabs[nuclide[0]] = None 569 if isinstance(this, openmc.Material): 570 for sab_name in this._sab: 571 sab = openmc.data.ThermalScattering.from_hdf5( 572 library.get_by_material(sab_name, data_type='thermal')['path']) 573 for nuc in sab.nuclides: 574 sabs[nuc] = library.get_by_material(sab_name, 575 data_type='thermal')['path'] 576 else: 577 if sab_name: 578 sab = openmc.data.ThermalScattering.from_hdf5(sab_name) 579 for nuc in sab.nuclides: 580 sabs[nuc] = library.get_by_material(sab_name, 581 data_type='thermal')['path'] 582 583 # Now we can create the data sets to be plotted 584 xs = {} 585 E = [] 586 for nuclide in nuclides.items(): 587 name = nuclide[0] 588 nuc = nuclide[1] 589 sab_tab = sabs[name] 590 temp_E, temp_xs = calculate_cexs(nuc, 'nuclide', types, T, sab_tab, 591 cross_sections) 592 E.append(temp_E) 593 # Since the energy grids are different, store the cross sections as 594 # a tabulated function so they can be calculated on any grid needed. 595 xs[name] = [openmc.data.Tabulated1D(temp_E, temp_xs[line]) 596 for line in range(len(types))] 597 598 # Condense the data for every nuclide 599 # First create a union energy grid 600 energy_grid = E[0] 601 for grid in E[1:]: 602 energy_grid = np.union1d(energy_grid, grid) 603 604 # Now we can combine all the nuclidic data 605 data = np.zeros((len(types), len(energy_grid))) 606 for line in range(len(types)): 607 if types[line] == 'unity': 608 data[line, :] = 1. 609 else: 610 for nuclide in nuclides.items(): 611 name = nuclide[0] 612 data[line, :] += (nuc_fractions[name] * 613 xs[name][line](energy_grid)) 614 615 return energy_grid, data 616 617 618def calculate_mgxs(this, data_type, types, orders=None, temperature=294., 619 cross_sections=None, ce_cross_sections=None, 620 enrichment=None): 621 """Calculates multi-group cross sections of a requested type. 622 623 If the data for the nuclide or macroscopic object in the library is 624 represented as angle-dependent data then this method will return the 625 geometric average cross section over all angles. 626 627 Parameters 628 ---------- 629 this : str or openmc.Material 630 Object to source data from 631 data_type : {'nuclide', 'element', 'material', 'macroscopic'} 632 Type of object to plot 633 types : Iterable of values of PLOT_TYPES_MGXS 634 The type of cross sections to calculate 635 orders : Iterable of Integral, optional 636 The scattering order or delayed group index to use for the 637 corresponding entry in types. Defaults to the 0th order for scattering 638 and the total delayed neutron data. 639 temperature : float, optional 640 Temperature in Kelvin to plot. If not specified, a default 641 temperature of 294K will be plotted. Note that the nearest 642 temperature in the library for each nuclide will be used as opposed 643 to using any interpolation. 644 cross_sections : str, optional 645 Location of MGXS HDF5 Library file. Default is None. 646 ce_cross_sections : str, optional 647 Location of continuous-energy cross_sections.xml file. Default is None. 648 This is used only for expanding an openmc.Element object passed as this 649 enrichment : float, optional 650 Enrichment for U235 in weight percent. For example, input 4.95 for 651 4.95 weight percent enriched U. Default is None 652 (natural composition). 653 654 Returns 655 ------- 656 energy_grid : numpy.ndarray 657 Energies at which cross sections are calculated, in units of eV 658 data : numpy.ndarray 659 Cross sections calculated at the energy grid described by energy_grid 660 661 """ 662 663 # Check types 664 cv.check_type('temperature', temperature, Real) 665 if enrichment: 666 cv.check_type('enrichment', enrichment, Real) 667 cv.check_iterable_type('types', types, str) 668 669 cv.check_type("cross_sections", cross_sections, str) 670 library = openmc.MGXSLibrary.from_hdf5(cross_sections) 671 672 if data_type in ('nuclide', 'macroscopic'): 673 mgxs = _calculate_mgxs_nuc_macro(this, types, library, orders, 674 temperature) 675 elif data_type in ('element', 'material'): 676 mgxs = _calculate_mgxs_elem_mat(this, types, library, orders, 677 temperature, ce_cross_sections, 678 enrichment) 679 else: 680 raise TypeError("Invalid type") 681 682 # Convert the data to the format needed 683 data = np.zeros((len(types), 2 * library.energy_groups.num_groups)) 684 energy_grid = np.zeros(2 * library.energy_groups.num_groups) 685 for g in range(library.energy_groups.num_groups): 686 energy_grid[g * 2: g * 2 + 2] = \ 687 library.energy_groups.group_edges[g: g + 2] 688 # Ensure the energy will show on a log-axis by replacing 0s with a 689 # sufficiently small number 690 energy_grid[0] = max(energy_grid[0], _MIN_E) 691 692 for line in range(len(types)): 693 for g in range(library.energy_groups.num_groups): 694 data[line, g * 2: g * 2 + 2] = mgxs[line, g] 695 696 return energy_grid[::-1], data 697 698 699def _calculate_mgxs_nuc_macro(this, types, library, orders=None, 700 temperature=294.): 701 """Determines the multi-group cross sections of a nuclide or macroscopic 702 object. 703 704 If the data for the nuclide or macroscopic object in the library is 705 represented as angle-dependent data then this method will return the 706 geometric average cross section over all angles. 707 708 Parameters 709 ---------- 710 this : openmc.Nuclide or openmc.Macroscopic 711 Object to source data from 712 types : Iterable of str 713 The type of cross sections to calculate; values can either be those 714 in openmc.PLOT_TYPES_MGXS 715 library : openmc.MGXSLibrary 716 MGXS Library containing the data of interest 717 orders : Iterable of Integral, optional 718 The scattering order or delayed group index to use for the 719 corresponding entry in types. Defaults to the 0th order for scattering 720 and the total delayed neutron data. 721 temperature : float, optional 722 Temperature in Kelvin to plot. If not specified, a default 723 temperature of 294K will be plotted. Note that the nearest 724 temperature in the library for each nuclide will be used as opposed 725 to using any interpolation. 726 727 Returns 728 ------- 729 data : numpy.ndarray 730 Cross sections calculated at the energy grid described by energy_grid 731 732 """ 733 734 # Check the parameters and grab order/delayed groups 735 if orders: 736 cv.check_iterable_type('orders', orders, Integral, 737 min_depth=len(types), max_depth=len(types)) 738 else: 739 orders = [None] * len(types) 740 for i, line in enumerate(types): 741 cv.check_type("line", line, str) 742 cv.check_value("line", line, PLOT_TYPES_MGXS) 743 if orders[i]: 744 cv.check_greater_than("order value", orders[i], 0, equality=True) 745 746 xsdata = library.get_by_name(this) 747 748 if xsdata is not None: 749 # Obtain the nearest temperature 750 t = np.abs(xsdata.temperatures - temperature).argmin() 751 752 # Get the data 753 data = np.zeros((len(types), library.energy_groups.num_groups)) 754 for i, line in enumerate(types): 755 if 'fission' in line and not xsdata.fissionable: 756 continue 757 elif line == 'unity': 758 data[i, :] = 1. 759 else: 760 # Now we have to get the cross section data and properly 761 # treat it depending on the requested type. 762 # First get the data in a generic fashion 763 temp_data = getattr(xsdata, _PLOT_MGXS_ATTR[line])[t] 764 shape = temp_data.shape[:] 765 # If we have angular data, then want the geometric 766 # average over all provided angles. Since the angles are 767 # equi-distant, un-weighted averaging will suffice 768 if xsdata.representation == 'angle': 769 temp_data = np.mean(temp_data, axis=(0, 1)) 770 771 # Now we can look at the shape of the data to identify how 772 # it should be modified to produce an array of values 773 # with groups. 774 if shape in (xsdata.xs_shapes["[G']"], 775 xsdata.xs_shapes["[G]"]): 776 # Then the data is already an array vs groups so copy 777 # and move along 778 data[i, :] = temp_data 779 elif shape == xsdata.xs_shapes["[G][G']"]: 780 # Sum the data over outgoing groups to create our array vs 781 # groups 782 data[i, :] = np.sum(temp_data, axis=1) 783 elif shape == xsdata.xs_shapes["[DG]"]: 784 # Then we have a constant vs groups with a value for each 785 # delayed group. The user-provided value of orders tells us 786 # which delayed group we want. If none are provided, then 787 # we sum all the delayed groups together. 788 if orders[i]: 789 if orders[i] < len(shape[0]): 790 data[i, :] = temp_data[orders[i]] 791 else: 792 data[i, :] = np.sum(temp_data[:]) 793 elif shape in (xsdata.xs_shapes["[DG][G']"], 794 xsdata.xs_shapes["[DG][G]"]): 795 # Then we have an array vs groups with values for each 796 # delayed group. The user-provided value of orders tells us 797 # which delayed group we want. If none are provided, then 798 # we sum all the delayed groups together. 799 if orders[i]: 800 if orders[i] < len(shape[0]): 801 data[i, :] = temp_data[orders[i], :] 802 else: 803 data[i, :] = np.sum(temp_data[:, :], axis=0) 804 elif shape == xsdata.xs_shapes["[DG][G][G']"]: 805 # Then we have a delayed group matrix. We will first 806 # remove the outgoing group dependency 807 temp_data = np.sum(temp_data, axis=-1) 808 # And then proceed in exactly the same manner as the 809 # "[DG][G']" or "[DG][G]" shapes in the previous block. 810 if orders[i]: 811 if orders[i] < len(shape[0]): 812 data[i, :] = temp_data[orders[i], :] 813 else: 814 data[i, :] = np.sum(temp_data[:, :], axis=0) 815 elif shape == xsdata.xs_shapes["[G][G'][Order]"]: 816 # This is a scattering matrix with angular data 817 # First remove the outgoing group dependence 818 temp_data = np.sum(temp_data, axis=1) 819 # The user either provided a specific order or we resort 820 # to the default 0th order 821 if orders[i]: 822 order = orders[i] 823 else: 824 order = 0 825 # If the order is available, store the data for that order 826 # if it is not available, then the expansion coefficient 827 # is zero and thus we already have the correct value. 828 if order < shape[1]: 829 data[i, :] = temp_data[:, order] 830 else: 831 raise ValueError("{} not present in provided MGXS " 832 "library".format(this)) 833 834 return data 835 836 837def _calculate_mgxs_elem_mat(this, types, library, orders=None, 838 temperature=294., ce_cross_sections=None, 839 enrichment=None): 840 """Determines the multi-group cross sections of an element or material 841 object. 842 843 If the data for the nuclide or macroscopic object in the library is 844 represented as angle-dependent data then this method will return the 845 geometric average cross section over all angles. 846 847 Parameters 848 ---------- 849 this : openmc.Element or openmc.Material 850 Object to source data from 851 types : Iterable of str 852 The type of cross sections to calculate; values can either be those 853 in openmc.PLOT_TYPES_MGXS 854 library : openmc.MGXSLibrary 855 MGXS Library containing the data of interest 856 orders : Iterable of Integral, optional 857 The scattering order or delayed group index to use for the 858 corresponding entry in types. Defaults to the 0th order for scattering 859 and the total delayed neutron data. 860 temperature : float, optional 861 Temperature in Kelvin to plot. If not specified, a default 862 temperature of 294K will be plotted. Note that the nearest 863 temperature in the library for each nuclide will be used as opposed 864 to using any interpolation. 865 ce_cross_sections : str, optional 866 Location of continuous-energy cross_sections.xml file. Default is None. 867 This is used only for expanding the elements 868 enrichment : float, optional 869 Enrichment for U235 in weight percent. For example, input 4.95 for 870 4.95 weight percent enriched U. Default is None 871 (natural composition). 872 873 Returns 874 ------- 875 data : numpy.ndarray 876 Cross sections calculated at the energy grid described by energy_grid 877 878 """ 879 880 if isinstance(this, openmc.Material): 881 if this.temperature is not None: 882 T = this.temperature 883 else: 884 T = temperature 885 886 # Check to see if we have nuclides/elements or a macrocopic object 887 if this._macroscopic is not None: 888 # We have macroscopics 889 nuclides = {this._macroscopic: (this._macroscopic, this.density)} 890 else: 891 # Expand elements in to nuclides with atomic densities 892 nuclides = this.get_nuclide_atom_densities() 893 894 # For ease of processing split out nuc and nuc_density 895 nuc_fraction = [nuclide[1][1] for nuclide in nuclides.items()] 896 else: 897 T = temperature 898 # Expand elements in to nuclides with atomic densities 899 nuclides = this.expand(100., 'ao', enrichment=enrichment, 900 cross_sections=ce_cross_sections) 901 902 # For ease of processing split out nuc and nuc_fractions 903 nuc_fraction = [nuclide[1] for nuclide in nuclides] 904 905 nuc_data = [] 906 for nuclide in nuclides.items(): 907 nuc_data.append(_calculate_mgxs_nuc_macro(nuclide[0], types, library, 908 orders, T)) 909 910 # Combine across the nuclides 911 data = np.zeros((len(types), library.energy_groups.num_groups)) 912 for line in range(len(types)): 913 if types[line] == 'unity': 914 data[line, :] = 1. 915 else: 916 for n in range(len(nuclides)): 917 data[line, :] += nuc_fraction[n] * nuc_data[n][line, :] 918 919 return data 920